NI MATRIXx Xmath Control Design Module ™ Xmath Control Design Module April 2007 370753C-01 TM
Support Worldwide Technical Support and Product Information ni.
Important Information Warranty TThe media on which you receive National Instruments software are warranted not to fail to execute programming instructions, due to defects in materials and workmanship, for a period of 90 days from date of shipment, as evidenced by receipts or other documentation. National Instruments will, at its option, repair or replace software media that do not execute programming instructions if National Instruments receives notice of such defects during the warranty period.
Conventions The following conventions are used in this manual: <> Angle brackets that contain numbers separated by an ellipsis represent a range of values associated with a bit or signal name—for example, DIO<3..0>. [] Square brackets enclose optional items—for example, [response]. » The » symbol leads you through nested menu items and dialog box options to a final action.
Contents Chapter 1 Introduction Using This Manual.........................................................................................................1-1 Document Organization...................................................................................1-1 Bibliographic References ................................................................................1-2 Commonly Used Nomenclature ......................................................................1-2 Related Publications ...............
Contents Chapter 3 Building System Connections Linear System Interconnection Operators ..................................................................... 3-1 Linear System Interconnection Functions ..................................................................... 3-4 afeedback( )..................................................................................................... 3-4 Algorithm.......................................................................................... 3-5 append( ) ....
Contents Bode Frequency Analysis ..............................................................................................5-7 bode( )..............................................................................................................5-10 margin( ) ..........................................................................................................5-12 nichols( )..........................................................................................................
Contents Appendix A Technical References Appendix B Technical Support and Professional Services Index Xmath Control Design Module viii ni.
1 Introduction The Control Design Module (CDM) is a complete library of classical and modern control design functions that provides a flexible, intuitive control design framework. This chapter starts with an outline of the manual and some user notes. It also contains a tutorial that presents several problems and uses a variety of approaches to obtain solutions. The tutorial is designed to familiarize you with many of the functions in this module.
Chapter 1 Introduction particular system properties or to change the format of a system. These topics include continuous/discrete system conversion, as well as finding equivalent transfer function state-space representations. • Chapter 3, Building System Connections, details Xmath functions that perform different types of linear system interconnections. It also discusses a number of simpler connections that have been implemented as overloaded operators on system objects.
Chapter 1 Introduction • H(s) is used to denote the frequency response, over some range of frequencies of a system where s is the Laplace variable. H(q) is used to indicate that the system can be continuous or discrete. • A single apostrophe following a matrix variable, for example, x', denotes the transpose of that variable. An asterisk following a matrix variable (for example, A*) indicates the complex conjugate, or Hermitian, transpose of that variable.
Chapter 1 Introduction Control Design Tutorial This tutorial illustrates the use of functions and commands provided in Xmath and the Xmath Control Design Module to solve control problems. The emphasis of the tutorial is on using a number of different approaches, not on any one “correct” way to solve a problem. It demonstrates the flexibility of Xmath’s tools and scripting language to customize your analysis in a way that is as straightforward and mathematically intuitive as possible.
Chapter 1 Introduction ssys (a state space system) = A -0.4 1 -1.4 0 -0.01 0 0 9.8 -0.02 B 6.3 0 9.8 C 0 0 1 D 0 X0 0 0 0 State Names ----------Pitch Rate Pitch Angle Horizontal v Input Names ----------Rotor Angle Output Names -----------Horizontal v System is continuous Use check( ) to convert the model to transfer-function form: [,Gs] = check(ssys,{tf,convert}) Gs (a transfer function) = 2 9.8(s - 0.5s + 6.3) ----------------------------------------2 (s + 0.656513)(s - 0.236513s + 0.
Chapter 1 Introduction Input Names ----------Rotor Angle Output Names -----------Horizontal v System is continuous The system has poles and zeros in the right half of the complex plane and therefore is open-loop unstable. Checking the pole and zero locations confirms this: ol_poles=poles(ssys) ol_poles (a column vector) = 0.118256 - 0.367817 j 0.118256 + 0.367817 j -0.656513 ol_zeros=zeros(ssys) ol_zeros (a column vector) = 0.25 + 2.4975 j 0.25 - 2.
Chapter 1 Introduction One approach to stabilizing this system is to try to cancel the pole at –0.656513 by adding a compensator, K1(s), with a zero at –0.656513. It is important to understand that this is primarily an academic exercise. Accurate pole-zero cancellations are impracticable in the real world, and the mode corresponding to that pole still exists internally.
Chapter 1 Introduction Figure 1-2. Locus of all Open-Loop and Closed-Loop Roots of Gs If you cannot move the slider so that the gain is exactly 2, click the box to the right of the slider and enter 2. To close the interactive root locus dialog box, select File»Exit. Xmath Control Design Module 1-8 ni.
Chapter 1 Introduction Close the loop using the single-input syntax of feedback( ), which implements direct unity-gain negative feedback, and obtain the system’s step response using step( ): Kc = 2; cl_syscomp1 = feedback(Kc*K1s*K2s*Gs); v = step(cl_syscomp1, 0:.2:25); plot(v,{xlab="Time", ylab="Horizontal Velocity"}) The resulting plot is shown in Figure 1-3. This result is not desirable.
Chapter 1 Introduction U(s) + – Y(s) G(s) Kc1K1(s) Kc2K2(s) Figure 1-4. Block Diagram of the Closed-Loop Controller This is a block diagram of the closed-loop controller with compensator Kc1K1(s) in the feedforward path and Kc2K2(s) in the feedback path. This time, instead of having all your gain Kc in the forward path of the closed-loop system, the system gain is split between the two compensators.
Chapter 1 Introduction Because cl_syscomp2 contains an internal pole-zero cancellation, you can rewrite it in minimal form and then check the closed-loop pole and zero locations: cl_syscomp2m = minimal(cl_syscomp2); The system has 1 uncontrollable state cl_poles = poles(cl_syscomp2m) cl_poles (a column vector) = -0.166518 -1.0325 + 1.16109 j -1.0325 - 1.16109 j -37.132 cl_zeros = zeros(cl_syscomp2m) cl_zeros (a column vector) = 0.25 - 2.4975 j 0.25 + 2.
Chapter 1 Introduction Figure 1-5. Helicopter Velocity Tracking Step Input at the Rotor You also can look at the gain and phase margins of the system. H = bode(cl_syscomp2m, {npts = 200, !wrap}); [gm,pm] = margin(H) There are no 0 dB gain crossings. gm (a pdm) = domain | ---------+---------0.250101 | 26.1694 ---------+---------pm is null The bode plot of the closed-loop system is shown in Figure 1-6. Xmath Control Design Module 1-12 ni.
Chapter 1 Introduction Figure 1-6. Closed-Loop System Bode Plot The domain of the gain and phase margin PDMs indicates the frequency (in hertz) at which the margin occurs. So the gain can be increased by about 26.1 dB before the system becomes unstable.
Chapter 1 Introduction You can verify that your system is controllable, then define the closed-loop poles you want and use poleplace( ) to find the feedback gains required given the system A and B matrices. [,,nuc] = controllable(Gs) nuc (a scalar) = 0 Because the number of uncontrollable states is zero, Gs is controllable. This means that you can use feedback through appropriately-sized gains to position the system’s closed-loop poles anywhere you want.
Chapter 1 Introduction Specify the observer poles at [–3 + 3j, –4] and call poleplace( ) again: op = [-3+3*jay, -4]; L = poleplace(A',C',op) L (a row vector) = 5.46645 4.67623 9.58 You connect the controller to the observer using lqgcomp( ). L needs to be a column vector, so you transpose it. sys_obc = lqgcomp(ssys, Kfsb, L'); You can use names( ) to modify the names of the state estimates to be more descriptive.
Chapter 1 Introduction You can choose to scale the system output here for zero steady-state error in the step response. This is accomplished in an intuitive manner, dividing the system sys_cl by the desired scaling factor. sys_cl = sys_cl/51.76; v_obc = step(sys_cl, 0:.1:10); plot (v_obc, {xlab = "Time", ylab = "Magnitude"}) In Figure 1-9 the step response shows zero-steady-state error, little overshoot, and a response time of less than seven seconds. Figure 1-9.
Chapter 1 Introduction When you create the estimator system sys_est, you use the original A matrix for the state-update equation, but you provide a zero external input (a B matrix of zero). The output matrix is an identity matrix passing back the three real state values and the three estimated state values as output, again with no external input values affecting the output. Here you use the optional system( ) keyword X0 to set the real state values to [1,2,3] and the estimated state values to [–1,–2,–3].
Chapter 1 Introduction Figure 1-10. Multiple Plots Showing Time Needed for States to be Correctly Tracked by Estimator, Given Incorrect Initial Values Helicopter Hover Problem: Discrete Formulation Discrete-time control systems are most frequently designed in one of two ways: either directly implemented in the discrete domain, or first solved as continuous problems—often deriving directly from differential equations of motion—and then discretized.
Chapter 1 Introduction You can use the default exponential discretization method with dt = 0.01 and compare frequency responses between the original system and the discretized system: ssysd = discretize(ssys, 0.01); f = freq(ssys,logspace(.001,10,200)); fd = freq(ssysd,logspace(.001,10,200)); In the following statements you compute the gain and phase of both systems and then plot them.
Chapter 1 Introduction Figure 1-12. Step Response of a Discrete System Using Discretized Observer-Based Controller As you discretize the compensator, form the closed-loop, scaled system, and simulate its response to a step input, you must ensure that the sampling interval is the same (dt = 0.01). sys_obcd = discretize(sys_obc, 0.01); sys_cld = feedback(ssysd,sys_obcd)/51.76; v_cld = step(sys_cld, 0:0.
Chapter 1 Introduction The linearized state-space equations, including the actuator and sensor dynamics, are as follows: · θ 0 0 ·x 0 0 = ·· 15.54 – 10.93 θ – 5.31 0 x·· 1 0 0 1 0 0 0 – 16.24 θ 0 x + 0 u · 0 θ ·x 1.96 θ x y = 57.29 0 0 0 0 29.9 0 0 θ· x· θ is the angle (in radians) the wedge makes with the vertical axis, x is the position of the sliding mass, and u is the control input voltage. The outputs are scaled to give the measured angle in degrees and the measured position in meters.
Chapter 1 Introduction Because this system is open-loop unstable and has fairly fast poles in both halves of the s-plane, you want to make sure it can bring the effect of an external disturbance (such as a sharp push to the cart) to zero as quickly as possible.
Chapter 1 Introduction Figure 1-13. Response of Full-State Feedback Controller to a Unit Step Disturbance Having established your regulator design, you build the estimator and simulate performance of the closed-loop system feeding back state estimates. You select the weights for the estimator based on the assumption that the state noise intensities corresponding to the wedge angle are smaller than those corresponding to the wedge position.
Chapter 1 Introduction Figure 1-14. Response of Observer-Based Controller to a Unit Step Disturbance Xmath Control Design Module 1-24 ni.
2 Linear System Representation Xmath provides a structure for system representation called a system object. This object includes system parameters in a data structure designed to reflect the way these systems are analyzed mathematically. Operations on these systems are likewise defined using operators that mirror as closely as possible the notation control engineers use.
Chapter 2 Linear System Representation Transfer Function System Models One way of representing continuous-time finite-dimensional linear time-invariant systems is with the transfer function: num ( s ) H ( s ) = -----------------den ( s ) with num(s) and den(s) being polynomials in s. They can be specified either by their roots or their coefficients. Transfer functions are defined using the Laplace transform operators for continuous time and the forward shift operator z for discrete time.
Chapter 2 Linear System Representation form the same transfer function as that derived in the preceding transfer function equation using known pole, zero, and gain values: 2 ( s – 0.5 ) H ( s ) = --------------------------------(s + 2)(s + 4) (2-4) The systems represented in Equations 2-3 and 2-4 can be represented using Xmath’s system objects, as shown in Example 2-1. The Xmath transfer function system object currently can be used to represent single-input, single-output systems only.
Chapter 2 Linear System Representation be used so long as a consistent choice of variable is used for both numerator and denominator polynomials. The transfer function in pole-zero-gain form from the preceding equation can be similarly implemented using the polynomial( ) function to specify the numerator and denominator by their roots. The / operator also can be used to create systems in transfer function form, as an alternative to using system( ). Note H4 = 2*polynomial(0.
Chapter 2 Linear System Representation State-Space System Models State-space models comprise the second category of linear system representations in Xmath. In state-space form, first-order differential (continuous-time) and difference (discrete-time) equations are represented as a set of state and output updates. The states are represented by a vector x; u and y are vectors with as many elements as there are inputs and outputs, respectively.
Chapter 2 Linear System Representation Again, you create the system using the system( ) function. This time you use the optional dt keyword to indicate that this system is discrete. A = [0,1;-0.75,0]; B = [1,0]'; C = [0,1]; D = 0; sys4 = system(A,B,C,D, {dt = 0.5}) sys4 (a state space system) = A 0 1 -0.75 0 B 1 0 C 0 1 D 0 X0 0 0 System is discrete, sampling at 0.5 seconds.
Chapter 2 Linear System Representation The system( ) function can create both the transfer-function and state-space forms of the system object. It requires four compatibly-sized matrices to create a state-space system, or a pair of polynomials to create a transfer function. You can use optional keywords to store additional information about your system. Assigning dt to a positive scalar value indicates that the system is discrete, with a sampling period equal to that value.
Chapter 2 Example 2-3 Linear System Representation Using system( ) to Change the Attributes of an Existing System sys4=system([0,1;-0.75,0],[1,0]',[0,1],[], {dt=0.5}); sys4 = system(sys4, {inputNames = "Current", outputNames = "Velocity", stateNames = ["Torque","Angle"]}) sys4 (a state space system) = A 0 1 -0.75 0 B 1 0 C 0 1 D 0 X0 0 0 State Names ----------Torque Angle Input Names ----------Current Output Names -----------Velocity System is discrete, sampling at 0.5 seconds.
Chapter 2 Linear System Representation done internally to return A, B, C, and D, though the format of the variable Sys itself remains unchanged. The transfer function must be proper. Using the systems defined in Examples 2-1 and 2-2, Example 2-4 illustrates the use of abcd( ). Example 2-4 Using abcd( ) to Extract the State-Space Matrices H3=makepoly([2,-1],"s")/makepoly([1,6,8],"s"); sys4=system([0,1;-0.75,0],[1,0]',[0,1],0, {dt=0.5}); You can extract the state-space matrices from each.
Chapter 2 Linear System Representation numden( ) [num,den] = numden(Sys) The numden( ) function returns the numerator and denominator polynomials comprising a single-input, single-output system in transfer function form. If the system is in state-space form, an internal conversion is performed to find the transfer function equivalent, but the format of the system variable itself remains unchanged. State-space systems used in conjunction with numden( ) must be single-input, single-output.
Chapter 2 Linear System Representation these polynomials into a transfer-function and uses period( ) to set the sampling interval to match that of sys4. Example 2-6 Using period( ) to Extract the Sampling Period [num,den]=numden(sys4); H4 = system(num,den,{dt = period(sys4)}) H4 (a transfer function) = -0.75 ----------(z2 + 0.75) System is discrete, sampling at 0.5 seconds.
Chapter 2 Linear System Representation outputNames = "Velocity", stateNames = ["Torque","Angle"]}); [,,stateNames] = names(sys5)? statenames (a row vector of strings)=Torque Angle Size and Indexing of Dynamic Systems The size of a system object is defined by how many outputs, inputs, and (in the case of a state-space system) states it has. You can use the size( ) function to find these dimensions.
Chapter 2 Linear System Representation check(sys, {stable}) ans (a scalar) = 0 check(sys, {discrete, ss}) ans (a scalar) = 1 [, tfsys] = check(sys, {tf, convert}) tfsys (a transfer function) = (z + 0.26) --------------------(z + 0.26)(z - 1.875) initial delay outputs 0 0 System is discrete, sampling at 0.001 seconds. Discretizing a System Many systems where behavior derives from physical equations of motion can be modeled most naturally as continuous processes, using differential equations.
Chapter 2 Linear System Representation discretization methods used based on the specification of each keyword are discussed in the following sections. Numerical Integration Methods: forward, backward, tustins Xmath provides three methods of numerical integration of a differential transfer function: the forward and backward rectangular rules, and Tustin’s rule (also called the bilinear or trapezoidal transform).
Chapter 2 Linear System Representation Table 2-2.
Chapter 2 Linear System Representation The exponential keyword assumes that the response value between samples is constant and can, therefore, be represented by a zero-order hold polynomial. When exponential is specified, the continuous-time step response is discretized using the Z-transform, then the result is divided by the Z-transform of a step z/(z – 1) to produce the desired transfer function.
Chapter 2 Linear System Representation gain_z = 20*log10(abs(freq(Hd_z,F))); gain_e = 20*log10(abs(freq(Hd_e,F))); and plot it (as shown in Figure 2-1). plot ([gainc,gain_f,gain_b,gain_t,gain_z,gain_e], {legend = ["Continuous", "Forward", "Backward", "Tustins", "Pole Zero", "Exponential"], x_log, xlab="Frequency (Hz)",ylab="Magnitude (dB)"}) Figure 2-1.
Chapter 2 Linear System Representation Many of the discretization techniques discussed in the Hold Equivalence Methods: exponential and firstorder section can be easily reversed to obtain a continuous equivalent to a discrete system. The makecontinuous( ) function implements these reverse algorithms based on the keyword you specify as shown in Example 2-10. Although makecontinuous( ) accepts an input system in any form, it returns the continuous-time system as a state-space system.
Chapter 2 Linear System Representation Now convert back to the continuous form: Hc = makecontinuous(Hd_f, {forward}); [num,den] = numden(Hc) num (a polynomial) = (s + 0.36) den (a polynomial) = 2 (s + 0.999998)(s + 1)(s + 0.79s + 0.16) Although makecontinuous( ) restores the continuous-time poles and zeros, it cannot match gains precisely.
3 Building System Connections Large system models are frequently built by connecting smaller models together. You can perform different types of linear system interconnections using the Xmath functions discussed in this chapter. MathScript allows operators (*,+, and so on) to be overloaded—given different behaviors when used with different objects.
Chapter 3 Building System Connections Table 3-1. Summary of Interconnection Operators (Continued) Diagram Description Cascade connection of Sys1 and Sys2 where the output of Sys is y2 and the input is u1. Sys = Sys2 * Sys1 u1 u y1 u2 y2 Sys2 Sys1 y Cascade connection of Sys1 and inverted Sys2 where Sys = Sys1 * inv(Sys2), u = u2, and y = y1. Sys = Sys1/Sys2 u u2 y2 u1 y1 Sys1 inv(Sys2) y Cascade connection of inverted Sys1 and Sys2 where Sys = inv(Sys1) * Sys2, u = u2, and y = y1.
Chapter 3 Building System Connections Table 3-1. Summary of Interconnection Operators (Continued) Diagram Description Sys = adj[Sys1] If Sys1 is in state-space form and comprises the matrices (A1,B1,C1,D1), Sys is the adjoint system and comprises (–A1',C1',B1',D1'). If Sys1 is a transfer function, it is converted internally to state-space form. p1/p2 Alternate method to create a system, where p1 and p2 are the numerator and denominator polynomials, respectively; does not allow the use of keywords.
Chapter 3 Building System Connections Linear System Interconnection Functions afeedback( ), append( ), connect( ), and feedback( ) connect dynamic systems in state–space or transfer–function form to produce a larger system in state-space form. The following restrictions apply to all of these functions: • Both systems must have the same sample rate. • Improper dynamic systems (systems with more zeros than poles) are not allowed.
Chapter 3 Example 3-1 Building System Connections • By default, feedback is defined to be negative. • The number of outputs from the first system must equal the number of inputs to the second system. • The number of outputs from the second system must equal the number of inputs in the first. • Both systems must have the same sample rate. • Improper dynamic systems (systems with more zeros than poles) are not allowed.
Chapter 3 Building System Connections if the condition estimate for either matrix is less than eps. For more information on this condition estimate, refer to the MATRIXx Help for the Xmath function rcond( ). Using Sys to denote the state-space representation for the complete feedback loop, you can express its component matrices A, B, C, and D as combinations of the component matrices you obtained from Sys1 and Sys2. The full matrices used with two input systems are shown in the next example.
Chapter 3 Building System Connections is created by appending the inputs, outputs, and states of Sys1 and Sys2. A larger number of systems can be appended by appending two at a time. • Both systems must have the same sample rate. • Improper dynamic systems—systems with more zeros than poles—are not allowed, because Sys is represented in state-space form. The output is a dynamic system in block form as shown in Figure 3-2. u1 Sys1 y1 u2 Sys2 y2 Sys Figure 3-2.
Chapter 3 Building System Connections A = A1 0 C = C1 0 0 A2 0 C2 B = B1 0 D = D1 0 0 B2 0 D2 append( ) performs a check to make sure both Sys1 and Sys2 have the same sample rate, and adopts this rate for the appended system. Any initial conditions on the states are also appended columnwise. connect( ) Sys = connect(Sys1,{K,M,N}) The connect( ) function performs a general interconnection around a system.
Chapter 3 Building System Connections • By default, feedback is defined to be positive. To enforce negative feedback, specify connect(Sys,-K). • A “selection matrix” has a single 1 in each row; the rest of the row contains zeros. This is useful for indicating the subset of system inputs and outputs to be used. In many cases, however, it is simpler to extract desired inputs and outputs through indexing. • Both systems must have the same sample rate.
Chapter 3 Building System Connections 2 C -1.5 1.5 D 0 X0 0 0 Algorithm For the feedback system shown in Example 3-3, you can write the following system equations: x· = A 1 x + B 1 u 1 y1 = C1 x + D1 u1 Combining these equations with the equation for the positive feedback input term: u 1 = Ky 1 + Mu and multiplying by the input and output gains M and N, you obtain the following state-space equations describing the entire system between input u and output y.
Chapter 3 Building System Connections feedback( ) Sys = feedback(Sys1,{Sys2}) The feedback( ) function connects two dynamic systems together in a feedback loop as shown in Figure 3-4. • By default, feedback is defined to be negative. • Both systems must have the same sample rate. • Improper dynamic systems (systems with more zeros than poles) are not allowed. • When only one system is specified, it must be square (it must have an equal number of inputs and outputs).
Chapter 3 Building System Connections -1 1 D 0 X0 0 0 State Names ----------Input Names ----------Input 1 Output Names -----------Output 1 System is continuous Algorithm The system used for the feedback loop, Sys2, is optional. If it is not specified, a default state-space system is used with NULL matrices for A2, B2, and C2 and an identity matrix for D2 so that unity gain feedback is implemented.
Chapter 3 Building System Connections The single system resulting from the feedback combination of Sys1 and Sys2 has u1 as its input, y1 as its output, and a state vector consisting of the appended states of Sys1 and Sys2. Using these five equations to find the state-space dynamics of the complete system results in the overall system description.
4 System Analysis This chapter discusses time-domain solutions of the equations underlying transfer functions and state-space system models, and what these solutions tell us about the stability of the system. Xmath provides a number of functions for performing this system analysis and computing the time-domain system response to both general and specific “standard” inputs.
Chapter 4 System Analysis The time-response of discrete systems is found directly as a summation of the information from preceding time points in the state and input histories.
Chapter 4 System Analysis and define the zeros of S(λ) as any values of λ for which the system matrix drops rank. For single-input single-output systems this is equivalent to the polynomial zeros of the transfer-function numerator. This definition is somewhat more complex for MIMO systems. In terms of the dynamic response associated with the poles and zeros of a system, a pole is said to be stable if the response it contributes decays over time.
Chapter 4 System Analysis directly from the roots of the transfer function numerator. If Sys is in state-space form, the definition of its zeros arises from the system matrix, S ( λ ) = λI – A B –C D (4-3) and its MIMO transfer function: –1 R ( λ ) = C ( λI – A ) B + D (4-4) Defining n as the number of states in the system, p as the number of outputs, and m as the number of inputs, the normal rank of S(λ) is n + min(m,p). If the rank of S(λ) equals the normal rank, the system is nondegenerate.
Chapter 4 System Analysis ans (a column vector) = -0.98875 + 2.4773 j -0.98875 - 2.4773 j k (a scalar) = 1.34P Algorithm The algorithm used for state-space zero computation creates a reduced-order S(λ), using Householder reflections to do the necessary orthogonal row and column compressions of the state-space matrices. The eigenvalues of this reduced matrix are then found using QZ. This method handles the degenerate case and systems with zeros at infinity [EmV82].
Chapter 4 Example 4-3 System Analysis Dynamic Response through Partial Fraction Expansion To illustrate how you can examine the stability and dynamic response of a system using Xmath, start with the open-loop transfer function system ( s + 0.5 ) G ( s ) = ---------------------------------------2 s ( s + 2 ) ( s + 10 ) You close a unity feedback loop around this system, as shown in Figure 4-1. + V(s) u(s) Y(s) G(s) – Gcl(s) Figure 4-1.
Chapter 4 System Analysis 0 0 Input Names ----------Input 1 Output Names -----------Output 1 System is continuous You can examine the stability of Gcl(s) by representing it as a sum of partial fractions, using the residue( ) function. residue(syscl) ans (a pdm) = Poles | -------------------------+-----------------------0.0177496 - 0.158936 j | Order 1 0.0180045 + ... -0.0177496 + 0.158936 j | Order 1 0.0180045 -... -1.95266 | Order 1 -0.0478224 -10.0118 | Order 1 0.
Chapter 4 System Analysis Figure 4-2. Transient Response of the Closed-Loop System as a Function of Time You also can calculate the impulse response directly with t = [0 : 0.1 : 350 ]; hi = impulse(syscl,t); plot(hi, { xlab = "Time (sec)", ylab = "Transient Response"}) Calculating the impulse response gives you the transient response shown in Figure 4-2.
Chapter 4 System Analysis and orders for which the residue(s) should be found. If a user-specified value for pls is not actually a pole of the system or if the order requested is greater than the multiplicity of the pole, the corresponding residue is returned as zero. C contains the value of the constant term. Example 4-4 uses the transfer function from Example 2-10, Verifying a Discretization Using makecontinuous( ). Example 4-4 Calculating the Residues of a System G= system(0.5*polynomial([-0.
Chapter 4 System Analysis (s + 1) (s + 0.79s + 0.16) initial integrator outputs 0 0 0 0 Input Names ----------Input 1 Output Names -----------Output 1 System is continuous Note G2 matches the system G where residues were computed in Example 4-4.
Chapter 4 System Analysis Often it is desirable to run several simulations with different inputs. In this case, you can define a PDM whose columns contain the input vectors for the different simulations. Then u will be ny × q × Nsamp, where q is the number of different simulations to be run. The resulting y will be ny × q × Nsamp, with each column of the PDM corresponding to a different simulation.
Chapter 4 System Analysis Figure 4-3. System Time Response to a Series of Step Signals The (system)*(PDM) construct for performing time-domain simulation is used analogously no matter how many inputs the system has. For a multi-input, multi-output system, the number of rows of the input PDM must match the number of inputs of the system. For an example of general time-domain simulation for a MIMO system, refer to Example 4-6.
Chapter 4 System Analysis Impulse Response of a System An impulse input to a system is defined somewhat differently depending on whether the system is discrete or continuous. For a continuous-time system, an impulse is a unit-area signal of infinite amplitude and infinitely small duration occurring at time t = 0, and having value zero at all other times.
Chapter 4 System Analysis A continuous system and its discrete-time equivalent (computed using the impulse-invariant z-transform) have impulse responses differing only by a factor of 1/dt. Note impulse( ) computes the impulse response by using the B matrix from the system’s state-space representation as the initial conditions. A system with ni inputs has ni initial conditions, each of which is set up as a column of the B matrix.
Chapter 4 System Analysis Figure 4-4. 15-Second Impulse Response deftimerange( ) tvec = deftimerange(Sys) deftimerange( ) computes a regular time vector (in units of seconds) that can be used in time-domain simulations to observe the effects of all or most of the input system’s dynamics, as indicated by pole and zero location. Within deftimerange( ), the poles of the system are obtained using poles( ).
Chapter 4 System Analysis The maximum time, Tmax, is computed as follows, with vP denoting the vector of scaled poles and dt the period: Tmax=abs(log(.05)/... max(real(vP(find(real(vP)<>0))))) If Tmax == null # all poles purely imaginary Tmax=100*dt endIf Tmax=max(Tmax,10*dt) tvec=0:dt:Tmax Though deftimerange( ) calls minimal( ) to remove any pole-zero cancellations, it does not consider the location of the system zeros in computing the time vector.
Chapter 4 System Analysis The simulation performed in initial( ) uses an input of zero for each point in the time vector. The output Y is a PDM where domain is the time vector. By varying the initial values of the states individually, you can determine which is the most sensitive. For an example using initial( ) to determine the sensitivity of the states, refer to Example 4-8. Example 4-8 Using initial( ) to Determine the Sensitivity of the States Sys = system([-2.3,0.01,5.1;0,-0.35,-2;0,2,-.35], [1,.
Chapter 4 System Analysis Step Response The response of a system to a unit step input is one of the most commonly used measures of how well a given control system’s output tracks the system input. A unit step is a time signal which is zero up until the beginning of the time period of interest, and one thereafter. This indicator is popular because it is easy to compute and interpret. It also is mathematically possible to calculate the response to any input if the response to a unit step is known.
Chapter 4 System Analysis From Figure 4-6 you see that the delay time (td) is about 0.5 seconds, the rise time (tr) is 0.8 seconds, the peak time (tp) is 1.6 seconds, the settling time (ts) is about 5.5 seconds, and the maximum overshoot (Mp) is about 24%. Mp ts tp tr td Figure 4-6. 15-Second Step Response, Showing Performance Measures You can compute these values from the 151-point step response data vector Y and substantiate your estimates.
Chapter 4 System Analysis Get indices of all values > 0.9 * final value: gt_9_10 = find(Y > 0.9*Yf); Subtract domain values to get time duration: tr = domain(Y(gt_9_10(1,1)))-... domain(Y(gt_1_10(1,1))) tr (a scalar) = 0.8 Get peak value of response: maxY = max(Y, {channels}); Index and time corresponding to peak value: maxtp = find(Y == maxY); tp = domain(Y(maxtp(1))) tp (a scalar) = 1.6 Convert peak value to a percentage > 1: Mp = round(10000*(maxY-1))/100 Mp (a scalar) = 23.
5 Classical Feedback Analysis The open-loop systems analyzed in Chapter 4, System Analysis, generally perform in a satisfactory manner only if the system model is very accurate and there are no external disturbances. These conditions usually are not met. Feedback presents an effective way to control the output of a system. The functions in this chapter address the problem of suitably controlling an open-loop plant through output feedback.
Chapter 5 Classical Feedback Analysis Because open-loop systems are generally easier to study and model than closed-loop systems, you want to design closed-loop systems based on information obtainable from the open-loop system. Root Locus In Chapter 4, System Analysis, you learned how the location of the system poles and zeros affects the stability of the system, so an effective feedback control design should take into account the closed-loop pole and zero locations.
Chapter 5 Classical Feedback Analysis rlocus( ) rlroots=rlocus(sys,K,{xmin,xmax,ymin,ymax,pattern,graph}) rlocus(sys,{xmin,xmax,ymin,ymax,pattern})# (interactive) The rlocus( ) function computes and draws root locus diagrams for continuous-time and discrete-time SISO systems. The first syntax, in which a vector of gain values is specified, generates a plot showing the closed-loop pole locations for each gain.
Chapter 5 Classical Feedback Analysis Figure 5-2. Root Locus of H for Gain K = 0.07 This syntax allows you to vary the root locus gain through an interactive form. Within this form, you can change the gain value through either a slider or an editable label where value corresponds to the current slider position. The slider range is automatically updated when the slider is moved to its maximum or minimum value, or when a gain value outside the current slider limits is entered into the editable label.
Chapter 5 Classical Feedback Analysis As the gain varies, small ✱’s appear on the locus indicating the closed-loop pole location for that choice of gain. The locus shown in Figure 5-2 shows that for small gain values the closed-loop system is stable, with all of its roots in the left half of the complex plane. Frequency Response and Dynamic Response The frequency response of a dynamic system is the output, or response, of a system given unit-amplitude, zero-phase sinusoidal input.
Chapter 5 Classical Feedback Analysis For discrete-time state-space systems with a sampling interval of T, the frequency response for each frequency point ω is shown in the following equation: H ( jw ) = C ( e jwT –1 I – A) B + D Algorithm The algorithm, based on [Lau86], uses a Hessenberg decomposition to simplify the previous equations and is quite robust. It finds matrices P and H such that A = PHP', where PP' = P'P = I and H is a Hessenberg matrix, and substitutes for A.
Chapter 5 Classical Feedback Analysis For an example of frequency response of a simple system, refer to Example 5-2. Given the single-input, single-output open-loop plant in Figure 5-3, where U(s) and Y(s) are the frequency domain input and output, respectively, you can examine its response characteristics and see how you can improve them using the frequency-response based control design functions in this chapter. U(s) s2(s (s + 0.5) + 2)(s + 10) Y(s) Figure 5-3.
Chapter 5 Classical Feedback Analysis because they can be used to assess the relative stability of a closed-loop system given the frequency response of the open-loop system. It should be noted that the open-loop system should be stable and minimum phase, having no right-half plane poles or zeros, for this type of analysis [Oga70].
Chapter 5 Classical Feedback Analysis Referring to the entire closed-loop system in Figure 5-1 as Gcl, the poles of Gcl are the roots of its denominator—that is, the values of s such that either of the following is true: 1 + G(s) = 0 G ( s ) = –1 The magnitude (absolute value) of G(s) is 1 at each pole of Gcl(s), and the phase (given by the four-quadrant arctangent) of G(s) is –180° at each pole of Gcl(s).
Chapter 5 Classical Feedback Analysis Referring to Figure 5-4, notice the additional lines drawn on the plots at the frequencies where the gain crosses the 0 dB line and where the phase crosses the 180° line. When the gain crosses the 0 dB line, the phase is about –168°, 12° away from –180°. So the phase margin is approximately 12°. Similarly, when the phase crosses the –180° line, the gain is about –44 dB (44 dB from the 0 dB line), and thus the gain margin is 44 dB.
Chapter 5 Classical Feedback Analysis gain margin phase margin Figure 5-4. Bode Plot Showing System Gain and Phase Margins These plots illustrate how the location of the system poles and zeros shapes the gain and phase curves. Each pole contributes a factor of –20 dB per decade (frequency interval from ω to 2ω). The two poles at zero cause the magnitude response of the system to start with a slope of –40 dB/decade. The zero at 0.5 radians/sec (about 0.08 Hz) contributes a factor of approximately 20 dB.
Chapter 5 Classical Feedback Analysis Each of these contributes a phase angle φ defined by: φ = atan ( ω ⁄ p n ) with ω and pn expressed in the same units, either radians per second or Hz, and using a four-quadrant arctangent function similar to that provided by atan2( ) in Xmath.
Chapter 5 Classical Feedback Analysis margin( ) loops over all the frequency points in the response and performs the following computation for phase and gain margins at each, denoting gain margin as Mg and phase margin as Mp: Δphase Δω Mp ( i ) = phase ( i ) + ------------------- ⎛ Δω – gain ( i + 1 ) ---------------⎞ Δω ⎝ Δgain⎠ ( gain ( i ) + Δgain ) Δω Mg ( i ) = – --------------------------------------------- ⎛ Δω – phase ( i + 1 ) -------------------⎞ ⎝ Δphase Δphase⎠ This loop finds all the frequen
Chapter 5 Classical Feedback Analysis Note margin( ) also returns the frequencies at which the phase crosses the –180° line and the gain crosses the 0 dB line. These results match the gain and phase margins shown graphically in Figure 5-4. nichols( ) [H,dB,Phase] = nichols(Sys,{F,keywords}) The nichols( ) function is another useful frequency domain tool for examining system performance in dynamic systems.
Chapter 5 Classical Feedback Analysis The result is shown in Figure 5-5. Figure 5-5. nichols( ) Gain-Phase Plot Nyquist Stability Analysis Nyquist analysis is a frequency domain method for examining system performance of dynamic systems. Nyquist plots typically consist of the real part of the frequency response plotted against the imaginary part of the response.
Chapter 5 Classical Feedback Analysis plant is open-loop stable, then there should be no encirclements. If the plant has one open-loop unstable pole, there should be one negative (counter-clockwise) encirclement. The stability criterion is most easily derived from the SISO transfer-function representation of a system. The Nyquist plot for a MIMO system consists of a set of plots, one for each output, each containing as many input frequency response curves as there are system inputs.
Chapter 5 Classical Feedback Analysis The Nyquist plot Xmath generates is complete only for the frequencies you specify. Ideally you would obtain a plot based on the frequency response from ω = 0 to ω = ∞. However, a good choice of frequency range usually comes close enough. When you have obtained the Nyquist plot from approximately w = 0 to ∞, you can reflect it about the real axis to get a complete plot of the open-loop frequency response from –∞ to +∞.
Chapter 5 Classical Feedback Analysis Figure 5-7. Nyquist Plot of the Open-Loop System for Frequencies from 0.01 Hz to 10 Hz Figure 5-8. Nyquist Plot of the Open-Loop System for Frequencies from 0.5 Hz to 5 Hz Xmath Control Design Module 5-18 ni.
Chapter 5 Classical Feedback Analysis By combining the information from the two plots, reflecting them across the real axis to account for the negative frequency response and augmenting them with a line closing the contour in the clockwise direction, you obtain the sketch of the encirclement pattern shown in Figure 5-9.
Chapter 5 Classical Feedback Analysis -0.52263 0.00336213 + 3.75217 j 0.00336213 - 3.75217 j -11.4841 Two of the poles of the closed-loop system are now unstable. Linear Systems and Power Spectral Density A key characteristic of the linear, time-invariant systems represented in Xmath is that the transfer function between a system input and a system output is just the Fourier transform of the response at that output to a delta impulse at that input.
Chapter 5 Classical Feedback Analysis values at which the power spectral density is to be computed and where dependent matrices are the input power spectral density matrix at each frequency. psd( ) computes a cross power spectral density matrix for each of a user–specified set of frequency values ω, returning them together in the PDM Yspec: Y spec = H × ( jw ) × U spec jw × H ( – jw ) psd( ) calls freq( ) internally to compute the frequency response, H, of the system.
6 State-Space Design The functions in this chapter are generally termed “modern control” tools. They are based on the state-space linear system representation, and employ methods which are generally applicable to both SISO and MIMO problems. For a review of the state-space system representation, refer to the State-Space System Models section of Chapter 2, Linear System Representation. The process of state-space control system design comprises several distinct steps.
Chapter 6 State-Space Design matrix B, then the mode of the system associated with the corresponding eigenvalue cannot be controlled with any input. You can think of this in the SISO transfer function case as a cancellation between a numerator and denominator root—where you cannot control the system in the direction of that root (mode).
Chapter 6 v + – u x = Ax + Bu y = Cx + Du State-Space Design y K Figure 6-1.
Chapter 6 State-Space Design X0 0 Input Names ----------Input 1 Output Names -----------Output 1 System is continuous T (a square matrix) = 2.22045e-16 0 -1 0 1 0 -1 0 2.22045e-16 nuc (a scalar) = 2 These results indicate that only the first state of the system corresponds to a controllable mode, and the remaining two are uncontrollable.
Chapter 6 State-Space Design Beginning with the basic state-space equations (the Du output term can be omitted without loss of generality): x· = Ax + Bu y = Cx you can obtain expressions for the successive derivatives of the output term, thus forming a complete description of the initial condition on the output: C y 0 0 0 u ·y = CA x + CB 0 0 u· 2 y·· CAB CB 0 u·· CA Generally, you term the following matrix C CA CA … 2 ... O = CA n–1 the observability matrix.
Chapter 6 State-Space Design x u x = Ax + Bu y Cx x x = A + Bu + Ly L y Cx + – y Figure 6-2. General Observer Block Diagram If the observability matrix is nonsingular, you will be able to put the eigenvalues (pole locations) of (A – LC), shown in Equation 6-4, anywhere you want. Thus, you can choose them to make x̃ decay to zero as quickly as possible.
Chapter 6 Example 6-2 State-Space Design Observability of a System A system is described by: A = [1,0,0.01;0,1,0;0,0,1]; B = [1,0,0]';C = [0.6,0.8,0];D = 0; Sys = system(A,B,C,D); Performing, [SysO,T,nuo] = observable(Sys); The system has 1 unobservable state This example indicates that one state of the system’s states corresponds to an unobservable mode, but that the other two are observable.
Chapter 6 State-Space Design minimal( ) [SysM,T,nuco] = minimal(Sys,{tol}) Because nonminimal systems are uncontrollable, unobservable, or both, you want to be able to compute the minimal realization for a given system. This comprises the controllable and observable parts of the dynamic system. The minimal( ) function calls both the controllable( ) and observable( ) functions, then extracts the part of the original system that is both controllable and observable.
Chapter 6 State-Space Design ans (a scalar) = 1 zeros(SysM) ans is null stair( ) [SysT,T,nc] = stair(Sys,tol) The stair( ) converts a dynamic system to staircase form. In the staircase form, the A and B system matrices are linearly transformed so that they are partitioned into controllable and uncontrollable parts. By duality, converting the transpose of a system into staircase form results in its being separated into observable and unobservable parts.
Chapter 6 State-Space Design Duality and Pole Placement The new state-update equation in the Controllability section and the Observability and Estimation section, the time-update equation in the Observability and Estimation section, along with the corresponding block diagrams in Figures 6-1 and 6-2, indicate how you can move the eigenvalues, or poles, of a minimal system through the choice of a feedback gain K (or L).
Chapter 6 State-Space Design poleplace( ) is unusual among Xmath’s modern control design functions in that only the A and B matrix variables are used as input, rather than a complete system variable. This is done because the other state-space matrices are not needed in the computation, and in many cases it is desirable to change or perturb the elements of the A and B matrices slightly to simulate actual conditions, without having to reformat the entire system.
Chapter 6 State-Space Design pairs as poles( ). For each pole value in poles( ), poleplace( ) forms a vector by subtracting the pole’s value from each diagonal element of S except for the last element (0). The resulting matrix is then divided by the corresponding value in the random complex vector. The complex value is padded with zeros to form a vector that is row compatible with the matrix.
Chapter 6 State-Space Design For continuous-time systems, x· = Ax + Bu the quadratic performance index takes the form: ∞ J = ∫ x' ( t ) u' ( t ) 0 R xx R xu x ( t ) dt R xu' R uu u ( t ) For the discrete case where the system is defined as a multistage process: x k + 1 = Ax k + Bu k the performance index is defined similarly except that a summation sign replaces the integral of the preceding quadratic performance index equation.
Chapter 6 State-Space Design The optimal estimator and regulator problems illustrate the principle of duality—that for any given system realization {A,B,C} there is a dual realization {A',C',B'} with related controllability and observability. Refer to the Duality and Pole Placement section.
Chapter 6 State-Space Design f u x Figure 6-4. Diagram of Plant for the Inverted Pendulum Problem Figure 6-4 shows the pendulum at φ = 0 and φ > 0. The distance of the cart from some initial reference point along the line of its motion is represented by the state variable x. You can measure the angle φ and the distance x easily—in fact, you will use measurements as your system outputs—but it is more difficult to obtain accurate measurements of the rate at which x and φ change.
Chapter 6 State-Space Design Ruu is a scalar because you have only one input for this particular model. [Kr,ev,P] = regulator(ipsys,Rxx,Ruu); Kr Kr (a row vector) = -348.778 Note -32.1056 -100 -27.3036 You will use this regulator gain later in designing a compensator. Linear Quadratic Estimator The LQR approach discussed in the preceding section is based on the assumption that the values of all the states are available.
Chapter 6 State-Space Design + x u x = Ax + Bu + Gw Cx + Du x x = Ax + Bu + Ke(y – y) + y Cx + Du Ke y + – e=y–y Figure 6-5. Diagram of the Estimator Representation estimator( ) inputs include the dynamic system Sys, and the noise intensity matrices Qxx , Qyy, or Qxy.
Chapter 6 State-Space Design numerical difficulties are encountered, the algorithm will attempt to determine whether or not the problem is well posed. Checks are made to determine the reachability and the positive definiteness or semipositive-definiteness of the covariance matrices.
Chapter 6 State-Space Design The discrete-time estimator follows from a similar system description, using the discrete-time difference equation representation of the system, as shown in the following equations. x k + 1 = Ax k + Bu k + Gω y k = Cx k + Du k + ν You obtain the discrete-time estimator by considering the state estimate at two separate stages. Begin with the assumption that an estimate of the state exists prior to each measurement of the output information.
Chapter 6 State-Space Design this measurement update, derived in [Kal60], are shown in the following equations. x̂ k = x k + K e ( y k – Cx k ) –1 –1 P k = ( M k + C'Q yy C ) –1 Substituting the system and noise matrices for the steady-state case, you solve the discrete Riccati equation to obtain P and thence Ke, as shown in Equation 6-7 and Equation 6-8.
Chapter 6 State-Space Design If you want the closed-loop system eigenvalues, compute them as the eigenvalues of A – KeC. For an example of how to design a state estimator for the inverted pendulum problem, refer to Example 6-6. Example 6-6 Designing a State Estimator for the Inverted Pendulum Problem Most systems have some level of internal process noise that affects the value of the states.
Chapter 6 State-Space Design a compensator is shown in Figure 6-6. This figure combines full-state regulator with gain Kr and state estimator with gain Ke.
Chapter 6 State-Space Design lqgcomp( ) SysC = lqgcomp(Sys,Kr,Ke,{direct}) The lqgcomp( ) function creates a dynamic compensator given a dynamic system having at least one state and the regulator and estimator gain matrices. The returned compensator SysC is always in state-space form. The regulator and estimator gains Kr and Ke need to have been calculated prior to the call to lqgcomp( ). You can use regulator( ) and estimator( ) to compute these gains if they are not already available.
Chapter 6 State-Space Design simulate the system’s response to a slow sine input, starting with the cart at rest and the pendulum initially held in the upright (φ = 0) position to obtain Figure 6-7: t = 0:0.
Chapter 6 State-Space Design Figure 6-7. f and x as a Function of Time, Starting from Zero, as a Result of a Sinusoidal Force Applied to the System Input Riccati Equation Riccati equations, which take one of two distinct forms, arise in a number of linear systems and controls problems. The best-known use is in the solution of the optimal regulator and estimator problems, as described in the Linear Quadratic Regulator section and the Linear Quadratic Estimator section.
Chapter 6 State-Space Design continuous-time Riccati equation, which is used if B and S are not specified, is: A′P + PA – PRP + Q = 0 Note The meaning of R is quite different in this case. riccati( ) [P,resid, Kr, ev] = riccati(A,Q,R,{B,tol,S,d}) Here, A can be either a matrix or a system object.
Chapter 6 0 0 0 State-Space Design -4 B 0 -2 0 1 C 0.211325 0.756044 0.000221135 0.330327 D 0 X0 0 0 0 0 System is continuous [P,resid] = riccati(A,Q,R); norm(A'*P+P*A-P*R*P+Q,1) ans (a scalar) = 2.53297e-13 R = 1e-5; [P,resid]=riccati(Sys,Q,R); norm(A'*P+P*A-P*B*inv(R)*B'*P +Q,1) ans (a scalar) = 2.52492e-13 The small residue indicates that the problem was well posed and the solution is reliable.
Chapter 6 State-Space Design 0 1 0 B 1 0 0 C 1 0 0 D 0 X0 0 0 0 System is discrete, sampling at 1 seconds. [P,resid]=riccati(Sys,Q,RD,B); norm(A'*P*A-P-A'*P*B*inv(RD+B'*P*B)*B'*P*A+Q,1) resid (a scalar) = 7.90593e-12 [P,resid]=riccati(Sys,Q,RD,B); norm(A'*P*A-P-A'*P*B*inv(RD+B'*P*B )* B'*P*A+Q,1) ans (a scalar) = 7.90593e-12 Steady-State System Response Using Lyapunov Equations The Lyapunov family of matrix equations are used in a number of control design problems.
Chapter 6 State-Space Design all the eigenvalues of the system A matrix are negative. The discrete Lyapunov equation is: AXA' + C = X (6-12) Analogously, the preceding equation has a unique solution X when λi(A)λj(A) ≠ 1 for all i and j. Again, this means that a unique X exists for a stable discrete-time system matrix A, because all eigenvalues of A have absolute value less than 1 in this case.
Chapter 6 State-Space Design or the following for the discrete case: y k = Cx k + Du k (6-16) Y = CXC' + DQD' These results derive from the Lyapunov method of stability analysis for linear systems. Steady state means that at some point the states no longer change. The derivative term x· approaches zero in the large for continuous systems, and xk+1 = xk for discrete-time ones. The state value vector x for which this is true is defined as the equilibrium state.
Chapter 6 State-Space Design of the diagonal elements of the Schur-decomposed A and B matrices sum to zero, a warning is given that the continuous equation solution may not be unique. A similar warning appears for the discrete equation solution if the product of any of the eigenvalues is 1. To solve the special Lyapunov equation, use the following syntax: lyapunov(A,C) Example 6-10 Lyapunov Equation Solutions The following examples each give results close to zero.
Chapter 6 State-Space Design Special Lyapunov Equation A × X + X × A' = – C A = [-4,10;2,7]; C = [.3,6;2,9]; X = lyapunov(A,C) X (a square matrix) = 1.1816 1.12431 -0.209028 -0.773611 A*X + X*A' + C; norm(A*X + X*A' + C,1) ans (a scalar) = 5.4956e-15 rms( ) [Yrms,Ycov] = rms(Sys,Ucov) The rms( ) function computes the root-mean-square response (average power at the system output) and the output covariance of a dynamic system driven by zero-mean white noise input.
Chapter 6 Example 6-11 State-Space Design rms( ) Response Sys = system([-2.3,0.01,5.1;0,-0.35,-2; 0,2,-.35],[1,.25,.25]',[1.34,0,0],0); w = logspace(0.01,1,50); Uspec = pdm(ones(w),w); [Ypsd,Yspec] = psd(Sys,Uspec); Balancing a Linear System Given a particular system model, the concept of model reduction centers on finding a lower-order model with similar input-output response characteristics. Typically this is assessed by comparing the impulse responses of the two systems [Moo81].
Chapter 6 State-Space Design For discrete-time systems, the integrals in the Wc and Wo equations are replaced by summation signs and the grammians are obtained as the solutions of the discrete-time Lyapunov equations: AW c A' + BB' = W c A'W o A + C'C = W o (6-19) The controllability grammian must be full-rank for the system to be completely controllable; similarly, the observability grammian must be full-rank for the system to be completely observable (refer to [Kai80]).
Chapter 6 State-Space Design and σ12 through σn2 are the singular values of the matrix H satisfying Σ2 = H'H. They are termed the Hankel singular values. The σk2 terms are ordered so that σ12 ≥ σ22 ≥ … ≥ σn2 ≥ 0. The balanced system essentially gives the best compromise between how well conditioned the system is with regard to controllability and observability.
Chapter 6 State-Space Design T is the transformation relating the states of the original system to the states of the balanced system x = Tx̂. Transforming to balanced coordinates can be useful in model reduction because the relative importance of the state to the system’s input/output performance is highlighted. balance( ) first finds the controllability and observability grammians using lyapunov( ) for both the discrete and continuous cases.
Chapter 6 State-Space Design and compare the condition numbers of the balanced system’s grammians: WcB=lyapunov(Ab,Bb*Bb'); WoB=lyapunov(Ab',Cb'*Cb); condition(WcB) ans (a scalar) = 12.7394 condition(WoB) ans (a scalar) = 12.7394 The condition numbers are now much smaller, and they are equal, indicating that the system is now equally well conditioned in terms of its controllability and observability.
Chapter 6 State-Space Design accurate. Given a variable Sys built from the matrices {A,B,C,D}, the modal decomposition SysMod is built from T–1 AT, T–1 B, CT, and D, where T is the transformation matrix to modal form. If you have complex poles, then T–1 AT is in block diagonal form. Initial conditions X0 also are transformed to T–1 X0.
Chapter 6 State-Space Design the derivative of x2 is set to zero, resulting in reduced-order state equations of the form: –1 –1 x· 1 = ( A 11 – A 12 A 22 A 21 )x 1 + ( B 1 – A 12 A 22 B 2 )u –1 –1 y = ( C 1 – C 2 A 22 A 21 )x 1 + ( D – C 2 A 22 B 2 )u In the discrete case, x2k + 1 is taken to be equal to x2k so that the state equations become: –1 x 1k + 1 = [ A 11 – A 12 ( A 22 – I ) A 21 ]x 1k + –1 [ B 1 – A 12 ( A 22 – I ) B 2 ]u k –1 y = [ C 1 – C 2 ( A 22 – I ) A 21 ]x 1k + –1 [ D – C 2 ( A 22 –
Chapter 6 State-Space Design 0 0 0 0.924711 B -0.00116788 0.00272531 0.00334243 -0.00162497 C 1 0.866186 -0.848754 -1.0118 D 0 X0 0 0 0 0 State Names ----------State 1 State 2 State 3 State 4 Input Names ----------Input 1 Output Names -----------Output 1 System is discrete, sampling at 0.2 seconds. T (a square matrix) = 1 0 0 0 0.866186 -0.668844 -0.641745 0.499722 -0.71998 -0.653294 0 -0.17991 -0.370059 0 0.043691 -0.15629 eig(A) ans (a column vector) = 0.37 0.52 0.665289 0.
Chapter 6 State-Space Design 0 0.52 1.94289e-16 6.66134e-16 0 0 0.665289 4.44089e-16 0 0 0 0.924711 SysMR = mreduce(SysM, [1,2,4]) SysMR (a state space system) = A 0.37 0 0 0 0.52 0 0 0 0.924711 B -0.00116788 0.00272531 -0.00162497 C 1 0.866186 -1.0118 D -0.00847566 X0 0 0 0 0 State Names ----------State 1 State 2 State 4 Input Names ----------Input 1 Output Names -----------Output 1 plot(step(SysM, 0:.2:10)) plot(step(SysMR, 0:.
Chapter 6 State-Space Design Figure 6-8. Modal System and Reduced Modal System Xmath Control Design Module 6-42 ni.
Technical References A [BeV88] T. Beelen, P. Van Dooren, “An improved algorithm for the computation of Kronecker’s canonical form of a singular pencil,” Linear Algebra and Applications, 105, pages 9–65, 1988. [BH75] A.E. Bryson and Y.C. Ho, Applied Optimal Control, Hemisphere Publishing Corporation, Washington, D.C., 1975. [ChB84] R.V. Churchill and J.W. Brown, Complex Variables and Applications, McGraw-Hill Book Company, New York, 1984. [DeS74] C.A. Desoer and J.D.
Appendix A Technical References [Kai80] T. Kailath, Linear Systems, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1980. [Kai81] T. Kailath, Lectures on Wiener and Kalman Filtering, Springer-Verlag, New York, 1981. [Lau79] A.J. Laub, “A Schur method for solving algebraic Riccati equations,” IEEE Transactions on Automatic Control, AC-24, pages 913–621, 1979. [Lau86] A.J.
Technical Support and Professional Services B Visit the following sections of the National Instruments Web site at ni.com for technical support and professional services: • Support—Online technical support resources at ni.
office Web sites, which provide up-to-date contact information, support phone numbers, email addresses, and current events.
Index A connect, 3-8 connection, 3-1 parallel, 3-1, 3-2 series, 3-2 using * operator, 1-16, 3-2 constant gain feedback, 3-8 constant magnitude and phase loci, 5-14 continuous equivalent to a discrete system, 2-18 continuous system, 4-13 analysis, 2-7 checking for, 2-12 continuous time Riccati equation, 6-13, 6-25 controllability, 6-4, 6-35 grammians, 6-30 matrix, 6-2 controllable, 6-3 controllable partition of a state-space system, 6-3 conventions used in the manual, iv coordinate transformation, 6-7 corne
Index estimator system, 1-17 examples (NI resources), B-1 expectation operator, 6-17 exponential discretization method, 1-19, 2-16 converting to continuous equivalent, 2-19 from a continuous system (example), 1-18 discrete-time Riccati equation, 6-25 discretize, 2-13 discretizing a system backward rectangular method, 2-14 forward rectangular method, 2-14 hold equivalence methods, 2-15 pole-zero matching, 2-15 trapezoid method, 2-14 Tustin’s method, 2-14 using ztransform for zero-order hold, 2-15 with expo
Index I M impulse, 4-13 input, 4-13 continuous time, 4-13 discrete time, 4-13 response, 4-13, 6-33 initial, 4-16 initial conditions, 4-14, 4-16 input disturbance matrix, 6-17 names, extracting, 2-11 instability, 5-9 instrument drivers (NI resources), B-1 internally balanced system (definition), 6-34 inverse, 3-3 inverted pendulum problem, 6-21 magnitude, 5-5, 5-8 makecontinuous, 2-17 verifying discretization with, 2-18 makepoly, 2-3 margin, 5-12 Markov parameters, 4-13 matrix controllability, 6-2 inputs
Index period, 2-10 phase, 5-5, 5-11 f, 5-8 margin, 5-8, 5-9, 5-12 rate of change, 5-12 tracking, 5-6 poleplace, 6-10 finding feedback gains, 1-14 finding observer gains, 1-15 poles, 1-1, 2-2, 4-3 open-loop, 5-15 placement, 6-1 problem, 6-10 unstable open-loop, 5-15 pole-zero cancellation, 1-11 matching, 2-15 power spectral density, 5-20, 5-21 process noise, 6-16 programming examples (NI resources), B-1 psd, 5-20 steady-state response to white noise, 6-29 white process, 6-23 nomenclature, 1-2 numden, 2-10
Index state-space system, 2-5 checking for, 2-12 controllable partition of, 6-3 convert to transfer function form, 1-5 decompose with abcd, 2-8 discrete, creating, 2-5 extracting transfer function polynomials with numden, 2-10 modes of, 6-37 observable partition of, 6-6 steady state, 6-30 system response, 6-28 step, 1-9, 4-18 step response, 4-18 support, technical, B-1 system, 2-6 analysis, 2-7 cascaded, 5-21 connections, 3-1 continuous, 4-13 controllability, 6-4 controllable, 6-1 discrete, 4-13 general in
Index T trapezoid method for discretization, 2-14 troubleshooting (NI resources), B-1 Tustin’s discretization method, 2-14 technical support, B-1 time domain simulation, general, 4-10 update, 6-19 training and certification (NI resources), B-1 transfer function checking for, 2-12 coefficients form, 2-3 converted to state space before decomposition, 2-9 creating, 2-3 form, definition, 4-2 formed from partial fractions, 4-9 pole-zero-gain form, 2-4 polynomials, extracting from state-space system, 2-10 syst
