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. National Instruments does not warrant that the operation of the software shall be
uninterrupted or error free.
A Return Material Authorization (RMA) number must be obtained from the factory and clearly marked on the outside of the package before any
equipment will be accepted for warranty work. National Instruments will pay the shipping costs of returning to the owner parts which are covered by
warranty.
National Instruments believes that the information in this document is accurate. The document has been carefully reviewed for technical accuracy. In
the event that technical or typographical errors exist, National Instruments reserves the right to make changes to subsequent editions of this document
without prior notice to holders of this edition. The reader should consult National Instruments if errors are suspected. In no event shall National
Instruments be liable for any damages arising out of or related to this document or the information contained in it.
EXCEPT AS SPECIFIED HEREIN, NATIONAL INSTRUMENTS MAKES NO WARRANTIES, EXPRESS OR IMPLIED, AND SPECIFICALLY DISCLAIMS ANY WARRANTY OF
MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. CUSTOMER’S RIGHT TO RECOVER DAMAGES CAUSED BY FAULT OR NEGLIGENCE ON THE PART OF NATIONAL
INSTRUMENTS SHALL BE LIMITED TO THE AMOUNT THERETOFORE PAID BY THE CUSTOMER. NATIONAL INSTRUMENTS WILL NOT BE LIABLE FOR DAMAGES RESULTING
FROM LOSS OF DATA, PROFITS, USE OF PRODUCTS, OR INCIDENTAL OR CONSEQUENTIAL DAMAGES, EVEN IF ADVISED OF THE POSSIBILITY THEREOF. This limitation of
the liability of National Instruments will apply regardless of the form of action, whether in contract or tort, including negligence. Any action against
National Instruments must be brought within one year after the cause of action accrues. National Instruments shall not be liable for any delay in
performance due to causes beyond its reasonable control. The warranty provided herein does not cover damages, defects, malfunctions, or service
failures caused by owner’s failure to follow the National Instruments installation, operation, or maintenance instructions; owner’s modification of the
product; owner’s abuse, misuse, or negligent acts; and power failure or surges, fire, flood, accident, actions of third parties, or other events outside
reasonable control.
Copyright
Under the copyright laws, this publication may not be reproduced or transmitted in any form, electronic or mechanical, including photocopying,
recording, storing in an information retrieval system, or translating, in whole or in part, without the prior written consent of National
Instruments Corporation.
National Instruments respects the intellectual property of others, and we ask our users to do the same. NI software is protected by copyright and other
intellectual property laws. Where NI software may be used to reproduce software or other materials belonging to others, you may use NI software only
to reproduce materials that you may reproduce in accordance with the terms of any applicable license or other legal restriction.
Trademarks
MATRIXx™, National Instruments™, NI™, ni.com™, and Xmath™ are trademarks of National Instruments Corporation.
Other product and company names mentioned herein are trademarks or trade names of their respective companies.
Members of the National Instruments Alliance Partner Program are business entities independent from National Instruments and have no agency,
partnership, or joint-venture relationship with National Instruments.
Patents
For patents covering National Instruments products, refer to the appropriate location: Help»Patents in your software, the patents.txtfile
on your CD, or ni.com/patents.
WARNING REGARDING USE OF NATIONAL INSTRUMENTS PRODUCTS
(1) NATIONAL INSTRUMENTS PRODUCTS ARE NOT DESIGNED WITH COMPONENTS AND TESTING FOR A LEVEL OF
RELIABILITY SUITABLE FOR USE IN OR IN CONNECTION WITH SURGICAL IMPLANTS OR AS CRITICAL COMPONENTS IN
ANY LIFE SUPPORT SYSTEMS WHOSE FAILURE TO PERFORM CAN REASONABLY BE EXPECTED TO CAUSE SIGNIFICANT
INJURY TO A HUMAN.
(2) IN ANY APPLICATION, INCLUDING THE ABOVE, RELIABILITY OF OPERATION OF THE SOFTWARE PRODUCTS CAN BE
IMPAIRED BY ADVERSE FACTORS, INCLUDING BUT NOT LIMITED TO FLUCTUATIONS IN ELECTRICAL POWER SUPPLY,
COMPUTER HARDWARE MALFUNCTIONS, COMPUTER OPERATING SYSTEM SOFTWARE FITNESS, FITNESS OF COMPILERS
AND DEVELOPMENT SOFTWARE USED TO DEVELOP AN APPLICATION, INSTALLATION ERRORS, SOFTWARE AND HARDWARE
COMPATIBILITY PROBLEMS, MALFUNCTIONS OR FAILURES OF ELECTRONIC MONITORING OR CONTROL DEVICES,
TRANSIENT FAILURES OF ELECTRONIC SYSTEMS (HARDWARE AND/OR SOFTWARE), UNANTICIPATED USES OR MISUSES, OR
ERRORS ON THE PART OF THE USER OR APPLICATIONS DESIGNER (ADVERSE FACTORS SUCH AS THESE ARE HEREAFTER
COLLECTIVELY TERMED “SYSTEM FAILURES”). ANY APPLICATION WHERE A SYSTEM FAILURE WOULD CREATE A RISK OF
HARM TO PROPERTY OR PERSONS (INCLUDING THE RISK OF BODILY INJURY AND DEATH) SHOULD NOT BE RELIANT SOLELY
UPON ONE FORM OF ELECTRONIC SYSTEM DUE TO THE RISK OF SYSTEM FAILURE. TO AVOID DAMAGE, INJURY, OR DEATH,
THE USER OR APPLICATION DESIGNER MUST TAKE REASONABLY PRUDENT STEPS TO PROTECT AGAINST SYSTEM FAILURES,
INCLUDING BUT NOT LIMITED TO BACK-UP OR SHUT DOWN MECHANISMS. BECAUSE EACH END-USER SYSTEM IS
CUSTOMIZED AND DIFFERS FROM NATIONAL INSTRUMENTS' TESTING PLATFORMS AND BECAUSE A USER OR APPLICATION
DESIGNER MAY USE NATIONAL INSTRUMENTS PRODUCTS IN COMBINATION WITH OTHER PRODUCTS IN A MANNER NOT
EVALUATED OR CONTEMPLATED BY NATIONAL INSTRUMENTS, THE USER OR APPLICATION DESIGNER IS ULTIMATELY
RESPONSIBLE FOR VERIFYING AND VALIDATING THE SUITABILITY OF NATIONAL INSTRUMENTS PRODUCTS WHENEVER
NATIONAL INSTRUMENTS PRODUCTS ARE INCORPORATED IN A SYSTEM OR APPLICATION, INCLUDING, WITHOUT
LIMITATION, THE APPROPRIATE DESIGN, PROCESS AND SAFETY LEVEL OF SUCH SYSTEM OR APPLICATION.
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The sequence File»Page Setup»Options directs you to
pull down the File menu, select the Page Setup item, and select Options
from the last dialog box.
This icon denotes a note, which alerts you to important information.
bold
Bold text denotes items that you must select or click in the software, such
as menu items and dialog box options. Bold text also denotes parameter
names.
italic
Italic text denotes variables, emphasis, a cross-reference, or an introduction
to a key concept. Italic text also denotes text that is a placeholder for a word
or value that you must supply.
monospace
Text in this font denotes text or characters that you should enter from the
keyboard, sections of code, programming examples, and syntax examples.
This font is also used for the proper names of disk drives, paths, directories,
programs, subprograms, subroutines, device names, functions, operations,
variables, filenames, and extensions.
monospace italic
Italic text in this font denotes text that is a placeholder for a word or value
that you must supply.
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Commonly Used Nomenclature ......................................................................1-2
Related Publications........................................................................................1-3
MATRIXx Help...............................................................................................1-3
Control Design Tutorial.................................................................................................1-4
Helicopter Hover Problem: State Feedback and Observer Design .................1-13
Chapter 2
Linear Systems Represented in Xmath..........................................................................2-1
Transfer Function System Models.................................................................................2-2
abcd( )..............................................................................................................2-8
discretize( )......................................................................................................2-13
Numerical Integration Methods: forward, backward, tustins ...........2-14
Pole-Zero Matching: polezero ..........................................................2-15
Z-Transform: ztransform...................................................................2-15
Hold Equivalence Methods: exponential and firstorder ...................2-15
makecontinuous( ) ...........................................................................................2-17
© National Instruments Corporation
v
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Contents
Chapter 3
Linear System Interconnection Operators..................................................................... 3-1
Linear System Interconnection Functions..................................................................... 3-4
Algorithm.......................................................................................... 3-5
append( ) ......................................................................................................... 3-6
Algorithm.......................................................................................... 3-10
Chapter 4
System Analysis
System Stability: Poles and Zeros................................................................................. 4-2
combinepf( ).................................................................................................... 4-9
Impulse Response of a System...................................................................................... 4-13
impulse( ) ........................................................................................................ 4-13
deftimerange( )................................................................................................ 4-15
initial( )............................................................................................................ 4-16
Chapter 5
Feedback Control of a Plant Model............................................................................... 5-1
Root Locus..................................................................................................................... 5-2
rlocus( ) ........................................................................................................... 5-3
Frequency Response and Dynamic Response ............................................................... 5-5
freq( )............................................................................................................... 5-5
Algorithm.......................................................................................... 5-6
Xmath Control Design Module
vi
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
margin( )..........................................................................................................5-12
nichols( )..........................................................................................................5-14
nyquist( )..........................................................................................................5-16
psd( )................................................................................................................5-20
Chapter 6
Controllability................................................................................................................6-1
observable( )....................................................................................................6-6
Minimal Realizations.....................................................................................................6-7
minimal( )........................................................................................................6-8
Duality and Pole Placement...........................................................................................6-10
poleplace( )......................................................................................................6-10
regulator( ).......................................................................................................6-14
Algorithm..........................................................................................6-30
rms( ) ...............................................................................................................6-32
Balancing a Linear System ............................................................................................6-33
balance( ) .........................................................................................................6-35
Modal Form of a System ...............................................................................................6-37
modal( ) ...........................................................................................................6-37
mreduce( )........................................................................................................6-38
© National Instruments Corporation
vii
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Contents
Appendix A
Appendix B
Technical Support and Professional Services
Index
Xmath Control Design Module
viii
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Using This Manual
This manual provides an overview of different aspects of linear systems
analysis, describes the Xmath Control Design function library, and gives
examples of how you can use Xmath to solve problems rapidly. It also
explains how you can represent and analyze linear systems in Xmath and
provides a brief syntax listing and supplementary algorithm information for
behavior are provided in the Xmath Help.
Document Organization
•
Chapter 1, Introduction, 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 2, Linear System Representation, describes the types of linear
systems that can be represented within Xmath. In addition, it discusses
the implementation of systems as objects–data structures
encompassing different information fields. The Xmath functions for
creating a system or extracting its components are part of the general
Xmath package and not exclusive to the Control Design Module, but
they are used so extensively that they warrant a detailed treatment here.
This chapter also discusses the functions you can use to check for
© National Instruments Corporation
1-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
particular system properties or to change the format of a system. These
topics include continuous/discrete system conversion, as well as
•
•
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.
system stability and time-domain analysis. These include poles, zeros,
and residue. The chapter moves from the discussion of time-domain
functions for obtaining impulse and step responses, as well as
examining system response to arbitrary initial conditions. In addition,
the General Time-Domain Simulation section discusses a
mathematically natural syntax for time-domain system simulation
with any input.
•
•
Chapter 5, Classical Feedback Analysis, discusses topics pertaining
to classical feedback-based control design. These include root locus
techniques and functions for frequency-domain analysis of
closed-loop systems, given open-loop system descriptions.
Chapter 6, State-Space Design, focuses on modern control. Beginning
with the topics of system controllability and observability, it covers
general pole placement, linear quadratic control, and system
balancing.
Bibliographic References
Throughout this document, bibliographic references are cited with
bracketed entries. For example, a reference to [DeS74] corresponds to
a document published by Desoer and Schulman in 1974. For a table of
bibliographic references, refer to Appendix A, Technical References.
Commonly Used Nomenclature
This manual uses the following general nomenclature:
•
•
Matrix variables are generally denoted with capital letters; vectors are
represented in lowercase.
G(s) is used to denote a transfer function of a system where s is the
Laplace variable. G(q) is used when both continuous and discrete
systems are allowed.
Xmath Control Design Module
1-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Related Publications
For a complete list of MATRIXx publications, refer to Chapter 2,
MATRIXx Publications, Online Help, and Customer Support, of the
MATRIXx Getting Started Guide. The following documents are particularly
useful for topics covered in this manual:
•
•
•
•
•
•
•
•
•
•
MATRIXx Getting Started Guide
Xmath User Guide
Control Design Module
Interactive Control Desing Module
Interactive System Identification Module, Part 1
Interactive System Identification Module, Part 2
Model Reduction Module
Optimization Module
Robust Control Module
Xμ Module
MATRIXx Help
Control Design Module function reference information is available in the
MATRIXx Help. The MATRIXx Help includes all Control Design functions.
Each topic explains a function’s inputs, outputs, and keywords in detail.
Refer to Chapter 2, MATRIXx Publications, Online Help, and Customer
Support, of the MATRIXx Getting Started Guide for complete instructions
on using the Help feature.
© National Instruments Corporation
1-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
The models in this tutorial are adapted from the studies in [ShH92], of the
equations presented in [FPE87], for the longitudinal motion of a helicopter
near hover, and in [HW91], for the inverted-wedge-balancing problem.
Helicopter Hover Problem: An Ad Hoc Approach
[FPE87] gives this state-space model for the longitudinal motion of the
helicopter:
·
q
–0.4 0 –0.01
q
θ
v
6.3
0
·
=
+
δ
1
0
0
θ
·
–1.4 9.8 –0.02
9.8
v
q
0 0 1 θ
v
y =
letting the state variables q, θ, and v represent the helicopter’s pitch rate,
pitch angle, and horizontal velocity, respectively. The input control to the
system is the rotor tilt angle, δ.
You can store the information that this model provides in an Xmath
state-space system object:
A = [-0.4,0,-0.01;1,0,0;-1.4,9.8,-0.02];
B = [6.3;0;9.8];
C = [0,0,1];
D = 0;
ssys = system(A,B,C,D,
outputNames="Horizontal v",
stateNames =["Pitch Rate", "Pitch Angle",
"Horizontal v"]})
Xmath Control Design Module
1-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
ssys (a state space system) =
A
-0.4 0 -0.01
1
0
0
-1.4 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.149274)
initial integrator outputs
0
0
0
© National Instruments Corporation
1-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.4975 j
Try to stabilize the system using feedback compensation. You have two
major performance goals to achieve through your controller design: first,
the system must be closed-loop stable, and second, you want the system
output to track a unit step input. To begin, put two compensators in the
feedforward path of the closed-loop system. Figure 1-1 is a closed-loop
block diagram of helicopter system H(s) with compensators K1(s) and K2(s)
in the feedforward path.
U(s)
Y(s)
+
–
G(s)
K1(s)
K2(s)
Figure 1-1. Block Diagram of Helicopter System H(s) with Compensators K1(s) and
K2(s) in the Feedforward Path
Xmath Control Design Module
1-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Note 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.
This compensator must have a pole for realizability, so you add one at –10,
which is far enough away that its effect on dynamic response will be small
compared to that of the system’s other modes. In addition, you need to add
a zero to the left of the positive (and unstable) poles to pull the closed-loop
system roots into the left half plane. Choose s = 0 for the zero location and,
again, select a corresponding pole at –10. Call this second compensator
K2(s). To create these two compensators:
K1s=polynomial(ol_poles(3),"s")/...
polynomial(-10,"s");
K2s=polynomial(0,"s")/polynomial(-10,"s");
You then can cascade them in series with the original system Gs(or ssys)
and examine the locus of closed-loop roots for varying total compensator
gain Kc. The poles out at –10 have a smaller effect on the system dynamics
than do the poles closer to the origin, so you can use the optional
rlocus( )keywords to zoom in on the part of the locus nearer the origin.
rlocus(K2s*K1s*Gs, {xmin = -2, xmax = 2})
The single input syntax activates interactive mode. A user interface lets
you change the gain through the Feedback Loop Gain slider and button.
The graphics window shows the closed-loop locus as a solid line, with the
open-loop poles shown as large xs and the open-loop zeros shown as Os.
Increase the gain by moving the slider; notice the asterisks (*) denoting
closed-loop pole location moving along the locus. The system is maximally
stable with total compensator gain Kc = 2 as shown in Figure 1-2. In this
figure, small xs denote the pole location for Kc = 2, and the root locus gain
window shows settings.
© National Instruments Corporation
1-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
Close the loop using the single-input syntax of feedback( ), which
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. You
want the output (the helicopter velocity) to track the step input provided as
the rotor tilt angle, not zero out its effects over time (which would be an
appropriate response if the input corresponded to a disturbance). This
results from the compensator zero at s = 0 in the forward path of the
feedback loop.
Figure 1-3. Helicopter Velocity Response to a Step Input at the Rotor
Instead, you now place K1(s) in the forward path and K2(s) in the feedback
path, so that the closed-loop system now has the configuration shown in
Figure 1-4.
© National Instruments Corporation
1-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
The gains Kc1 and Kc2 are defined such that Kc = 2 = Kc1Kc2 and the
closed-loop transfer function Tc1(s) is unity at s = 0(DC).
The closed-loop transfer function is represented by:
Kc1K1(s)G(s)
Tcl(s) = ------------------------------------------------------------------
1 + Kc1Kc2K1(s)K2(s)G(s)
You can find the values of the individual transfer functions at s = 0 using
freq( ), and then substitute to solve the previous equation:
a = makematrix(freq(K1s*Gs,0));
b = makematrix(freq(K1s*K2s*Gs,0));
Solving:
Kc1 = (1+2*b)/a
Kc1 (a scalar) = 0.0241778
Kc2 = 2/Kc1
Kc2 (a scalar) = 82.7206
You now call feedback( )again, this time using its second input
argument to indicate that the outputs of the first input system (forward path)
are fed back as the inputs to the second system (feedback path) in a
negative-feedback loop.
cl_syscomp2 = feedback(Kc1*K1s*Gs, Kc2*K2s);
Xmath Control Design Module
1-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
Because cl_syscomp2contains 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.4975 j
-10
Now, examine the step response as shown in Figure 1-5.
vcomp = step(cl_syscomp2m, 0:.1:25);
plot (vcomp, {xlab = "Time",
ylab = "Horizontal Velocity"})
© National Instruments Corporation
1-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Helicopter Hover Problem: State Feedback and Observer Design
The approach taken in the Helicopter Hover Problem: An Ad Hoc
be used in practice because uncertainty in modeling generally precludes
exact knowledge of the location of the pole one plans to cancel.
Another approach is to feed the information obtained from the states back
to the inputs through gains calculated to relocate the closed-loop poles.
Refer to the Controllability section of Chapter 6, State-Space Design,
for more information. For this approach, you first need to verify that your
system is controllable and observable. When you have confirmed that it
is—that there are no hidden modes—you can design a full-state feedback
control law that will place the system eigenvalues at values that will yield
a stable system. Because the system is observable, you then can design an
estimator to yield estimates for the missing states. Again, you will require
that your system track a step input.
© National Instruments Corporation
1-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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, Gsis controllable. This
means that you can use feedback through appropriately-sized gains to
position the system’s closed-loop poles anywhere you want. If you choose
the three poles to be moved to –1 j and –2, you get the following set of
gains:
clp = [-1+jay, -2];
Kfsb = poleplace(A,B,clp)
Note poleplace( )does not require you to list both poles in a conjugate pair.
If you assume that the outputs of the system are just the values of all the
states, you can draw the open-loop system block diagram as shown in
Figure 1-7. In this figure, the feedback path is shown in dotted lines and the
open-loop system in solid lines.
r
u
x
+
–
⋅
x = Ax + Bu
Kfsb
Figure 1-7. Full-State Feedback Regulator
Because you do not have access to all three states—only one, the horizontal
velocity, is returned as an output—you need to estimate the other states,
thus implementing an observer-based controller. The block diagram for the
observer and controller together is shown in Figure 1-8.
Estimator
Control
Kfsb
Plant
u
y
x
⋅
x = Ax + Bu
y = Cx + Du
x = (A – LC)x + Bu + Ly
Figure 1-8. Complete Controller and Estimator
Xmath Control Design Module
1-14
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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( ). Lneeds 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. To distinguish the estimated states from the “true” states,
you can use the +operator to append the string estto the estimated state
names, as shown in this example.
[,,osNames] = names(sys_obc);
estNames = osNames + [" (est)"," (est)"," (est)"];
estNames'?
ans (a column vector of strings) =
Pitch Rate (est)
Pitch Angle (est)
Horizontal v (est)
You can append these modified names to sys_obc:
sys_obc = system(sys_obc,{stateNames=estNames});
Then you close the loop and verify that the closed-loop poles are all in the
left-half plane.
sys_cl = feedback(ssys, sys_obc);
poles(sys_cl)
ans (a column vector) =
-1 + 1 j
-1 - 1 j
-2
-4
-3 + 3 j
-3 - 3 j
© National Instruments Corporation
1-15
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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_clby 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. Step Response for Observer-Based Design
The system response is quite good, implying that your state estimates were
satisfactory. You can do some further simulation, this time returning all the
states directly from the original plant, and get a graphical picture of how the
estimates track the actual states. First, you need to create the closed-loop
system with all states available.
The abcd( )function extracts the A, B, C, and D matrices from a system
object. When you call it here, all you are interested in is the closed-loop
A matrix, so you do not need to extract the other state-space matrices.
A_cl = abcd(sys_cl);
Xmath Control Design Module
1-16
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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 X0to set the real state values to
[1,2,3] and the estimated state values to [–1,–2,–3].
By simulating with a general input over two seconds, you can see how long
it takes for the state values provided by the estimator to correct the incorrect
initial conditions and track the real state values.
[,,allStates]=names(sys_cl);
sys_est = system(A_cl, zeros(6,1), eye(6,6), ...
zeros(6,1),{x0 = [1,2,3,-1,-2,-3], ...
stateNames = allStates});
state_resp = sys_est*pdm(ones(100,1), 0:(2/99):2);
Plot the results, referring to Figure 1-10:
plot(state_resp,{strip=2,xlab="Time",
legend=["State","State estimate"]})
Even in the relatively short time span of this simulation, the estimates and
the real states quickly converge.
© National Instruments Corporation
1-17
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
Figure 1-10. Multiple Plots Showing Time Needed for States to be Correctly Tracked
by Estimator, Given Incorrect Initial Values
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. Here you take the second
approach with the problem solved in the Helicopter Hover Problem: State
Feedback and Observer Design section.
A guideline for choosing a sample rate for a system to be discretized is that
it be significantly less than the smallest time constant of the continuous
system divided by π.
Look at the open-loop pole magnitudes of your original open-loop
continuous-time system ssys:
max(abs(ol_poles))
ans (a scalar) = 0.656513
Xmath Control Design Module
1-18
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
db = 20*log10(abs(f)); ph = (180/pi)*atan2(f);
plot([db;ph;dbd;phd],{strip=2,xlog,
ylab = ["Gain (dB)";"Phase (deg)"],
x_lab = "Frequency (Hz)",
legend = ["ssys";"ssysd"]})
In Figure 1-11 you can see the frequency responses match closely,
indicating that this discretization method captures the continuous system’s
dynamics accurately.
Figure 1-11. Frequency Response of ssys and Its Discrete Equivalent ssysd
© National Instruments Corporation
1-19
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.01:10);
plot (v_cld, {xlab = "Time", ylab = "Magnitude"})
The resulting response is shown in Figure 1-12.
Inverted Wedge-Balancing Problem: LQG Control
[HW91] discusses an approach to balancing an inverted wedge by
controlling the location of a sliding mass along the inside of the wedge.
This example illustrates use of optimal control with a multi-input,
multi-output (MIMO) system. This approach is based on minimizing a
quadratic performance index with weight values based on the natural
constraints of the system.
Xmath Control Design Module
1-20
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 1
Introduction
The linearized state-space equations, including the actuator and sensor
dynamics, are as follows:
·
θ
x
θ
0
0
0
0
1
0
0
1
0
0
0
·
x
=
+
u
·
··
15.54 –10.93 0
0
θ
θ
·
··
–5.31
0
0 –16.24
1.96
x
x
θ
x
57.29
0
0
0 0
y =
·
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.
A = [0,0,1,0;0,0,0,1;
15.54,-10.93,0,0;
-5.31,0,0,-16.24];
B = [0,0,0,1.96]';
C = [57.29,0,0,0;0,29.9,0,0];
D = [0;0];
states = ["Angle", "Mass Position",
"Angular Velocity","Mass Velocity"];
wsys = system(A,B,C,D,
{inputNames= "Voltage",
stateNames = states,
outputNames=["Measured Angle","Measured
Position"]});
You need to ensure that you have no uncontrollable or unobservable modes
of the system:
[,,nuco] = minimal(wsys)
nuco (a scalar) = 0
Because there are no uncontrollable or unobservable states, you can
proceed with the design of a regulator and estimator. The weighting matrix
used here in designing the regulator reflects the desire to bring the value of
the first state, the angle with the vertical, to zero as quickly as possible.
© National Instruments Corporation
1-21
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
[Kr,EVr,Pr] = regulator(wsys,diag([1e8,1,1,1]),1);
You then can verify that the regulator gain Krcan be used with full-state
feedback to control this system by using an identity matrix for C to feed
back the states:
[no, ni, ns] = size(wsys);
augwsys = system(A,B,eye(ns, ns),[]);
creating the compensator (which is a system object, though it has no states
and thus has NULLA, B, and C matrices) with the gains Kr:
comp = system([],[],[],Kr);
and feeding back the states:
wsysreg = feedback(augwsys, comp);
You then can observe the system response to a sustained disturbance by
simulating a five-second step response:
stepreg = step(wsysreg, 0:0.01:5);
plot (stepreg, {legend=states,
xlab="Time",ylab="Magnitude",
title="System Step Response with "+...
"Full State Availability",!grid})
The resulting plot is shown in Figure 1-13.
Xmath Control Design Module
1-22
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The output weight matrix
reflects your higher priority on the wedge angle than position.
The following steps generate the plot shown in Figure 1-14:
[Ke,EVe,Pe] = estimator(wsys,
diag([1e-3,1,1e-3,10]), diag([14,0.01]));
wcomp = lqgcomp(wsys,Kr,Ke);
wlqg = feedback(wsys,wcomp);
resp = step(wlqg,0:0.01:3);
plot (resp, {legend = names(wlqg),
xlab = "Time",ylab = "Magnitude",
title = "Observer-Controller System "+...
"Step Response",!grid})
© National Instruments Corporation
1-23
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. This chapter outlines the
types of linear systems the system object represents and then discusses the
implementation of a system within Xmath. The functions used to create a
system object and to extract data from this object are an intrinsic part of the
object and are also described. Finally, this chapter discusses the functions
check( ), discretize( ), and makecontinuous( ), which use
information stored in the system object to convert systems from one
representation to another.
Linear Systems Represented in Xmath
Xmath handles finite-dimensional, linear, and time-invariant linear
systems in both discrete and continuous time. These systems take one
of the forms shown in Table 2-1.
Table 2-1. Summary of Linear Systems
System Type
State-spec
Continuous Time
Discrete Time
xk + 1 = Axk + Buk
·
x = Ax + Bu
y = Cx + Du
yk = Cxk + Duk
Transfer function
H(s) = C(sI – A)–1B + D
H(z) = C(zI – A)–1B + D
The transfer function representation can be used to describe single-input,
single output (SISO) systems only; there are no restrictions on the number
of input and outputs that can be specified for a state-space system. All of
these systems can be created using the Xmath system( )function.
© National Instruments Corporation
2-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. Both forms of transfer functions are written
with positive coefficients, each higher order terms having successively
larger coefficients.
Xmath does not automatically perform cancellations of polynomial roots
appearing in both the numerator and the denominator of a transfer function.
If you want to cancel common roots in a transfer function, use the function
cancel( ). For state-space systems, refer to the minimal( )function.
For more information, refer to the Minimal Realizations section of
Chapter 6, State-Space Design.
To illustrate how you arrive at a particular transfer function, if you have a
system differential equation that takes the form:
··
·
·
y + 6y + 8y = 2u – u
(2-1)
Laplace-transforming equation (assuming zero initial conditions) yields:
s2Y(s) + 6sY(s) + 8Y(s) = 2sU(s) – U(s)
(2-2)
Collecting terms, you can find the transfer function from U(s) to Y(s), H(s):
Y(s)
2s – 1
s2 + 6s + 8
----------- = H(s) = -------------------------
(2-3)
U(s)
The roots of the numerator polynomial are the zeros of the transfer
function, and the roots of the denominator are its poles. In some
circumstances, you might want to construct a transfer function based on
where you know the pole and zero locations to be. For example, you can
Xmath Control Design Module
2-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
Linear System Representation
function equation using known pole, zero, and gain values:
2(s – 0.5)
H(s) = --------------------------------
(2-4)
(s + 2)(s + 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. State-space form can
be used to describe systems with multiple inputs or outputs. For more
information, refer to the State-Space System Models section.
Example 2-1
Creating Transfer Functions
The polynomials in the numerator and denominator of the transfer function
in Equation 2-3 are both in coefficients form, (described using just
coefficients, not roots). makepoly( )creates two polynomials and passes
them to the system( )function:
num3 = makepoly([2,-1],"s");
den3 = makepoly([1,6,8],"s");
H3 = system(num3,den3)
This displays as:
H3 (a transfer function) =
2s - 1
----------
s2 + 6s + 8
initial integrator outputs
0
0
Input Names
-----------
Input 1
Output Names
------------
Output 1
System is continuous
The three statements used to create the transfer function could be more
the transfer function is optional. Any variable, including the default x, can
© National Instruments Corporation
2-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Note The /operator also can be used to create systems in transfer function form, as an
alternative to using system( ).
H4 = 2*polynomial(0.5,"s")/polynomial([-2,-4],"s")
which displays as:
H4 (a transfer function) =
2(s - 0.5)
--------------
(s + 2)(s + 4)
initial integrator outputs
0
0
Input Names
-----------
Input 1
Output Names
------------
Output 1
System is continuous
In both of these cases you have created a continuous system. Systems
created in Xmath contain sample rate information as well as the numbers
representing system dynamics. However, unless a sample rate is explicitly
given as a keyword to system( ), it defaults to zero and the system is
continuous. For an illustration of how to create a discrete system, refer to
Example 2-2. The full discussion of the system( )function in the
system( ) section contains a listing of all the keywords associated with
system( ).
Xmath Control Design Module
2-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. This system model is useful for representing multi-input,
multi-output (MIMO) systems.
continuous time:
·
x = Ax + Bu
y = Cx + Du
discrete time:
xk + 1 = Axk + Buk
yk = Cxk + Duk
A straightforward mathematical transformation from the state-space form
to the transfer function form is as follows:
H(q) = C(qI – A)–1B + D
All of the forms represented in these equations can be represented using
Xmath’s system objects, as shown in Example 2-2.
Example 2-2
Creating a Discrete State-Space System
Suppose you have a system which you describe in state-space form as:
0
1
1
0
xk + 1
=
xk +
uk
–0.75 0
yk =
xk
0 1
and you know that the sample period of the system is 0.5 seconds between
samples—that is, the states and outputs are updated at every discrete
interval k, consisting in this case of 0.5 seconds.
© National Instruments Corporation
2-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
Linear System Representation
Again, you create the system using the system( )function. This time you
use the optional dtkeyword 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
0
1
D
X0
0
0
System is discrete, sampling at 0.5 seconds.
Although five lines of MathScript were used to be as explicit as possible in
creating this system, the call to system( )can encompass all of them.
Note When you create a system object, its inputs (A, B, C, D) are no longer needed.
Basic System Building Functions
The functions discussed in the following sections are available with the
general Xmath package. However, Control Design Module users will find
these functions an intrinsic part of their work, warranting this discussion.
The Xmath Help provides additional details about these functions and
examples of their use.
system( )
Sysd=system(A,B,C,D,{dt,inputNames,
outputNames,stateNames,X0})
Sys = system(num,den,{dt,inputNames,outputNames})
Sys = system(Sys,{keywords})
Xmath Control Design Module
2-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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 dtto a positive scalar value indicates that the system is
discrete, with a sampling period equal to that value. If dtis not specified,
the system is continuous, with a sampling period defaulting to zero.
Because information indicating whether the system is continuous or
discrete is encapsulated within the system object itself, Xmath does not
Systems can be recognized by Xmath’s functions as discrete or continuous
using the check( )function and handled accordingly. For more
information, refer to the Using check( ) with System Objects section.
The capability to assign a discrete sample rate does not actually discretize
a continuous-time system, however. For information on discretizing a
system, refer to the Discretizing a System section.
A shortcut for creating state-space systems with an all-zero D matrix is
to use a null-matrix specifier ([]) for the D matrix instead of entering an
appropriately sized zero matrix. This will automatically set the D matrix to
be a zero matrix with row size equal to the row size of C, and column size
equal to the column size of B.
In addition, descriptive names for the inputs and outputs of a system can
be specified as vectors of string names and assigned to the inputNames,
outputNames, and stateNameskeywords. stateNamesis valid only
when used in conjunction with a state-space system, as is the keyword X0,
When you have created a system, you can modify it by changing the values
of any of the keywords discussed in this section by calling system( )with
Examples 2-1 and 2-2 illustrate how system( )can be called to create a
transfer function and state-space system, respectively. system( )also can
be used to change the attributes of an existing system.
Note In Example 2-3, the []notation indicates that the D matrix should be an
appropriately sized (in this case, scalar) zero matrix.
© National Instruments Corporation
2-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
Linear System Representation
Example 2-3
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
0
1
D
X0
0
0
State Names
-----------
Torque Angle
Input Names
-----------
Current
Output Names
------------
Velocity
System is discrete, sampling at 0.5 seconds.
abcd( )
[A,B,C,D,X0] = abcd(Sys)
The abcd( )function extracts the component A, B, C, and D matrices
described in equations from a state-space system object as shown in the
State-Space System Models section. In addition, it returns the initial
conditions on the states if a fifth output argument is requested.
abcd( )can be called on systems in either state-space or transfer function
form. If the system is a transfer function, the conversion to state-space is
Xmath Control Design Module
2-8
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
done internally to return A, B, C, and D, though the format of the variable
Sysitself 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.
Note For the transfer function H3, an internal conversion is performed.
[A3,B3,C3,D3] = abcd(H3)?
A3 (a square matrix) =
-2 1.58114
0 -4
B3 (a column vector) =
0
2
C3 (a row vector) = -1.58114 1
D3 (a scalar) = 0
[A4,B4,C4,D4,X0] = abcd(sys4)
A4 (a square matrix) =
0
1
-0.75 0
B4 (a column vector) =
1
0
C4 (a row vector) = 0 1
D4 (a scalar) = 0
X0 (a column vector) =
0
0
© National Instruments Corporation
2-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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
conjunction with numden( )must be single-input, single-output.
As noted in the Transfer Function System Models section, common roots in
the numerator and denominator polynomials are not canceled.
Example 2-5 uses the state-space system from Example 2-2 to illustrate the
use of numden( ).
Example 2-5
Using numden( ) to Extract the Transfer Function Polynomials
sys4=system([0,1;-0.75,0],[1,0]',[0,1],0,
{dt=0.5});
[num,den] = numden(sys4)?
num (a polynomial) =
-0.75
den (a polynomial) =
(z2 + 0.75)
Because numand denare polynomial objects and not a complete system,
the discrete sampling time is not explicitly saved. You can use check( )
with the convertkeyword to map the two internal representations to each
other, as described in the Using check( ) with System Objects section.
However, notice that z was used as the polynomial variable, indicating that
these numerator and denominator polynomials were obtained from a
discrete-time system. Had the system been continuous, s would have been
used instead of z.
period( )
dt = period(Sys)
The period( )function extracts the sample period (in seconds) of a
In Example 2-5, you found the numerator and denominator polynomials
corresponding to the discrete state-space system. Example 2-6 combines
Xmath Control Design Module
2-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
check( )provides a more concise means of converting between
state-space and transfer function form, as described in the Using check( )
with System Objects section, but this example illustrates how the output of
one function can be specified directly as keyword input to another.
names( )
[outputNames,inputNames,stateNames] = names(Sys)
The names( )function extracts matrices of strings representing the input,
output, and (if the system is in state space form) state names of a system.
names( )also can be used to extract information from the PDM and
polynomial objects. More information on these functions can be found in
the MATRIXx Help.
When you create a system without specifying any names, a default set of
names are assigned to it. Unlike user-specified names, these default names
are not displayed in the Xmath Commands window. However, all, or any
subset of the names you select to store with the system still can be extracted
using names( )as shown in Example 2-7.
Example 2-7
Using names( ) to Extract the Variable Names Associated with a System
H3 = system(makepoly([2,-1],"s"),
makepoly([1,6,8],"s"));
[outputNames, inputNames] = names(H3)
outputNames (a string) = Output 1
inputNames (a string) = Input 1
sys5=system([0,1;-0.75,0],[1,0]',[0,1],0,{dt=0.5,
inputNames = "Current",
© National Instruments Corporation
2-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
You can index into a dynamic system to create a new dynamic system
which has a subset of the original inputs and outputs:
Sys = Sys1(i,j) is defined to be a system such that y = y1(i) and u = u1(j).
i and j can both be vectors as well, in which case multiple inputs and
outputs will be extracted.
The previous definition of indexing was designed with the traditional
definition of a transfer function in mind.
y(q) = Sys(q) × u(q)
Using check( ) with System Objects
Several common attributes of systems can be easily determined using
Xmath’s ability to distinguish between object types and characteristics.
You can use the check( )function with systems, as shown in
Example 2-8, to determine whether a system is in transfer function or
state-space form, discrete, continuous, or stable. In addition, you can use
check( )with the convertkeyword to change a system’s representation
between SISO state-space and transfer-function forms.
Example 2-8
Using check( ) with a System
a = [1.875,0;0,-0.26];
b = [1;0];
c = [0.5,1];
d = 0;
sys = system(a,b,c,d, {dt = 0.001});
Because this system is discrete and has a pole where magnitude exceeds 1,
it is not stable.
Xmath Control Design Module
2-12
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. Therefore, you often choose to discretize these models for use
with a digital controller. A number of mathematical methods have been
developed to approximate the behavior of a continuous system in a
discrete-time representation with an appropriately fast sampling rate.
Xmath provides two functions, discretize( )and
makecontinuous( ), which encompass a range of these techniques.
discretize( )converts a system from its representation as a continuous
function in the s-domain to a discrete-time z-domain function.
makecontinuous( )does the reverse, transforming a discrete system to
its continuous form.
discretize( )
SysD = discretize(Sys,{dt,exponential,forward,backward,
tustins,ztransform,polezero,firstorder})
The discretize( )function has a number of keywords that correspond
to the different methods of continuous-to-discrete conversion that are
implemented within Xmath. The sampling interval (in seconds) for the
is specified, a default of 0.5 seconds is used. The default discretization
method used is the exponential (step-invariant) transform. The different
© National Instruments Corporation
2-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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).
To convert the system description from a continuous differential equation
to a discrete difference equation, you approximate the value of the
derivative in the continuous equation over each dtseconds of time, then
find the area of the geometric region having width dtand height equal to
the derivative. You can do this in a number of ways, as discussed in
[FPW90].
For the forward rectangular method, you assume the incremental area term
between sampling times k * dt and (k + 1) * dt to be a rectangle having
width dt and height equal to the integral form of the differential equation at
time (k + 1) * dt. In essence, you get your amplitude estimate for each
rectangle by looking forward, hence the name. The backward rectangular
method arises similarly, except that you get the rectangle’s height by
looking backward and taking the value of the integral at k * dt. The forward
rectangular approach tends to overestimate the incremental area somewhat
and the backward approach tends to underestimate it (though with a
sufficiently small sampling interval, this may not pose a large problem).
The trapezoid rule strikes a balance between these two methods by taking
and using that value as the incremental area in approximating the difference
equation.
These approaches can be summarized as substitutions between the
continuous-time Laplace-transform operators and the discrete z-transform
operator z as shown in Table 2-1.
Table 2-2. Mapping Methods for discretize( )
Method of Approximation
Continuous to Discrete
z – 1
dt
Forward rectangular rule:
----------
s →
Xmath Control Design Module
2-14
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
Linear System Representation
Table 2-2. Mapping Methods for discretize( ) (Continued)
Method of Approximation
Continuous to Discrete
z – 1
zdt
Backward rectangular rule:
Keyword: backward
----------
s →
2(z – 1)
Tustin’s rule:
Keyword: tustins
---------------------
s →
dt(z + 1)
Pole-Zero Matching: polezero
The pole-zero matching method of discretizing a continuous system
follows from the relation between the continuous s and discrete z frequency
domains:
z = esT
where T is the sampling interval to be used for the discrete system.
Continuous-time poles and finite zeros are mapped to the z-plane using this
relation. Zeros at infinity are mapped into z = 0, where they do not affect
the frequency response.
After the poles and zeros have been mapped, the algorithm tries to
make sure the system gains are equivalent at some critical frequency.
If the systems have no poles or zeros at DC(s = 0, z = 1), the discrete-time
gain is selected such that the system gains match at DC. Alternatively,
if the systems have no poles or zeros at the Nyquist frequency
(s = p * j/T, z = –1), the gains are equalized at that frequency. In the
event that neither of these gains can be matched, no gain is chosen.
Z-Transform: ztransform
This method is a direct Z-transform of the continuous-time transfer
function, which corresponds to the Z-transform of the impulse response of
the system. If ztransformis used, you will match the impulse responses
of the continuous and discrete systems. The responses may differ slightly
due to round off error.
Hold Equivalence Methods: exponential and
firstorder
The discretization methods for exponentialand firstorderboth rely
on the approximation that the discrete-time response can be represented as
a hold on the sampled values of the continuous-time response.
© National Instruments Corporation
2-15
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 2
Linear System Representation
The exponentialkeyword assumes that the response value between
samples is constant and can, therefore, be represented by a zero-order hold
polynomial. When exponentialis 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.
The firstorderkeyword assumes extrapolation between samples
(connecting sample to sample in a straight line). If firstorderis
specified, the continuous-time ramp response is discretized using the
Z-transform and then the result is divided by the Z-transform of a ramp
z * dt/(z – 1)2 to produce the desired transfer function.
In each of these cases, the appropriate response (impulse, step, or ramp)
round off error.
While no one method of discretization will always perform best for all
systems and all sampling times, it is often a good idea to compare the
frequency response resulting from different discretized models to the
continuous response. Example 2-9 applies the forward, backward, tustins,
exponential, and matched pole-zero discretization methods.
Example 2-9
A Comparison of Several Discretization Methods
H = system(0.5*polynomial([-0.36]),
makepoly([1,2.79,2.74,1.11,0.16]));
Create a logspaced vector for the frequency range of the response:
F = logspace(.001,5,200);
Perform the discretization using the different algorithms:
Hd_f = discretize(H,0.1,{forward});
Hd_b = discretize(H,0.1,{backward});
Hd_t = discretize(H,0.1,{tustins});
Hd_z = discretize(H,0.1,{polezero});
Hd_e = discretize(H,0.1,{exponential});
Now you can calculate the magnitude response as a function of frequency,
gainc = 20*log10(abs(freq(H,F)));
gain_f = 20*log10(abs(freq(Hd_f,F)));
gain_b = 20*log10(abs(freq(Hd_b,F)));
gain_t = 20*log10(abs(freq(Hd_t,F)));
Xmath Control Design Module
2-16
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. Comparison of Different Frequency Response Techniques
Although most of the discretizations used would give acceptable
approximations to the continuous-time response, notice that most of them
diverge greatly at higher frequencies. You may find it illustrative to run this
example with larger and smaller sampling intervals to see how the choice
of sampling rate, as well as the choice of method, affects the accuracy of
the discretized frequency response.
makecontinuous( )
Sys=makecontinuous(SysD,{exponential, forward,
backward,tustins, ztransform})
© National Instruments Corporation
2-17
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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
makecontinuous( )accepts an input system in any form, it returns the
continuous-time system as a state-space system.
The forward, backward, and Tustin methods for mapping from the s-plane
to the z-plane can be easily reversed using the equivalencies shown in
Table 2-3.
Table 2-3. Mapping Methods for makecontinuous( )
Method of Approximation
Discrete to Continuous
z → 1 + s(dt)
Forward rectangular rule:
Keyword: forward
1
Backward rectangular rule:
Keyword: backward
--------------------
z →
1 – s(dt)
1 + s(dt) ⁄ 2
Tustin’s rule:
Keyword: tustins
----------------------------
z →
1 – s(dt) ⁄ 2
Discrete-to-continuous algorithms using matrix logarithms (to reverse the
exponential operations involved in doing the z-transform for the impulse
invariant zero-order hold) are available for the exponential
(step-invariant) transformation and the ztransform(impulse-invariant)
methods. A limitation of these methods, however, is that they will not return
a meaningful continuous equivalent to a discrete system that has pure
delays (1/z terms), because the logarithm of zero is infinite.
Example 2-10 Verifying a Discretization Using makecontinuous( )
Create a system:
H = 0.5*polynomial([-0.36])/...
polynomial([-1,-1,-0.395+0.06305*jay,
-0.395-0.06305*jay]);
Form the discrete equivalent using the forward approximation:
Hd_f = discretize(H,0.1, {forward});
Xmath Control Design Module
2-18
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
© National Instruments Corporation
2-19
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. A number of simple
types of connections have been implemented as overloaded operators on
systems, while more complex connections are available through
specialized functions.
Linear System Interconnection Operators
Overloaded operators provide a quick way to perform different types of
basic connections between systems. Table 3-1 illustrates these operations
on a pair of systems Sys1 and Sys2 with outputs y1 and y2 and inputs u1
and u2, respectively.
Table 3-1. Summary of Interconnection Operators
Diagram
Description
Sys = Sys1 + Sys2
Parallel connection where y = y1 + y2. The inputs are
tied together where u = u1 = u2.
u1
y1
Sys1
u
y
+
Sys2
u2
y2
Sys = Sys1 – Sys2
Parallel connection where y = y1 + y2. In the unary
case, Sys=–Sys1 where y=–y1. The inputs are tied
together where u = u1 = u2.
u1
y1
Sys1
u
y
–
Sys2
u2
y2
© National Instruments Corporation
3-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
Building System Connections
Table 3-1. Summary of Interconnection Operators (Continued)
Diagram
Sys = Sys2 * Sys1
Description
Cascade connection of Sys1 and Sys2 where the output
of Sysis y2 and the input is u1.
u1
y1
u2
y2
y
u
Sys2
Sys1
Sys = Sys1/Sys2
Cascade connection of Sys1and inverted Sys2where
Sys = Sys1 * inv(Sys2), u = u2, and y = y1.
u2
y2
u1
y1
y
y
u
Sys1
inv(Sys2)
Sys = Sys1\Sys2
Cascade connection of inverted Sys1 and Sys2 where
Sys = inv(Sys1) * Sys2, u = u2, and y = y1.
u2
y2
u1
y1
u
inv(Sys1)
Sys2
Sys = [Sys1;Sys2
Parallel connection where y = [y1;y2]. The inputs are
tied together where u = u1 = u.
u1
y1
Sys1
Sys2
u
y
u2
y2
Sys = [Sys1,Sys2
Parallel connection where y = y1 + y2 and
u =[u1;u2].
u1
y1
Sys1
y
+
u
Sys2
u2
y2
Sys = Sys1'
If Sys1is in state-space form and comprises the matrices
(A1,B1,C1,D1), Syscomprises (A1',C1',B1',D1').
If Sys1 is a transfer function, it is converted internally
to state-space form.
Xmath Control Design Module
3-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
Building System Connections
Table 3-1. Summary of Interconnection Operators (Continued)
Diagram
Description
Sys = adj[Sys1]
If Sys1is in state-space form and comprises the matrices
(A1,B1,C1,D1), Sysis the adjoint system and comprises
(–A1',C1',B1',D1'). If Sys1is a transfer function, it is
converted internally to state-space form.
p1/p2
Alternate method to create a system, where p1and p2
are the numerator and denominator polynomials,
respectively; does not allow the use of keywords.
Sys = inv(Sys1)
The inverse (pseudoinverse) of a system can be found
using inv(Sys1). If Sys1 is a transfer function,
inv(Sys1)is the reciprocal of the transfer function.
If Sys1 is a state-space system (A1,B1,C1,D1), then
Sys= system(A,B,C,D)where A,B,C,Dare defined
as follows:
D = pinv(D1)
A = A1-B1*D*C1
B = B1*D
C = -D*C1
Dynamic systems also can be flexibly combined with scalars and
compatibly sized matrices using the operators in Table 3-1. A compatibly
sized matrix is one having the same dimensions as the dynamic system’s
D matrix (row size equal to the number of outputs and column size equal
to the number of inputs).
Operations performed with a dynamic system and a matrix M as the
system([],[],[],M).
The *operator can be used with a system and a PDM to find the time
response of the system to the general input data stored in the PDM. For
a detailed description of time simulation in Xmath, refer to the General
Time-Domain Simulation section of Chapter 4, System Analysis.
© National Instruments Corporation
3-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
•
If the systems to be connected are in transfer-function form, they must
be expressed in the same dependent variable.
In describing the algorithms used in these connection functions, we will
often refer to the component matrices of a state-space system Sys1 as A1,
B1, C1, and D1.
afeedback( )
loop, and obtains a single system representation for the complex loop.
Sysis organized as shown in Figure 3-1. Additional external inputs to the
feedback path are included with outputs from the feedforward path.
Figure 3-1 illustrates that outputs of the feedback path system are included
with forward path outputs. For an example of how to use afeedback( ),
refer to Example 3-1.
e1
u1
y1
Sys1
Sys2
+
–
+
+
y2
u2
e2
Sys
Figure 3-1. afeedback System Configuration
Xmath Control Design Module
3-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
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.
•
When only one system is specified, it must be square (it must have an
equal number of inputs and outputs).
Example 3-1
Using afeedback( ) to Connect Two Systems
Sys1 = system([.5,1;0,2],[1,0]',[0,1],0);
Sys2 = system([1,-.2;1,0],[1,0]',[1,1],0);
saf = afeedback(Sys1,Sys2);
Algorithm
If only one system input (Sys1) is provided to afeedback( ), the second
input (Sys2) defaults to a zero-state system with unity gain. This is
analogous to a state-space system with NULLvalues for the A, B, and C
matrices, and with an identity matrix for D. Notice that you use the Xmath
definition of a non-square identity matrix. In this case, the row dimension
of D equals the number of inputs to Sys1, and the column dimension equals
the number of outputs of Sys1. In the following discussion, you denote the
state-space matrices of Sys1by A1, B1, C1, and D1, and you follow the same
convention for Sys2.
The two systems are first internally converted to a state-space form, if
necessary, and subdivided into the A, B, C, and D state-space matrices.
Scaling matrices S1 and S2 are computed for Sys1 and Sys2as follows:
S1 = I + D1D2
S2 = I + D2D1
Additionally, you define:
B
1s = B1/S2 and D1s = D1/S2
B2s = B2/S1 and D2s = D2/S1
Matrix right-division problems must be well-posed, with the scaling
matrices S1 and S2 nonsingular. afeedback( )displays an error message
© National Instruments Corporation
3-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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 Systo 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. In the case of a constant-gain feedback, A, B, C, and D are
computed using only the matrix partitions shown in bold type.
The initial conditions for the systems are appended to each other
columnwise.
A1 – B1sD2C1
–B1sC2
A =
B =
C =
D =
B2sC1
A2 – B2sD1C2
B1s –B1sD2
B2sD1 B2s
S1 ⁄ C1 –D1sC2
D2sC1 S2C2
D1s –D1D2s
D2D1s D2s
afeedback( )cannot be used with improper transfer functions—systems
having more zeros than poles—because this algorithm is strictly
state-space.
append( )
Sys = append(Sys1,Sys2)
The append( )function appends two dynamic systems in a form suitable
for use with the connect( )function (refer to the connect( ) section). Sys
Xmath Control Design Module
3-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
Building System Connections
is created by appending the inputs, outputs, and states of Sys1and 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 Sysis represented in state-space form.
The output is a dynamic system in block form as shown in Figure 3-2.
u1
y1
Sys1
u2
y2
Sys2
Sys
Figure 3-2. Output of a Dynamic System
For an example of how to use append( ), refer to Example 3-2.
Using append( )
Example 3-2
s3 = system(makepoly([1,2]), makepoly([1,3,5,0]));
s4 = system(makepoly([1]), makepoly([1,4,4]));
sap = append(s3,s4);
In the following discussion, the component state-space matrices of Sys1are
denoted by A1, B1, C1, and D1, and you follow the same convention for Sys2.
The algorithm for append( )is done strictly using the state-space
representations of Sys1 and Sys2. For this reason, Sys1 and Sys2 cannot
be improper transfer functions (transfer functions having more zeros than
poles). The component A, B, C, and D matrices of Sys1 and Sys2 are
extracted using the abcd( )function. The A, B, C, and D matrices
comprising Sysare obtained as shown in the following example, where the
zero matrix elements span as many rows and columns as necessary.
© National Instruments Corporation
3-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
Building System Connections
A1
0
B1
0
A =
B =
0 A2
0 B2
C1
0
D1
0
C =
D =
0 C2
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. This provides two basic capabilities:
•
•
constant gain feedback
general input–output interconnection
In its simplest form, connect( )can be used to wrap constant gain
also can be used to specify which outputs are fed back to the input of the
system. By specifying the optional keywords Mand N, you also can specify
input and output gains.
General input–output interconnection is applicable to the block form
system provided by append( ), as described in the append( ) section.
Parameters used in the connectcommand are illustrated in Figure 3-3.
Notice that feedback is defined with a positive sign.
u
u1
y1
y
M
Sys1
K
N
+
+
uk
Sys
Figure 3-3. Parameters Used with the connect Command
Xmath Control Design Module
3-8
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
not allowed.
The number of outputs in the combined system is the sum of the number of
outputs from the two systems you are appending. For an example of how to
use gains for the input, the output, and the fed-back data, refer to
Example 3-3.
Example 3-3
Using connect( ) to Perform a General Output-Input Connection
tfsys=system(makepoly([1,2]),makepoly([1,3,0]))
tfsys (a transfer function) =
x + 2
-------
2
x + 3x
initial integrator outputs
0
0
Input Names
-----------
Input 1
Output Names
------------
Output 1
System is continuous
connect(tfsys,0.12,2,1.5)
ans (a state space system) =
A
-3 1
-0.12 0.12
B
0
© National Instruments Corporation
3-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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 = A1x + B1u1
y1 = C1x + D1u1
Combining these equations with the equation for the positive feedback
input term:
u1 = Ky1 + 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. If you do not specify any values for the gain matrices,
K defaults to zero (no feedback) and M and N default to appropriately-sized
identity matrices (unity gain on the input and output).
–1
–1
·
x = (A1 + B1(I – KD1) KC1)x + B1(I – KD1) Mu
y = N(I – KD1)–1C1x + ND1(I – KD1)–1Mu
This algorithm assumes that the closed-loop system is well posed to ensure
that Syswill be proper. The (I – KD1) term must be invertible, and a
warning appears if the condition estimate of the term (refer to rcond) is
less than eps.
Xmath Control Design Module
3-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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).
e1
u1
y1
Sys1
+
–
Sys2
Sys
Figure 3-4. Feedback System Configuration
For an example of how to implement unity gain feedback using the
feedback( )function, refer to Example 3-4.
Example 3-4
Implementing Unity Gain Feedback Using feedback( )
tfsys2 = polynomial(-1)/polynomial([-2,-3]);
feedback(tfsys2) # Note that this implements
# negative unity gain feedback
ans (a state space system) =
A
-2 1
1 -4
B
0
1
C
© National Instruments Corporation
3-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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 NULLmatrices for A2,
B2, and C2 and an identity matrix for D2 so that unity gain feedback is
implemented.
From Figure 3-4 and the state-space definitions of the systems, you derive
the following equations:
·
x1
=
=
=
=
=
A1x1 + Be1
C1x1 + D1e1
A2x2 + B2y1
C2x2 + D2y1
u1 – y2
y1
·
x2
y2
e1
Xmath Control Design Module
3-12
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 3
Building System Connections
The single system resulting from the feedback combination of Sys1and
Sys2has u1 as its input, y1 as its output, and a state vector consisting of the
appended states of Sys1and Sys2. Using these five equations to find the
state-space dynamics of the complete system results in the overall system
description.
·
x1
·
x =
•
·
x2
A1 – B1(I + D2D1)–1D2C1
–B1(I + D2D1)–1C2
A2 – B2(I + D1D2)–1D1C2
=
•
B2(I + D1D2)–1C1
x1
x2
B1(I + D2D1)–1
B2(I + D1D2)–1D1
+
u1
x1
(I + D1D2)–1C1 –D1(I + D2D1)–1C2
y1 =
+
x2
D1(I + D2D1)–1u1
This algorithm assumes that the closed-loop system is well constructed
(the (I + D2D1) and (I + D1D2) terms must be invertible). This condition
ensures that the output system Syswill be proper.
© National Instruments Corporation
3-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Time-Domain Solution of System Equations
Given the state-space equations:
·
x = Ax + Bu
y = Cx + Du
you obtain:
t
x(t) = eAtx0 + eAτBu(t – τ)dτ
∫
0
letting x0 denote any initial conditions on the system states. The integral
term in the preceding equation defines a convolution integral. Using *to
represent the convolution operator, the time-domain system output for all
time t ≥ 0 is:
h(t) = CeAtB + Dδ(t)
*
(4-1)
The response Y(s) of the system (with zero initial conditions) to a unit
impulse input δ(t) is H(s), the transfer function representation from the
Transfer Function System Models section of Chapter 2, Linear System
of H(s), the impulse response.
© National Instruments Corporation
4-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Using *to indicate discrete convolution, you can express the time domain
output as a function of the discrete impulse response:
yk
hk
=
=
CAkx0 + (hk u )
*
k
(4-2)
CAk – 1
B
(k > 0)
System Stability: Poles and Zeros
After you have general expressions for the response of a system over time,
according to Equations 4-1 and 4-2, you can assess the stability of the
system. For the purposes of system analysis within Xmath, you define a
stable system as one where output does not grow without bound for any
bounded input or initial condition. A necessary and sufficient condition for
this type of bounded-input bounded-output (BIBO) stability is:
∞
h(t) < M < ∞
∫
0
Continuous systems are BIBO stable if and only if all poles of the system
are in the left half of the complex plane; discrete systems are BIBO stable
if and only if all poles are within the unit circle in the complex plane.
For a coprime transfer function H(q) (one having no root cancellations
between the numerator and the denominator), the poles are the roots of the
denominator of H(q). H(q) is infinite at these values. Values of q for which
the numerator of H(q) is zero are termed the zeros of the system.
The poles of a system in transfer-function form are identical to the
eigenvalues of the A matrix in that system’s equivalent state-space
representation.
For systems in transfer-function form, zeros are easily defined as the
polynomial roots of the numerator. You define the system matrix for a
state-space system as
λI – A B
S(λ) =
–C
D
Xmath Control Design Module
4-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. If the response becomes larger over time, the pole is said to be
unstable. If the response remains unchanged over time, you describe the
pole that causes it as neutrally stable. All the closed-loop poles of a system
must be stable to describe the system as stable.
poles( )
p = poles(Sys)
If the input system Sysis in transfer-function form, poles( )obtains the
poles from the roots of the transfer function’s denominator (which are
automatically stored if the system is in zero-pole format or if the roots have
been previously calculated). If Sysis in state-space form, the poles are
computed as the eigenvalues of the A matrix. To see how to use poles( )
with a system in transfer function form, refer to Example 4-1.
Example 4-1
Using poles( ) with a System in Transfer Function Form
H = 0.5*polynomial([-0.36])/...
makepoly([1,2.79,2.74,1.11,0.16]);
poles(H)
ans (a column vector) =
-0.395 + 0.0630476 j
-0.395 - 0.0630476 j
-1
-1
zeros( )
[z,k] = zeros(Sys)
The zeros( )function finds the invariant zeros, the values of λ at which
R(λ) = 0 and S(λ) lose rank, and gain is returned only for SISO systems
© National Instruments Corporation
4-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
directly from the roots of the transfer function numerator. If Sysis in
state-space form, the definition of its zeros arises from the system matrix,
λI – A B
S(λ) =
(4-3)
–C
D
and its MIMO transfer function:
R(λ) = C(λI – A)–1B + 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.
of the system. For degenerate cases in which the normal rank of S(λ) is less
than n + r, the zeros are defined analogously.
If a system is minimal (that is, no other system with lower order and the
same R(λ) exists), these invariant zeros are termed transmission zeros.
When the matrix in Equation 4-4 loses rank for some value λ = λ0, there
exists a vector [x0' u0']' of initial states and inputs such that:
x0
u0
λ0I – A B
–C
= 0
D
Thus, there exists an initial state x0 such that the output y is zero for all
values of the input function defined over time t as u0eλt. Such zeros (λ0)
transmission of the system input to the output.
Note zeros = system zeros = {invariant zeros} ∩ {transmission zeros}.
For an example using zeros( )with a state-space system, refer to
Example 4-2. For more details on this topic, refer to [Kai80] and [DeS74].
Example 4-2
Using zeros( ) with a State-Space System
Sys=system([-2.3,0.01,5.1;0,-0.35,-2;0,2,-.35],
[1,.25,.25]',[1.34,0,0],0);
[z,k] = zeros(Sys)
Xmath Control Design Module
4-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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].
Note zeros( )also can be used as a matrix building function when used with scalar or
matrix input. For more details on this usage, refer to the MATRIXx Help.
Partial Fraction Expansion
By inverse Laplace or z-transforming a transfer function, you can identify
the impulse response based on knowledge of the system pole and zero
locations. The most convenient form to use in doing this is the partial
fraction expansion of the transfer function. Each term of the partial fraction
expansion has a constant numerator—the residue—and a pole term
denominator, as shown in the following equation, where p2 is a repeated
pole, and p4 and p5 are a conjugate pair:
r1
H(q) = C + -------------- + --------------------- + -------------- + -------------------------------------- + … + --------------
q + p1 q + p3 (q + p4)(q + p5) q + pn
r2
r3
α4q + α5
rn
(q + p2)2
Each pk represents a pole of the system, and the corresponding rk is the
residue at that pole. If pk is a repeated pole, it has M residues, where M is
the multiplicity of the pole. Complex pole pairs have complex residue pairs.
If the transfer function contains a constant (or feedthrough) term, this term
is represented by the scalar value C in the preceding equation. The values
of the residues give the magnitude of the response from the inverse
transform of the respective partial fraction terms. For an example of
dynamic response with partial fraction expansion, refer to Example 4-3.
[Oga70] provides a good reference on partial-fraction expansion for
different orders of complex and real poles. [ChB84] contains a thorough
mathematical treatment of residues.
© National Instruments Corporation
4-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
Example 4-3
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) = ----------------------------------------
s2(s + 2)(s + 10)
You close a unity feedback loop around this system, as shown in Figure 4-1.
+
–
u(s)
V(s)
Y(s)
G(s)
Gcl(s)
Figure 4-1. Constructing theClosed-Loop System Gcl(s) from the Open-Loop System
G(s), with Input U(s) and Output Y(s)
You can derive the expression for the closed-loop transfer function Gcl(s):
V(s) = U(s) – Y(s)
Y(s)
G(s)
⎛
⎝
⎞
⎠
Y(s) = G(s)V(s)
Gcl(s) = ----------- = --------------------
U(s) 1 + G(s)
Calculate the closed-loop transfer function.
Note You convert the state-space system returned by feedback( )to a transfer function
using check( ).
sys = polynomial(-0.5)/polynomial([0,0,-2,-10]);
syscl=feedback(sys);
[,syscl] = check(syscl,{tf, convert})
syscl (a transfer function) =
(s + 0.5)
-------------------------------------------------
2
(s + 1.95266)(s + 10.0118)(s + 0.0354992s + 0.02...
initial integrator outputs
0
0
Xmath Control Design Module
4-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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
-10.0118
| Order 1 -0.0478224
| Order 1 0.0118134
residue( )returns a PDM with the poles as the domain elements, and the
associated dependent matrices being the residue at each pole. It also can be
expressed in the following form:
–0.0478
(s + 1.95) (s + 10.012)
0.0118
Gcl(s) = ----------------------- + ----------------------------- +
0.036(s + 0.01775) + 0.16065(0.1589)
----------------------------------------------------------------------------------------------
(s + 0.01775)2 + (0.1589)2
Using a table of inverse Laplace transforms to convert this expression to the
transient time response rather than a complex frequency response, you can
rewrite the time response G(t) as:
G(t) = –0.0478e–1.95t + 0.0118e–10.012t
–0.018t( 0.036 cos(0.1589t) + 0.160 sin(0.1589t))
+
e
Notice from this example that because all the poles are in the left half plane,
the response each contributes is an exponential which decays with time, so
this closed-loop system is stable.
© National Instruments Corporation
4-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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
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.
Notice that this response actually takes quite a while to die out because of
the small time constants, which correspond to small pole values, in the
exponential terms. This is why poles with a small magnitude are frequently
called “slow” poles, whereas poles with a large magnitude contribute a
response which decays quickly and thus are called “fast” poles.
residue( )
[rp,C] = residue(sys,pls,ordr,{isVector,tol})
The residue( )function calculates the nth-order residue of a
transfer-function form system at any of its poles, including Infinity.
It returns a PDM rpwhere domain contains the pole locations and where
dependent matrices contain the residues corresponding to each pole.
plsand ordrare optional inputs allowing you to specify the pole values
Xmath Control Design Module
4-8
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
and orders for which the residue(s) should be found. If a user-specified
value for plsis 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. Ccontains 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.36]),
polynomial([-1,-1,-0.395+0.06305*jay,
-0.395-0.06305*jay]));
Rp=residue(G)
Rp (a pdm) =
Poles |
-------------------+-----------------------------
-0.395 - 0.06305 j | Order 1 0.738493 - 0.2277 j
| 2
0
-------------------+-----------------------------
-0.395 + 0.06305 j | Order 1 0.738493 + 0.2277 j
| 2
-------------------+-----------------------------
-1 | Order 1 -1.47699
| 2 -0.864864
0
-------------------+-----------------------------
combinepf( )
Sys = combinepf(Rp,C,{var})
combinepf( )reverses the operation performed by residue( ),
combining partial fractions into a single transfer function. It expects a PDM
of the form shown in Example 4-4 as input.
Use combinepf( )to convert partial fractions to a transfer function.
Using the variable Rpyou obtained in Example 4-4:
G2=combinepf(Rp, {var = "s"})
G2 (a transfer function) =
0.5s + 0.18
---------------------------
2 2
© National Instruments Corporation
4-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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 G2matches the system Gwhere residues were computed in Example 4-4.
General Time-Domain Simulation
When modeling a dynamic system and trying to determine its response to
the input values it is likely to receive in use, you generally will want to
simulate system behavior with more general input signals than the zero or
step inputs used in initial( ), impulse( ), and step( ). To do this,
use the *operator between a dynamic system object and a PDM containing
the input data you want to use in the simulation.
Borrowing from the standard frequency response notation for a system
where:
y(s) = H(s) × u(s)
Xmath defines the operation system*PDM as a time domain simulation.
Thus for any dynamic system Sys(continuous or discrete) and for a PDM
u representing the external stimulus as a function of time, the operation
y = Sys× u creates a PDM y which contains the outputs of the system as
a function of time.
For a dynamic system with ny outputs and nu inputs, the input vector is
defined to be ny × 1 and the output vector is ny × 1. Thus the input PDM
u should be m × 1 × p, where p is the number of time points in u.
Xmath Control Design Module
4-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
The input PDM must have a regular domain—that is, the interval between
each domain value and the one succeeding it must be the same over all
points in the domain. If the system is discrete, the domain intervals must be
equal to the system’s sampling period. If the system is continuous, it is
discretized using the exponential (zero-order hold) method, with the
Note For accurate results, you need to make sure this sampling interval is small enough
that discretization effects are negligible.
The next step is to create a general signal and store it as a PDM where
domain is time as shown in Example 4-5. Because you are using a SISO
system, this input is a single-channel PDM.
Example 4-5
Performing a General Time-Domain Simulation
t = 0:0.1:15;
osig = ones(1,30);
sig = [0*osig,0.5*osig,osig,0.5*osig,0*ones(1,31)];
U = pdm(sig,t);
Create the system:
Sys = system([-2.3,0.01,5.1;0,-0.35,-2;0,2,-.35],
[1,.25,.25]',[1.34,0,0],0);
and perform the simulation:
Y = Sys*U;
To see how well the system tracks the input signal, plot the input,
as follows, and the system’s response, shown in Figure 4-3.
plot ([U,Y], {legend = ["Input Signal",
"System Response"],line_color = "black",
xlab = "Time (sec)", ylab = "Amplitude"})
© National Instruments Corporation
4-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Example 4-6
General Time-Domain Simulation for a MIMO System
sys = system([0,1,0;0,0,1;-2,-4,-3],
[0,1;1,0;-1,1],[0,1,-1;1,2,1],[]);
u = pdm([sin(-3*pi:0.01:3*pi);
cos(-3*pi:0.01:3*pi)],{rows=2,columns=1});
y = sys*u;
Xmath Control Design Module
4-12
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
For a discrete system, an impulse can be thought of as a physical pulse
all other times.
The Laplace transform of the continuous-time impulse—often referred to
as δ(t)—is 1. Thus, the Laplace transform of a output of a system to a unit
impulse is merely its transfer function H(s), as discussed in the
Time-Domain Solution of System Equations section.
A similar definition, using the z-transform, can be made for the
discrete-time impulse response. However, the values of the impulse
response of a discrete system also have the property that they define the
Markov parameters for that system. Based on the state-space representation
of the system, these parameters are defined as the values
hi = {CAi – 1B,i = 1,2,…}
These parameters are uniquely determined by the transfer function of the
system [Kai80]:
∞
H(z) = C(zI – A)–1B =
hiz–i
∑
i = 1
and they also are the terms of the discrete impulse response.
impulse( )
y = impulse(Sys,t)
The impulse( )function computes the impulse response of a dynamic
system. The time vector, t, is an optional input. If not specified, a default
time range will be computed using deftimerange( ). Refer to the
response is calculated at each point in the time vector. For a discrete
system, the first n Markov parameters are returned, where n is the length
of the time vector (which must be regularly spaced).
© National Instruments Corporation
4-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
Note 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.
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. The impulse response is then a time-domain simulation of
the system using an appropriately-sized zero input.
The output yis a PDM where domain is the time vector t. Each dependent
columns as there are inputs of Sys. Thus the (i,j,k) element of yis the
impulse response at output ifrom input jat time k. In Figure 4-4, where
all the poles of this continuous system are stable (in the left half of the
complex plane), the impulse response eventually dies out to zero. For an
example of a 15-second impulse response of a stable state-space system,
refer to Example 4-7.
Example 4-7
15-second Impulse Response of a Stable State-Space System
Sys = system([-2.3,0.01,5.1;0,-0.35,-2;0,2,-.35],
[1,.25,.25]',[1.34,0,0],0);
Yt = impulse(Sys,0:0.1:15);
plot (Yt, {xlab = "Time (sec)",
ylab = "Amplitude"})
Xmath Control Design Module
4-14
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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( ). For continuous-time systems, the poles are scaled by a factor
of 1/2π (to convert from radians) and the time constant (in seconds) is
obtained as the reciprocal of four times the value of the pole with the
maximum absolute value (the “fastest” pole). For discrete-time systems,
the logarithm of the poles is taken and scaled by the sampling interval. The
sampling interval is automatically used as the step size for the tvectime
vector. If all system poles are integrators, the step size defaults to 0.01.
© National Instruments Corporation
4-15
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
The maximum time, Tmax, is computed as follows, with vPdenoting the
vector of scaled poles and dtthe 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. As a result, if Syshas zeros that are more than
a decade beyond its maximum or minimum poles, the effects of these zeros
may not be apparent in a time response calculated using tvec. You should
supply your own time vector to impulse( ), initial( ), and step( )
in these cases.
System Response to Initial Conditions
It is often assumed that the states of a system have zero initial conditions,
and the X0field of a state-space system object correspondingly defaults to
zero. In many cases, however, you need to examine the system response to
a given set of nonzero initial conditions; a common system design goal is
that this response become zero (or negligibly small) as quickly as possible.
The Xmath function initial( )allows you to do this. You also can use
superposition to calculate forced initial condition response.
initial( )
Y = initial(Sys,T,{X0})
The initial( )function computes the unforced response of a system
from its initial conditions. By default it uses the initial conditions stored
with the input system itself, but an alternate set of initial conditions can be
specified as the keyword X0. A time vector (with spacing equal to the
sampling period of the system if it is discrete) also can be specified, or
initial( )can compute a default time vector internally using
Xmath Control Design Module
4-16
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 4
System Analysis
point in the time vector. The output Yis 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,.25,.25]',[1.34,0,0],0);
ic1 = initial(Sys, 0:.1:15, {X0 = [1,0,0]});
ic2 = initial(Sys, 0:.1:15, {X0 = [0,1,0]});
ic3 = initial(Sys, 0:.1:15, {X0 = [0,0,1]});
plot ([ic1,ic2,ic3], {xlab = "Time (sec)",
legend = ["State 1", "State 2", "State 3"],
ylab = "Amplitude"})
Figure 4-5. 15-Second System Response to Unity Nonzero Conditions
at Each of the States
In Figure 4-5, notice that the value of the second state has the highest
maximum value and takes the longest to “zero out.”
© National Instruments Corporation
4-17
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The performance measures associated
with the step response are as follows:
•
Delay time (td)—The time required for the response to reach half its
final value.
•
•
•
•
Rise time (tr)—The time required for the response to rise from 10%
of its final value to 90% of its final value.
Peak time (tp)—The time required for the response to reach the peak
value of its first overshoot.
unity, expressed as a percent.
Settling time (ts)—The time required for the response to reach 5% of
its final value.
These performance measures are obtained easily with a few lines of
MathScript, as demonstrated in Example 4-9. For a plot of these
performance measures, refer to Figure 4-6.
step( )
Y = step(Sys,T)
The step( )function computes the unit step response of a dynamic
system over a time period which can be specified with the optional time
vector T. If Tis not specified, step( )computes a default time vector
using deftimerange( ).
The output, Y, is a PDM where domain is the time vector and dependent
matrices have row size equal to the number of inputs and column size equal
to the number of outputs.
Example 4-9
Performance Measurements for a Step Response
Y = step(Sys, 0:.1:15);
plot (Y,{x_lab = "Time (sec)",
ylab = "Amplitude", xinc=1, yinc=.1})
Xmath Control Design Module
4-18
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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
Yand substantiate your estimates.
First, you find the final value of the response:
Yf = makematrix(Y(151));
Get indices of all values >half the final value:
gt_half = find(Y > 0.5*Yf);
Time corresponding to first index in gt_half:
td = domain(Y(gt_half(1,1)))
td (a scalar) = 0.5
Get indices of all values > 0.1 * final value:
gt_1_10 = find(Y > 0.1*Yf);
© National Instruments Corporation
4-19
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.21
Reformat the step response in reverse time:
Yrev = Y(151:-1:1);
Response values within 0.05 * final value:
gt_05 = find(Yrev <= 0.95*Yf | Yrev >= 1.05*Yf);
Time value of first point within bound:
ts = domain(Yrev(gt_05(1,1)))
ts (a scalar) = 5.1
Xmath Control Design Module
4-20
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. They are most often applied
to single-input, single-output (SISO) systems. With the exception of
rlocus( )and bode( ), these functions also can be used with
multi-input, multi-output (MIMO) systems.
The key principle of feedback is that the output of a system be fed back,
compared to a reference or “desired” output value, and then the error
between the two terms used to correct the system’s output so that it matches
the reference. The basic diagram of a feedback control system is shown in
Figure 5-1.
G(s)
R(s)
U(s)
Y(s)
+
–
H(s)
K
Gcl(s)
Figure 5-1. Feedback Control System Block Diagram
The output of the open-loop system is KH(s); the output of the closed-loop
system shown in Figure 5-1 is given by:
Y(s)
KH(s)
---------- = -------------------------
R(s) 1 + KH(s)
© National Instruments Corporation
5-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 5
Classical Feedback Analysis
Because open-loop systems are generally easier to study and model than
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. If you represent the open-loop transfer function H(s) as the
quotient of the numerator and denominator as follows:
H(s) = N(s) ⁄ D(s)
you can rewrite the characteristic equation of the closed-loop system as
follows:
1 + KH(s) = D(s) + KN(s) = 0
This restates the fact that the open-loop system poles (which correspond
to K = 0) are the roots of the transfer function denominator, D(s). As K
becomes large, the roots of the previous characteristic equation approach
the roots of N(s)—the zeros of the open-loop system—or infinity. For a
closed-loop system with a nonzero, finite gain K, the solutions to the
preceding equation are given by the values of s where both of the following
are true:
KH(s) = 1
∠ H(s) = (2k + 1)π
(k = 0, 1, …)
of s that correspond to pole locations for all gains K, starting at K = 0
(the open-loop poles) and ending at K = ∞.
Root locus plots provide an important indication of what gain ranges can
be used while keeping the closed-loop system stable. As discussed in the
System Stability: Poles and Zeros section of Chapter 4, System Analysis,
continuous-time systems are stable as long as their poles are in the left half
of the s-plane (have a negative real part) and discrete-time systems are
stable as long as their poles remain within the z-plane unit circle.
The Xmath root locus-plotting utility exists for SISO systems only, though
either state-space or transfer function models can be specified.
Xmath Control Design Module
5-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. In the Graphics window, the
complete locus is drawn as a solid line, with os marking the location of
zeros and xs delineating open-loop pole locations. The second syntax
brings up a window through which you can interactively modify the
closed-loop gain and see the corresponding pole locations change on
the locus.
y values can be used to restrict the range of the selected s- or z-plane. These
can be changed interactively if the interactive syntax is used. Click
RECOMPUTE to activate rate changes.
Example 5-1 shows how to plot the rool locus created in Example 2-9,
A Comparison of Several Discretization Methods.
Example 5-1
Plotting a Root Locus
H = system(0.5*polynomial([-0.36]),
makepoly([1,2.79,2.74,1.11,0.16]));
real-imaginary plane as follows:
rlocus(H, {xmin=-2, xmax=0, ymax=0.5, ymin=-0.5})
These functions give the results shown in Figure 5-3. The large xs on the
plot correspond to the open-loop pole locations you found for this system
in Example 4-1, Using poles( ) with a System in Transfer Function Form,
and the zeros correspond to the single zero at –0.36.
© National Instruments Corporation
5-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Xmath Control Design Module
5-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. A sinusoidal
input with unit amplitude and zero phase, and frequency ω produces the
following sinusoidal output:
H(jω) = A(ω)ejφ(ω)
where A is the magnitude of the response as a function of ω, and φ is the
phase. The magnitude and phase of the system output will vary depending
on the values of the system poles, zeros, and gain. In many practical
engineering applications, the system poles and zeros are not precisely
known. Because the frequency response can be determined experimentally,
undesirable parts of the system’s frequency response then can be improved
by adding known compensation to the system.
freq( )
H=freq(Sys,F,{Fmin,Fmax,npts,track,delta})
The freq( )function calculates the frequency response of a system in
several different ways, depending on the system representation. For
continuous-time transfer functions, the frequency response H(ω) at a given
frequency ω is obtained by substituting the complex frequency value jω for
qin the following equation. For discrete-time transfer functions, the value
ejwT, with T the system sampling interval, is substituted for q instead.
(q + z1)…(q + zm)
Sys(q) = ---------------------------------------------
(q + p1)…(q + pn)
For continuous-time state-space systems, the basic method for finding
frequency response is to substitute different frequency values, represented
by ω, into the following equation:
H(jw) = C(jwI – A)–1B + D
© National Instruments Corporation
5-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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(ejwTI – A)–1B + 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. Because H is zero only below the first subdiagonal,
the number of operations needed to evaluate the response expression is
proportional to the square of the size of A.
freq( )allows you to prespecify frequency ranges of interest, or it can
generate a representative frequency range from minimum and maximum
frequencies you specify. It then evaluates the complex frequency response
over those frequencies, using specialized algorithms to do this efficiently.
You can specify either a complete set of frequency points (the optional
input F) or a range of frequency points (the keyword pair Fminand Fmax)
at which to evaluate the response. The trackkeyword indicates that phase
tracking will be used to determine the values of the frequencies between
Fminand Fmax. The number of intermediate frequency points produced
using trackvaries depending on the system and the Fminand Fmaxyou
choose. Alternately, you can use the nptskeyword to specify the exact
number of logarithmically-spaced frequency points you want computed.
Specifying trackinvokes an algorithm which tracks the phase of the
frequency response to make sure that all peaks and valleys are included in
the computed response. The deltakeyword indicates the amount of phase
change (measured in degrees) to which the response evaluation should be
sensitive. If phase change between two adjacent frequency points exceeds
this delta, closer frequencies are used until either the phase change is less
than deltaor a maximum number of iterations is reached. Evaluation is
forced at key frequency points which include the poles and the points lying
halfway between adjacent poles.
freq( )returns a PDM having the frequency range as its domain. The
dependent matrices of the frequency response PDM have as many rows as
the system has outputs, and as many columns as the system has inputs. For
MIMO systems, the (i,j) element of a dependent matrix is thus interpreted
as the frequency response from input j to output i. This frequency response
forms the core of the classical control design tools discussed in this chapter.
Xmath Control Design Module
5-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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)
Y(s)
(s + 0.5)
s2(s + 2)(s + 10)
Figure 5-3. Representation of the Open-Loop System
Frequency Response of a Simple System
Example 5-2
You can create the system directly in transfer function form:
sys = system(polynomial(-0.5),
polynomial([0,0,-2,-10]));
and then obtain the frequency response directly:
H = freq(sys,{Fmin = 0.01,Fmax = 10,npts = 150});
freq( )also can be called with a predefined vector of frequency points,
or you can specify that phase tracking be used to compute frequency points
between the minimum and maximum frequencies. The number of
frequency points used with tracking will vary. To illustrate:
H = freq(sys,{Fmin=0.01,Fmax=10,track,delta=.5});
size(H)
ans (a row vector) = 1
1
335
The dynamics of this system are adequately reflected in both frequency
responses. However, systems having more closely-placed pole and zero
locations are good candidates to use with the trackkeyword.
Bode Frequency Analysis
While freq( )provides you directly with the frequency response, other
tools in the Control Design module can give you more insight into what the
open- and closed-loop frequency responses of a system imply about the
system behavior. Bode plots of system frequency response are useful
© National Instruments Corporation
5-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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].
The following complex frequency response:
H(ω) = A(ω)ejφ(ω)
can be separated into two parts, which are both functions of the frequency:
•
•
ω: the magnitude, A(ω)
the phase, φ
The magnitude can be obtained as the absolute value of the response,
whereas the phase is obtained from the four-quadrant arctangent of the
response.
The standard Bode format comprises two subplots:
•
The upper plot shows the decibel gain (the common logarithm of
the magnitude, multiplied by 20) plotted against the logarithm of the
frequency. Logarithmic (decibel) plots are a particularly useful tool
for indicating magnitude response because the multiplication of
magnitudes is shown as the sum of their logarithms, thus allowing
you to determine the system response with varying gains quickly.
•
The lower plot shows the phase, in degrees, as a function of the
logarithm of the frequency. For both the gain and phase plots,
logarithmic frequency scaling is used because it allows a wide range
of frequency-dependent behavior to be displayed simultaneously.
Because the gain and phase plots are additive for systems cascaded in
series, Bode plots of an open-loop plant and potential compensators can be
added to determine the frequency-response characteristics of the complete
which the output magnitude is reduced by three decibels, or attenuated to
approximately 70.7% (a factor of ( 2 ⁄ 2) of its original value.
Bode plots also provide an important aid to evaluate how stable—or,
more specifically, how close to instability—a closed-loop system is. As
discussed in the System Stability: Poles and Zeros section of Chapter 4,
System Analysis, for the continuous case, the closed-loop poles of a stable
system lie in the left half of the complex plane.
Xmath Control Design Module
5-8
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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
of Gcl(s). For any neutrally stable system, the frequency response
magnitude will be equal to 1 (or 0 dB) and the phase will be –180° at the
frequency at which the closed-loop roots fall on the imaginary axis.
This analysis is often applied to systems where G(s) consists of a gain, K,
and a dynamic model, H(s), in series (as shown in Figure 5-1). For cases
in which increasing the gain leads to system instability, the system will be
stable for a given value of K if the magnitude of KH(s) is less than 1 at any
how close a system is to instability by examining the value of the magnitude
and phase at these critical values. These measures are termed the gain
margin and the phase margin. These are important because real-life models
are prone to uncertainties and changes in gain or phase. Typically, systems
become unstable with gains that are too high or have too much phase lag.
Refer to Example 5-4.
The gain margin indicates by how much the gain can be raised before the
closed-loop system becomes unstable. This critical gain value at which
instability results can be thought of in several ways. As described
previously, this gain value results in the closed-loop poles of the system
being located on the imaginary axis. In terms of the Nyquist stability
criterion, discussed in more detail under nyquist( ), this is the gain value
at which the Nyquist plot crosses the negative real axis, where the phase is
–180 degrees. The gain margin itself is the reciprocal of this value,
expressed in decibels.
The Bode plot provides a clear visual interpretation of the gain margin as
the number of decibels by which the gain exceeds zero when the phase
equals –180 degrees. The phase margin is the difference between the phase
at the point where the response crosses the unit circle (has unit magnitude,
or a gain of 0 decibels) and –180 degrees. These margins provide a measure
of how near the closed-loop system roots are to instability. Depending on
the complexity of the system, there may be multiple gain and/or phase
margins.
© National Instruments Corporation
5-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
bode( )
[H,dB,Phase] = bode(Sys,{F,keywords})
The bode( )function uses freq( )to compute the frequency response
of a system. By default, the freq( )keyword trackis on, but it can be
overridden. Refer to the freq( ) section for more details. When the
frequency response His found the decibel magnitude and the phase angle
in degrees are computed as follows:
dB=20*log10(abs(H); phase=(180/pi)*atan(H)
bode( )then produces the standard Bode format plots showing response
magnitude and phase as functions of frequency. Unlike freq( ), bode( )
does not require a frequency range or a pair of maximum and minimum
frequencies; if no range is specified, it uses deffreqrange( )to
calculate a default frequency range.
bode( )often generates more than one set of plots. For MIMO systems, a
plot is made for each output with multiple curves, one per input. If there are
multiple outputs, a menu will appear which allows you to select an input to
view.
If you want to see the response of the system from Example 5-2 to input
frequencies ranging from 0.01 Hz to 10 Hz, you can analyze a frequency
response using bode( ), as shown in Example 5-3.
Example 5-3
Analyzing a Frequency Response Using bode( )
sys = polynomial(-0.5)/polynomial([0,0,-2,-10]);
[H,dB,phase]=bode(sys,
{Fmin = 0.01,Fmax=10,npts = 300,!wrap})
You obtain the gain and phase plots as shown in Figure 5-4.
Xmath Control Design Module
5-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. These gain magnitude factors add, so the slope of the
gain plot changes from –40 dB/decade to about –20 dB/decade until you
begin to see the influence of the poles at 2 radians/sec. (0.318 Hz) and
10 radians/sec (1.59 Hz), each of which contribute another –20 dB/decade
to the slope of the magnitude plot.
The phase is a function only of the pole and zero locations. Notice that in
creating the phase plot with bode( ), you specified the !wrapkeyword.
This created a phase plot where range goes down to the full angle value of
the phase, rather than wrapping the phase between +180°. Each pole at zero
contributes –90° of phase.
The remaining poles are called first-order poles because they are of the
following form:
s + pn
© National Instruments Corporation
5-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 5
Classical Feedback Analysis
Each of these contributes a phase angle φ defined by:
φ = atan(ω ⁄ pn)
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. Thus the amount of phase contributed by a first order
pole at the frequency
ω = pn
(generally termed the corner frequency, because the asymptotes used to
draw different portions of the response intersect and form a corner) is –45°.
At frequencies beyond the corner frequency, the phase angle contributed by
that pole comes closer and closer to –90°. First-order zeros contribute phase
angle in the same manner except that the sign of the angle is positive.
margin( )
[gnMargin,phMargin,dPdF,dGdF] = margin(H)
The margin( )function is a useful tool for evaluating the stability margin
of a given system based on its frequency response. It returns PDMs
indicating the gain margin and the phase margin, as well as the rate of
change of gain and phase.
margin( )is defined for SISO systems only. It takes as input either a
single PDM representing frequency response or a pair of PDMs containing
gain information in decibels and phase information in degrees. In either
case, the domain of the input is the set of frequency points, ω.
Within margin( ), as within bode( ), the frequency response is
converted to decibel magnitude and degree phase. All angles are converted
to four-quadrant angles between 0° and 360°. Use the following notation
for each point i in the frequency range:
Δx = x(i + 1) – x(i)
Xmath Control Design Module
5-12
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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
Δω
⎛
⎝
⎞
⎠
------------------
---------------
Δω – gain(i + 1)
Mp(i) = phase(i) +
Δω
Δgain
(gain(i) + Δgain)
Δω
⎛
⎞
⎠
---------------------------------------------
------------------
Δω – phase(i + 1)
Mg(i) = –
⎝
Δphase
Δphase
This loop finds all the frequency intervals within the response which
contain –180° phase crossings and 0 decibel gain crossings. margin( )
then interpolates to find more exact frequency values for the crossings.
A gain margin value is returned for every pair of phases between which
a –180° phase value must occur, and a phase margin is returned for each
pair of gains between which a zero-decibel gain value must occur.
margin( )also computes the frequency-rate of change for both the phase
and the gain of the response.
the frequency response are too widely spaced.
You can examine the gain and phase margins of your open-loop system
quickly using margin( ), without having to draw the bode gain and phase
response plots first. Input to margin is the frequency response Hof the
system. Referring to the system defined in Example 5-3, you can see that
you already have Has the output from bode( ), but you can calculate it
explicitly using freq( )as shown in Example 5-4.
Example 5-4
Obtaining Gain and Phase Margin Using margin( )
H = freq(sys,{Fmin=0.01,Fmax=10,npts=300});
[Gm,Pm] = margin(H)
Gm (a pdm) =
- domain |
----------+----------
0.595514 | 44.5062
Pm (a pdm) =
domain
|
----------+----------
0.0257558 | 12.3814
© National Instruments Corporation
5-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The open-loop
frequency response is calculated and plotted against the gain in the standard
Nichols format (gain in decibels versus phase in degrees). Different points
on the plot thus correspond to different values of ω.
nichols( )plots are particularly useful as a means of obtaining the
closed-loop frequency response of a system from the open-loop response.
Nichols plots are frequently augmented with curves, or loci, of constant
magnitude or phase. These curves are drawn when the patternkeyword
is specified. Notice that each point on the open-loop response curve
corresponds to the response of the system at a given frequency, and the
closed-loop magnitude response at that frequency can be read off the
Nichols plot by noting the value of the magnitude locus which the point on
the curve intersects. The closed-loop phase can be determined in a similar
manner by noting the phase locus which the open-loop curve crosses.
[Oga70]
nichols( )is implemented in a manner very similar to that used for
calls freq( )internally, and converts the complex frequency response
to magnitude gain in decibels and phase in degrees. bode( )and
nichols( )differ only in the plots they produce. For MIMO systems,
nichols( )will produce plots with as many curves as there are system
inputs. A menu presents a selection of output responses. To generate a
Nichols plot, use the syntax shown in Example 5-5.
Example 5-5
nichols( ) Plot
A = [2, 0, -0.01; 2,-2,0; -1.4, 3, 0];
B = [3; 5; -1];
C = [1, 0, 4];
nsys = system(A,B,C,1);
H=nichols(nsys,
{Fmin=.01,Fmax=5,npts=300, pattern,!wrap})
Xmath Control Design Module
5-14
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. Nyquist plots are particularly useful in that they indicate the
stability of a closed-loop system, given an open-loop system which
includes a gain, K (it may be unity).
Nyquist’s stability criterion derives from Cauchy’s principle, which states
that a contour integral of a complex function will evaluate to zero as long
as the contour does not contain a singularity of that function [ChB84]. The
frequency response is the complex function in this case, and the contour
over which it is evaluated and plotted is determined by the frequency range
of the response.
Nyquist’s stability criterion states that the number of clockwise
encirclements of the –1 point on the real axis by the plot is equal to the
number of unstable closed-loop poles minus the number of unstable
open-loop poles. This criterion can be used to determine how many
encirclements are required for closed-loop stability. For example, if the
© National Instruments Corporation
5-15
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
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. You can
derive any plot from a context menu. If you close a feedback loop around a
SISO system in transfer function format, you obtain a closed-loop system
as shown in Figure 5-6.
U(s)
Y(s)
+
–
num(s)
den(s)
H(s) =
K
Figure 5-6. Closed-Loop System Containing a Variable Gain K
You obtain the following closed-loop transfer function from Y(s) to U(s):
Y(s)
KH(s)
----------- = -------------------------
U(s) 1 + KH(s)
Thus, the closed-loop roots are the roots of the equation 1 + KH(s) = 0.
The complex frequency response of KH(s), evaluated for s = jω in
continuous time and ejωT for discrete systems, will encircle (–1,0) in the
complex plane if 1 + KH(s) encircles (0,0). If you are examining the
Nyquist plot of H(s), you will notice that an encirclement of (–1/K,0) by
H(s) is the same as an encirclement of (–1,0) by KH(s). This fact allows you
to use one Nyquist plot to determine the stability of a system for any and
all values of K.
nyquist( )
H = nyquist(Sys,{F,keywords})
The nyquist( )function is structured very similarly to bode( )and
nichols( )in that it is largely a wrapper on the freq( )function to
obtain the system’s frequency response. The output His just the output from
the call to freq( ). The main difference from the other two functions is
system’s response. It generates the Nyquist plot by plotting the real part of
each point of the response against the imaginary part.
Xmath Control Design Module
5-16
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 5
Classical Feedback Analysis
The Nyquist plot Xmath generates is complete only for the frequencies you
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 +∞. Extend
the resulting curve, traveling clockwise, until the contour is closed. Refer
the augmented plots in Example 5-6. When you have done this, you can use
roots, Z, given the number of clockwise encirclements of the (–1/K,0) or
(–1,0) point and the number of unstable (right-half plane) poles of the
open-loop system.
For an example of how to use Nyquist plots to determine stable gains for
the closed-loop system, refer to Example 5-6.
Example 5-6
Using Nyquist Plots to Determine Stable Gains for the Closed-Loop System
By examining the Nyquist plot for your open-loop system
(s + 0.5)
G(s) = ----------------------------------------
s2(s + 2)(s + 10)
you can tell for what multiplicative gain values K the closed-loop system
will be unstable.
H = nyquist(sys,{Fmin=0.01,Fmax=10,npts=300});
but the plot gives more information than you need about the low frequency
response and not enough about the response at higher frequencies. Refer to
Figure 5-7.
You do another Nyquist plot, this time examining the high-frequency
response more closely. Refer to Figure 5-8.
H2= nyquist(sys,{Fmin=.5,Fmax=5,npts = 150})
© National Instruments Corporation
5-17
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. In this figure,
Nyquist contour is formed by drawing the system’s Nyquist plot for all
positive frequencies, reflecting it about the real axis to show plot for
negative frequencies, and completing the closed contour in a clockwise
direction.
Imaginary
Real
Figure 5-9. Nyquist Contour Formed by Drawing the System’s Nyquist Plot
for All Positive Frequencies
Because your open-loop system has no unstable (right-half plane) poles,
the number of unstable closed-loop poles for a given gain K will be equal
to the number of times the contour encircles the (–1/K, 0) point.
Referring back to the Nyquist plots, you see that for K > 168
(or –1/K > –0.006), you have two encirclements of the (–1/K, 0) point,
and thus two unstable closed-loop poles. For gain values less than 167,
the closed-loop system is stable. You can verify this with a small
experiment, using a value of 169 for K:
sysu = 169*sys;
sysucl = feedback(sysu);
poles(sysucl)
ans (a column vector) =
© National Instruments Corporation
5-19
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The power spectral density of a time series is defined
as the Fourier transform of the autocorrelation function of the series.
Given these two concepts, you can obtain the power spectral density of the
output of a linear, time-invariant system just by knowing the power spectral
density of the input and the system’s transfer function [Leo89], [GrD86].
Representing the transfer function by H(q) and the power spectral densities
of the input and output as SU(q) and SY(q), respectively:
SY(q) = H(q) 2SU(q)
You also can obtain the cross-power spectral densities:
SYU(q) = H(q)SU(q)
∗
SUY(q) = H (q)SU(q)
These results indicate that you can shape the spectrum of a linear system’s
output by using an input with an appropriate spectrum. Alternatively, you
can choose a system to give you the output spectrum you want, given a
fixed set of input data. When you use linear systems in transfer-function
form for such applications, you generally refer to them as filters rather than
systems.
psd( )
[Ypsd,Yspec] = psd(Sys,{Uspec})
cross-spectral density of a system’s outputs as a function of frequency,
given the frequency-dependent input power spectral density matrices. The
input parameter Uspecis a PDM where domain contains the frequency
Xmath Control Design Module
5-20
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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:
Yspec = H × (jw) × Uspecjw × H(–jw)
psd( )calls freq( )internally to compute the frequency response, H, of
the system. It uses the frequency range specified by the domain of Uspec.
The power spectral density of the output as a function of frequency, given in
Ypsd, is obtained from the real parts of the diagonal terms of the dependent
matrices in Yspec.
Some background information on power spectral density may be useful.
Given a time-domain input series U(t), the power spectral density of U(t)
is the Fourier transform of the autocorrelation of U(t). For a system with
q inputs each input spectral density dependent matrix within Uspecis a
square Hermitian matrix of size q. A Hermitian matrix is a square matrix
equal to its complex conjugate transpose. If Uspecis constant for all
frequencies (when the spectrum is white) then Uspeccan be specified as
a single matrix.
series, the output power spectral density of the first system can be used
as the input power density to the second system in a subsequent use of
psd( ).
For an example of how to verify the response of a system to white noise
input, refer to Example 5-7.
Example 5-7
Verifying the Response of a System to White Noise Input
You can easily generate the power spectral density of an input white noise
process.
sys = polynomial(-0.5)/polynomial([0,0,-2,-10]);
w = logspace(0.01,1,50);
Uspec = pdm(ones(w),w);
You then use psd( )to obtain the output power spectral density and the
cross-spectral density as a function of frequency.
[Ypsd,Yspec] = psd(sys,Uspec);
© National Instruments Corporation
5-21
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. First, you need to assess the controllability and observability of the
system. The designs discussed in this chapter are based on systems that are
both controllable and observable. When you have determined the
controllability and observability of the system, you can design a feedback
control law based on the set of state values. Next, you design an estimator
that estimates the state variable values based on the measured output.
Finally, you combine the controller and estimator to obtain a complete
compensator for the system.
want to optimize for a given system. This performance index is a quadratic
function reflecting the physical constraints of the system and the
characteristics of any noise that may be present. When this performance
index is a quadratic, you solve mathematically for the optimal control law
and estimator as discussed in the Linear Quadratic Regulator section and
the Linear Quadratic Estimator section.
This chapter concludes with a discussion of system balancing. The
controllability and observability grammians provide a measure of how
controllable and observable a system is. They also can be used to transform
a system to its internally balanced form.
Controllability
Controllability is the property of being able to move the states of a system
arbitrarily in a finite time, given some control input to the system. Although
a particular physical system may be controllable by this definition, not all
state-space models describing that system may be controllable. For
example, if there exists a system eigenvector orthogonal to the input
© National Instruments Corporation
6-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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).
It can be shown (refer to [Kai80]) that for a continuous-time system with
the state update equation:
·
x = Ax + Bu
(6-1)
you can define the controllability matrix for both continuous and discrete
systems as:
C = [B AB A2B…An – 1B]
(6-2)
For all modes of the system to be controllable, the controllability matrix C
must contain a linearly independent column vectors for each system mode.
Thus, with A an n × n matrix, C must have rank n for the system to be
controllable.
In the context of gain-state feedback, a system’s controllability determines
whether you may be able to change the effective dynamics of the system to
ones that yield a more desirable response.
Using full-state feedback, as shown in Figure 6-1, so that u = v – Kx.
Working through the system equations, you obtain
·
x = (A – BK)x + Bv
(6-3)
for the new state-update equation. If the system is controllable, you can
relocate the eigenvalues of the closed-loop system to any value by choosing
the vector of state gains K appropriately. Conversely, the eigenvalues
associated with uncontrollable modes remain unchanged, no matter what
value you choose for K.
Xmath Control Design Module
6-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
v
u
y
+
–
x = Ax + Bu
y = Cx + Du
K
Figure 6-1. Full-State Feedback Being Used to Relocate the Eigenvalues
of a Controllable System Based on the Value of the Gain K
controllable( )
[SysC,T,nuc]=controllable(Sys,{tol})
controllable in a given system?” The controllable( )function returns
the controllable partition of a state-space system, the number of
uncontrollable states in the original system, and a linear transformation
matrix which can be used to partition the states into controllable and
uncontrollable sets. For an example of how this is done, refer to
Example 6-1.
controllable( )uses the staircase algorithm, which is discussed in
more detail in the stair( ) section.
Example 6-1
Controllability of a System
Perform controllable( )on 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);
[SysC,T,nuc] = controllable(Sys)
The system has 2 uncontrollable states
SysC (a state space system) =
A
1
B
-1
C
-0.6
D
0
© National Instruments Corporation
6-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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
0
1
0
-1
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.
Similarly, if you form the controllability matrix for this system,
[,states] = size(A);
Con = B;
For i = 1:states-1;
Con = [B, A*Con];
endFor
det(Con)
ans (a scalar) = 0
confirming the results from controllable( ).
Observability and Estimation
As described in the Controllability section, the term controllability
describes whether or not a system’s states can be affected, and the system
eigenvalues relocated, by changes to the system input. The analogous
concept of observability describes whether it is possible to determine the
value of an individual state at a particular time by observing the system
outputs for a finite amount of time. In essence, an observable system is one
for which you can “observe” state values by knowing the output of the
system.
Xmath Control Design Module
6-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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
·
·
=
x +
CA
CA2
y
CB 0 0 u
··
··
y
CAB CB 0 u
Generally, you term the following matrix
C
CA
CA2
O =
…
CAn – 1
the observability matrix. If the rank of O is equal to n (where A is an n × n
matrix), you can always find a state vector which will realize the initial
conditions you want on the output, presuming that you know the initial
conditions on the input. If, however, the observability matrix loses rank
(is singular, in the SISO case), you will not be able to find states that give
you particular output conditions if those conditions lie in the null space of
the observability matrix.
The observability of a system is of particular importance when you want
system dynamics and the input and output values at a given time. If a system
values of the states and use the information available to zero the state
estimate error as quickly and accurately as possible.
Referring to Figure 6-2 and tracing through the system equations, you can
ˆ
˜
obtain the time-update equation for the state estimate error x = x – x.
Figure 6-2 is a general observer block diagram where the output estimate
ˆ
˜
˜
error y is defined as y = y – y.
© National Instruments Corporation
6-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
x
x
u
y
y
x = Ax + Bu
x = A + Bu + Ly
L
Cx
+
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.
·
˜
(6-4)
The problem of finding the eigenvalues of (A – LC) can be equivalently
posed as that of finding the eigenvalues of (A' – C'L'). This statement can
state-feedback controller (refer to the new state-update equation in the
Controllability section), with A, B, and K replaced by A', C', and L',
respectively. Notice that these two representations correspond to a
state-space system and its transpose. This illustrates the principle of duality
between the controller and estimator forms. For more information, refer to
the Duality and Pole Placement section.
observable( )
[SysO,T,nuo] = observable(Sys,{tol})
The observable( )function is the analogue to controllable( ).
controllable, its transpose {A',C',B',D'} is observable. observable( )
returns the observable partition of a state-space system, the number of
matrix which can be used to partition the states into observable and
unobservable sets. For an example of how to use the observable( )
function, refer to Example 6-2.
observable( )uses the staircase algorithm, which is described in more
detail in the stair( ) section.
Xmath Control Design Module
6-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
Example 6-2
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.
Similarly, if you form the observability matrix for this system,
[,states] = size(A);
Obs = C;
For i = 1:states-1;
Obs = [C; Obs*A];
endFor
det(Obs)
ans (a scalar) = 0
you see that the observability matrix is singular (its determinant is zero),
confirming the results you saw from observable( ).
Minimal Realizations
All state-space systems have an infinite number of realizations. All systems
have a minimum number of states needed to express the system dynamics,
but can be described using any number of states greater than or equal to this
minimum number. If a system has more states than are needed to express a
given transfer function, it will have unobservable and/or uncontrollable
modes corresponding to eigenvalues of the A matrix that are not poles of
the transfer function.
All minimal realizations of the same system are related by a coordinate
transformation.
© National Instruments Corporation
6-7
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
minimal( )is implemented directly as a wrapper on controllable( )
and observable( ). The controllable subsystem is extracted first, then
the observable part of the subsystem is returned. For an example of how to
find a minimal realization for a system with uncontrollable or unobservable
parts, refer to Example 6-3.
Example 6-3
Finding a Minimal Realization for 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);
Notice that the system has a number of zero-pole pairs which cancel each
other out:
poles(Sys)
ans (a column vector) =
1
1
1
zeros(Sys)
ans (a column vector) =
1
1
To find the minimal part of the system:
SysM = minimal(Sys);
The system has 2 uncontrollable states
poles(SysM)
Xmath Control Design Module
6-8
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The matrices are
partitioned as shown in the following example:
Auc
0
0
Astair
=
Bstair
=
Bc
Acuc Ac
where Ac is controllable, Auc is uncontrollable, and Acuc represents the
coupling from the uncontrollable states to the controllable part of the
representation. There is no coupling from the controllable states to the
uncontrollable ones.
T is a transformation matrix between the original system, Sys, and the
system in staircase form, SysT.
The transformations are as follows:
Astair = T –1AT
Bstair = T –1B
Cstair = CT
The optional tolerance, tol, indicates the threshold value beneath which
numbers in the transformed matrices should be rounded to zero.
The staircase algorithm used to partition the system derives from the
and [BeV88].
© National Instruments Corporation
6-9
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
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). Given the system’s state-space representation and any
desired set of closed-loop poles, you can solve an eigenvalue problem to
poles of a minimal system can be moved to any value, this approach does
not guarantee that the resulting gain is small or physically practical—just
that it is finite.
The similarity between the new state-update equation (Equation 6-1)
and the time-update equation (Equation 6-4) for controller and observer
feedback brings up the principle of duality with respect to the
controllability and observability of a system. Briefly, for a given state-space
system Sys1 with system matrices {A,B,C,D}, there exists a dual system
Sys2 described by {A*,C*,B*,D}, using * to denote a complex conjugate
transpose. If Sys1 is controllable, Sys2 will be observable, and vice versa.
This can be quickly verified by constructing the controllability and
observability matrices for both. Thus, the gain value that yields a set of
desired closed-loop poles for feedback control of a system also yields an
observer with the same pole locations for the system’s dual.
K = poleplace(A,B,poles)
The poleplace( )function solves the problem
eig(A – B * K) = polesfor single-input systems. This is essentially
the problem posed in Figure 6-1 and the new state-update equation in the
Controllability section. If you know where you want the system poles to
be located, poleplace( )returns the value of the gain vector Kthat will
move the closed-loop poles to the desired locations.
The syntax poleplace(A',C',poles)can be used for regulator
problems, or by duality, for estimator problems. In general, the system
{A,B} must be reachable (all unstable poles controllable) for controller
design and the system {A',C'} must be stabilizable (all unstable poles
observable) for estimator design. The current poleplace( )
implementation is limited to single-input systems.
Xmath Control Design Module
6-10
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
poleplace( )is unusual among Xmath’s modern control design
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. For an example of an arbitrary pole placement for a controllable
system, refer to Example 6-4.
Example 6-4
Arbitrary Pole Placement for a Controllable System
A = [0,1,0,0;21,0,0,0.8;0,0,0,1;0,0,0,-4];
B = [0,-2,0,1]';
C = eye(4,4);
D = zeros(4,1);
ipsys = system(A,B,C,D);
If you want to place the system poles in a Butterworth pattern:
Kc = poleplace(A,B,[-5+8.66*jay, -8.66+5*jay]);
You then can use this new gain vector as a feedback gain to create a new
system,
ipsysfb = feedback(ipsys, system([],[],[],Kc));
and verify that the poles of this new system are at the designated locations:
poles(ipsysfb)
ans (a column vector) =
-8.66 + 5 j
-8.66 - 5 j
-5 + 8.66 j
-5 - 8.66 j
You need specify only one complex pole in a conjugate pair of desired pole
locations. poleplace( )checks for conjugate pairs and adds conjugates
as necessary to the input poles( )vector.
Then the system matrix S is created:
A B
S =
r 0
r is a random row vector with as many rows as A has columns. You then
create a random complex row vector with as many elements and conjugate
© National Instruments Corporation
6-11
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. poleplace( )then divides the last element of this quotient vector
by the negative of the first element of the quotient vector, and the result is
the gain required to move that pole value. This sequence of steps is
performed as a matrix operation so that the complete gain vector is
computed immediately. rcond( )is called to examine the condition of the
matrix formed by all the quotient vectors. If the condition number returned
is less than eps× (the row size of A), poleplace( )displays a warning
message indicating that the eigenvectors of the closed-loop system are
ill-conditioned.
Linear Quadratic Regulator
A regulator is a feedback controller designed to drive the states of a
controllable system using acceptable amounts of control and keeping the
states within acceptable levels (where the designer can mathematically
define what constitutes “acceptable” in both cases). Figure 6-3 shows a
continuous-time regulator where the design presumes availability of all
states, feeding them back through the optimal gain array Kr to drive the
system so that the states return to zero as quickly as possible in the presence
of a disturbance or noise, represented by ω.
u
x
x = Ax + Bu
–Kr
Figure 6-3. Continuous-Time Regulator
In designing a regulator, the goal is to find a controller that minimizes the
effects of disturbances on the states of the system. Xmath’s linear quadratic
regulator function, regulator( ), uses a quadratic performance index to
establish the trade-off between the permissible state fluctuation and the
available energy, or amount of control, required to move the states.
Note In designing a regulator, assume that all the states of the system are available as
outputs.
Xmath Control Design Module
6-12
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
For continuous-time systems,
·
x = Ax + Bu
the quadratic performance index takes the form:
∞
Rxx Rxu
x(t)
J =
dt
x'(t) u'(t)
∫
Rxu' Ruu u(t)
0
For the discrete case where the system is defined as a multistage process:
xk + 1 = Axk + Buk
the performance index is defined similarly except that a summation sign
replaces the integral of the preceding quadratic performance index
equation.
Rxx is a real, symmetric, positive-semidefinite matrix indicating the
weighting of the cost on the elements of the state vector x. Ruu is a real,
symmetric, positive-definite matrix indicating the weighting of the cost on
the control inputs given by the vector u. Rxu is a real matrix indicating the
cross-weighting of the cost between states and inputs; for many
applications, it will consist of all zeros if the control and states are
uncorrelated.
Bryson and Ho showed in [BH75] that the optimal control which
minimized this quadratic performance index is a linear feedback
combination of the states, u = Krx, for both the continuous and discrete
cases.
For the continuous case, Kr is defined as follows, with P solving the
continuous-time Riccati equation:
–1
Kr = R (B'P + Rxu)
uu
–1
uu
Rxx + PA + A′P – (PB + Rxu)R (Rxu′ + B′P) = 0
and for the discrete case, P solves the discrete Riccati equation.
Kr = (Ruu + B'PB)–1(B'PA + Rxu)
A'PA – (A'PB + Rxu)(Ruu + B'PB)–1(B'PA + R′xu) + Rxx = P
© National Instruments Corporation
6-13
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
regulator( )
[Kr,ev,P] = regulator(Sys,Rxx,Ruu,{Rxu})
The regulator( )function calculates the optimal gain matrix Kr for a
given dynamic system with specified state weighting, control weighting,
and (optionally) cross-weighting matrices.
Alternatively, Kr can be obtained through a call to riccati( ):
[P,resid,Kr,ev]=riccati(Sys,Rxx,Ruu,{S=Rxu})
The syntax for riccati( )is discussed in the Riccati Equation section.
As shown in the diagram of a continuous-time regulator in Figure 6-3, the
state equation for the regulator is the following:
·
x = (A – BK)rx
If you want the closed-loop system eigenvalues, compute them as the
eigenvalues of (A – BKr).
If numerical difficulties are encountered, the algorithm will attempt
to determine whether or not the problem is well posed. Checks are made
to determine stabilizability and the positive definiteness or
semipositive-definiteness of the cost functionals.
The most important design parameters are Rxx and Ruu, which need to be
chosen to reflect the real limitations on how much control can be provided,
or how problematic large state values can be. For an example of how to
design a regulator for the inverted pendulum, refer to Example 6-5.
Example 6-5
Designing a Regulator for the Inverted Pendulum
A classic control design problem, the inverted pendulum, consists of a rod
(the pendulum) hinged to the top of a cart which can be moved freely in
either direction along a line. The goal of the controller is to supply an input
in Figure 6-4).
Xmath Control Design Module
6-14
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
f
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.
·
·
Designating [φ φ x x]′ as the state vector, you can set up the system in
Xmath:
A = [0,1,0,0;21,0,0,0.8;0,0,0,1;0,0,0,-4];
B = [0,-2,0,1]';
C = [1,0,0,0;0,0,1,0];
D = [0,0]';
ipsys = system(A,B,C,D);
You design a regulator with the assumption that all four states are available.
·
·
Recalling that you defined the state vector as [φ φ x x]′ , you can decide
the weighting you want to associate with each state. Refer to the quadratic
performance index equation in the Linear Quadratic Regulator section for
more information. For this particular problem, your most important
performance goal is that the pendulum stay upright—that is, that φ be
tightly controlled to stay as close to zero as possible. You also might prefer,
though to a lesser extent, that you not have to move the cart over too great
a distance. Physical limitations, such as the size of the room in which the
experiment is conducted, should be considered.
If you are not particularly concerned about the speed of the cart across the
floor or that rate of change of the angle, you might define,
Rxx = diagonal([1,0,0.1,0]);
with the larger values in the matrix corresponding to states whose values
you care most about. Presuming you are not too worried about the size of
the input u you impart to the cart:
Ruu = 1e-5;
© National Instruments Corporation
6-15
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
Ruuis 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 -32.1056 -100 -27.3036
Note 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. In the real world,
only the output values are generally available and they are frequently
section that you can obtain an estimate of the states using an observer if the
system is reachable. The problem solved with the optimal estimator
function estimator( )is that of finding the best estimate of the states,
given certain assumptions about the noise associated with the output.
As shown for the continuous case in Figure 6-5, the plant system is
augmented with an estimator—an observer used in conjunction with a
noisy system. The estimator supplies estimates of all the system states and
feeds back the difference between the estimated and the actual outputs
through the optimal estimator gain Ke. estimator( )calculates the
constant, optimal state-estimator gain matrix Ke for a dynamic system. The
estimator gain is derived by minimizing the expected mean square of the
ˆ
error between the measured output y and the output from the estimator, y.
This model takes into account that there may be some process noise within
the system model (plant) itself as well as some noise inherent in the device
used to measure the outputs.
Xmath Control Design Module
6-16
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
+
y
x
x
u
x = Ax + Bu + Gw
x = Ax + Bu + Ke(y – y)
Ke
Cx + Du
Cx + Du
+
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. For a linear–time–invariant process
described by:
Sys = system(A,B,C,D)
The following equation describes the complete plant:
·
x = Ax + Bu + Gw
y = Cx + Du + n
A, B, C, and D are directly from the previous state-space system
representation where ω is the input disturbance, G is the input disturbance
matrix and ν is the measurement noise. The noise intensity matrices are
defined as,
E(v(t)v′(τ)) = Qyyδ(t – τ)
E(Gω(t)ω′(τ)G') = Qxxδ(t – τ)
E(Gω(t)v′(τ)) = Qxyδ(t – τ)
where E is the expectation operator and δ is the delta function.
The noises ω and ν are assumed to be white and zero mean. Qyy has matrix
dimensions equal to the number of plant outputs and must be positive
definite, while Qxx has matrix dimensions equal to the number of plant
states and must be positive semi-definite. In many cases the input
disturbances and output noises are uncorrelated so that Qxy = 0. If
© National Instruments Corporation
6-17
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Because not all the values in the state vector are directly available from
measurements, your goal is to find an estimate of the state vector which
minimizes, in a least-squares sense, the error between the actual state vector
ˆ
and the estimated state vector. This estimated vector is denoted by x.
Because you want to minimize the error between this estimate and the
actual state values, the quadratic expression to be minimized becomes:
J =
•
ˆ
ˆ
(x(t) – x(t))′ (y(t) – y(t))′
∫
0
ˆ
Qxx Qxy
(x(t) – x(t))
dt
ˆ
Qxy′ Qyy (y(t) – y(t))
For the case of a discrete-time system, this quadratic expression is
evaluated as a summation rather than as an integral. No additional
information is provided by the inclusion of the Du term, so it can be
omitted without loss of generality.
A derivation of the differential equation for the continuous-time state
ˆ
vector estimate, x, can be found in [Kai81]. In the limit, this differential
equation, which provides the values for the continuous-time optimal
estimator, is
·
ˆ
ˆ
x = (A – KeC)x + Bu + Key
(6-5)
where Ke = (PC' + Qxy')Qyy–1 and where the matrix P is obtained by solving
the algebraic Riccati differential equation:
–1
yy
0 = PA' + AP – (PC' + Qxy)Q (Qxy' + CP) + Qxx
The two preceding equations describe the continuous-time Kalman-Bucy
filter [KaB61].
Xmath Control Design Module
6-18
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
xk + 1 = Axk + Buk + Gω
yk = Cxk + Duk + ν
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. This
pre-existing estimate is called xk . The estimated state value after each
ˆ
measurement update is denoted by xk . This method takes into account the
fact that the system’s states change between measurements due to the
system dynamics. The optimization problem, then, consists of minimizing
the estimate error covariance M after each measurement update. This
minimization is performed in [Kai81]. This problem is expressed in the
same manner as in the preceding quadratic expression (for J), except that a
summation sign replaces the integral as you are working with discrete data
and you replace the variable J with M, to denote that this covariance follows
each measurement update.
In determining the state values xk from each measured yyk, consider the time
just prior to a new measurement for yk. At this point xk and Mk are the
current estimates for the state and covariance. xk is derived from the
ˆ
previous measurement xk – 1, and Mk is derived from the previous
post-measurement error covariance, Pk – 1, as shown in Equation 6-6.
ˆ
xk = Axk – 1 + Gk
Mk = APk – 1A' + GQxxG'
(6-6)
This is referred to as the time update, because you are propagating the state
forward in time until the next measurement arrives.
Then, at the time immediately following the measurement, you effect the
measurement update. This reflects the new information in an improved
ˆ
state estimate xk and a somewhat smaller covariance, Pk. The equations for
© National Instruments Corporation
6-19
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
this measurement update, derived in [Kal60], are shown in the following
equations.
ˆ
Pk = (Mk–1 + C'Q–1C)–1
yy
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.
A'PA – A'PC'(Qyy + CPC')–1CPA + Qxx = P
(6-7)
where
–1
A = A' – C'Q Qxy
yy
Qxx = Qxx – Qxy'Q–yy1Qxy
and the discrete feedback gain Ke is given by
K′e = (Qyy + CPC′)–1(CPA′ + Qxy)
(6-8)
estimator( )
[Ke,ev,P] = estimator(Sys,Qxx,Qyy,{Qxy})
The estimator( )function calculates the optimal gain matrix Ke
for a given dynamic system with specified process, measurement, and
(optionally) cross-weighting noise matrices.
Alternatively, Ke can be obtained through a call to riccati( ):
[P,resid,Ke,ev]=riccati(Sys',Qxx,Qyy,{S=Qxy})
The syntax for riccati( )is described in the Riccati Equation section.
As shown in the estimator diagram in Figure 6-4, the state equation for the
estimator is:
·
ˆ
ˆ
x = (A – KeC)x + (B – D)u + Gω
Xmath Control Design Module
6-20
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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. Returning to the inverted-pendulum plant of
Example 6-5, assume that internal disturbances enter the system with the
inputs. You thus can define a Qxx, which is a function of the input matrix:
Qxx = 0.25*B*B';
Similarly, allow for some disturbance noise affecting the accuracy of the
outputs you measure from the system. You have two output measurements
for this system, thus two separate sources of noise. Assume that the noise
affecting one output measurement does not affect that other, and that the
effects of measurement noise are rather small for this instance.
Qyy = diagonal([1e-6,3e-6]);
[Ke,ev,P] = estimator(ipsys,Qxx,Qyy);
Ke
Ke (a rectangular matrix) =
42.6012 -6.21395
965.35 -159.818
-18.6419
-401.88
4.67597
68.8522
Now that you have access to a set of augmented states for the system (found
with the differential equation for the continuous state vector shown in
Equation 6-5), you can find the optimal controller based on the assumption
of full-state feedback.
Linear Quadratic Gaussian Compensation
Many real-world control system design problems lend themselves to
solutions using a regulator, except that not all the states are available as
directly measured or computed outputs.
A compensator combines your ability to control a system using full state
feedback with our ability to estimate the system states given the system
output. You can design the controller and estimator separately and then
combine them to make the system respond as desired, based on the
measured output. The combination of system, controller, and estimator into
© National Instruments Corporation
6-21
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Combining the plant, or system, equations with those of the regulator and
estimator, you can simplify the system equations for the compensator as
follows:
·
ˆ
ˆ
ˆ
x = Ax + Bu + Ke(y – (Cx + Du))
·
ˆ
ˆ
ˆ
x = Ax – KeCx – (KeD – B)u + Key
ˆ
ˆ
x = [A – KeC – (B – KeD)Kr]x + Key
ˆ
u = – Krx + (0)y
(6-9)
Controller
Estimator
x = Ax + Bu + Ke(y – (Cx + Du))
w
n
x
Plant
x = Ax + Bu
y = Cx + Du
y
Regulator
–Kr
u
+
+
u
y
Figure 6-6. Linear Quadratic Gaussian Compensator (in the Bold Rectangle)
Equation 6-9 describes the state-space equations for both the
continuous-time compensator and the discrete-time compensator if no unit
delay is used between the time at which an input arrives at the system and
the time at which the new output appears. However, if you are working with
a real-time system which enforces a unit delay between the measurement
and the control update, you will need to create a “direct” compensator
in predictor form. With this direct implementation, the system output
equations become the same as the state update equations, multiplied by
a factor of the regulator gain Kr.
Xmath Control Design Module
6-22
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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 SysCis always in state-space
form.
The regulator and estimator gains Kr and Ke need to have been calculated
estimator( )to compute these gains if they are not already available.
These functions give you the option to incorporate the presence of the white
process and measurement noises ω and ν in your model, as shown in
Figure 6-6.
Example 6-7 uses the inverted pendulum problem estimator and regulator
gain vectors (obtained in the preceding two examples) to form a
compensator. Notice that for this example you use the optional fields of the
system( )function to store information that will be useful when you
combine and simulate the compensator-plant system.
Example 6-7
Combining the Regulator and Estimator into a Full Compensator
A = [0,1,0,0; 21,0,0,0.8; 0,0,0,1; 0,0,0,–4];
B = [0,–2,0,1]';
C = [1,0,0,0;0,0,1,0];
D = [0;0];
ipsys = system(A,B,C,D, {inputNames = "Force (u)",
outputNames = ["phi","x"],
stateNames = ["phi","d(phi)","x","d(x)"]});
Kr=regulator(ipsys,diagonal([1,0,.1,0]),1e-5);
Ke=estimator(ipsys,0.25*B*B',
diagonal([1e-6,3e-6]));
ipsysc=lqgcomp(ipsys,Kr,Ke);
Now that you have the compensator, you need to connect it to the original
plant:
ipsyscl = afeedback(ipsys, ipsysc);
The system now has eight states (four from the compensator, four from the
original plant), three individual inputs (u and y, where y comprises φ and x)
and three outputs (again, u and y). Here you have not added additional
inputs to the system for process and measurement noise. Now you can
© National Instruments Corporation
6-23
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.01:15;
u = pdm([sin(t/2); zeros(t);zeros(t)],t);ycl =
ipsyscl*u;
[outNames] = names(ycl);
Set plot attributes for all three plots:
p1=plot ({xlab = "Time",ylab = "Amplitude",
columns = 1,rows = 3, hold})
for i = 1:3
p1=plot (ycl(i,1), {graph_number = i,
legend = outNames(i)});
endfor
plot(p1,{!hold})
Note The different y-axis scaling for each subplot is shown in Figure 6-7. This figure
shows φ and x as a function of time, starting from zero, as a result of a sinusoidal force
applied to the system input.
Xmath Control Design Module
6-24
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
Figure6-7. fandx asaFunctionofTime,StartingfromZero,asaResultofaSinusoidal
Force Applied to the System Input
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.
The continuous-time Riccati equation is given by:
A′P + PA – (PB + S) × inv(R)(B′P + S′) + Q = 0
The discrete-time Riccati equation is given by:
A′PA – P – (A′PB + S) × inv(R + B′PB)(B′PA + S′) + Q = 0
This function can be used to solve the optimal regulator problem, and by
duality, the optimal estimator problem. An alternative form of the
© National Instruments Corporation
6-25
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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, Acan be either a matrix or a system object.
If Ais a matrix, riccati( )solves the continuous Riccati equation in the
Riccati Equation section unless the B matrix is present; then it solves the
discrete form (also shown in the Riccati Equation section).
If Ais a system object, the solution method depends on whether the system
is continuous or discrete. The B matrix is unnecessary to distinguish
continuous from discrete.
The algorithms used are based on [Lau79] and [PLS80]. For the continuous
case, an ordinary Schur solver is used. For the discrete-time case, the
solution uses a generalized eigenvalue solver. For an example of a
continuous Riccati equation, refer to Example 6-8. For an example of
a discrete Riccati equation, refer to Example 6-9.
Example 6-8
Continuous Riccati Equation
A'P + PA – (PB + S)R–1(B'P + S') + Q = 0
You can use riccati( )to find the Riccati solution and gain for the
optimal regulator problem posed in Equation 6-9:
A = [0,1,0,0;21,0,0,0.8;0,0,0,1;0,0,0,-4];
B = [0,-2,0,1]';
Q = diagonal([1,0,0.1,0]);
R = B*((1e-5)\B');
Sys=system(A,B,rand(B'),[])
ys (a state space system) =
A
0
21
0
1
0
0
0
0
0
0
0.8
1
Xmath Control Design Module
6-26
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
0
0
0
0
-4
B
-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.
Example 6-9
Discrete Riccati Equation
A'PA – P – (A'PB + S)(R + B'PB)–1(B'PA + S') + Q = 0
A = [0,0,0;.3,-0.1,0;0,1,0];
Q = [2,0,0;0,0,0;0,0,10];
B = [1,0,0]';
RD = .25 ;
Sys=system(A,B,B',[],{dt=1})
Sys (a state space system) =
A
0
0
0
0
0.3
-0.1
© National Instruments Corporation
6-27
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. The general continuous Lyapunov equation is
AX + XB = –C
The special form of the continuous Lyapunov equation replaces B with A':
AX + XA' = –C (6-11)
(6-10)
These continuous Lyapunov equations have a unique solution X when
λ(A) + λj(B) ≠ 0 for any eigenvalues λi,λj, as proved in [Kai80]. This also
means that for a stable continuous-time system, X will be unique because
Xmath Control Design Module
6-28
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
You can use the Lyapunov equation to compute the state covariance matrix
of a stable system with white noise input, as illustrated in [BH75]. For a
continuous-time state-space system described by
·
x = Ax + Bu
(6-13)
and supplied with zero-mean white noise ω(t) having covariance Q:
E[ω(t)ω'(τ)] = Qδ(t – τ)
the state covariance X = E[xx'] is given by the differential Lyapunov
equation:
·
X = AX + XA' + BQB'
(6-14)
For the discrete-time system described by
xk + 1 = Axk + Buk
the white noise input covariance is defined as in the continuous case, using
a Kronecker rather than a Dirac delta.
For this case, the state covariance matrix X arises from the solution of the
discrete Lyapunov equation:
X = AXA' + BQB'
(6-15)
After you have obtained the state covariance, you can obtain the output
covariance Y easily. Whether you are using the following equation for the
continuous case:
y = Cx + Du
© National Instruments Corporation
6-29
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
or the following for the discrete case:
yk = Cxk + Duk
Y = CXC' + DQD'
(6-16)
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. As stated in [Oga70],
is described as asymptotically stable if the state values approach the
equilibrium state over time, no matter what value of x one started with.
Such systems will always satisfy the following: for any positive-definite
matrix Q, a positive definite matrix X can be found satisfying in the
(Equation 6-16) for the discrete case.
Lyapunov equations also can be used to compute system controllability
and observability grammians, which play an important role in internal
balancing and model reduction. This application will be discussed further
in the Balancing a Linear System section.
lyapunov( )
X = lyapunov(A,B,{C, discrete})
continuous-time Lyapunov equations. When called with three inputs
(A,B,C), it solves the general continuous Lyapunov equation
(Equation 6-10); when called with two inputs (A,C), it solves the special
Lyapunov equation (Equation 6-11). When called with two inputs (A,B)
and the {discrete}keyword, it solves the discrete Lyapunov equation
(refer to Equation 6-12). For examples of discrete, continuous, and special
Lyapunov equation solutions, refer to Example 6-10.
Algorithm
The algorithm for lyapunov( )uses the Schur decomposition to convert
A and B to upper triangular form, then finds the Lyapunov equation
solution a column at a time by solving. lyapunov( )warns the user if the
eigenvalues of (A + eye(A)) are close to –1, in which case singularity may
occur and cause the function to terminate. Furthermore, if any combination
Xmath Control Design Module
6-30
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Discrete Lyapunov Equation
A × X × A' + C = –X
A = [1, 2;-3,.4];
C = [-1,3;6,2];
X = lyapunov(A,C, {discrete})
X (a square matrix) =
-0.0829686 0.946549
0.390993 -0.418771
norm(A*X*A'+C -X,1)
ans (a scalar) = 2.58127e-15
Continuous Lyapunov Equation
A × X + X × B = –C
A = [1,-3;2,5];
B = [-4,3;2,1];
C = [1,3;-6,2];
X = lyapunov(A,B,C)
X (a square matrix) =
2.62963 -2.11111
-3.7037
2.22222
A*X + X*B + C;
norm(A*X + X*B + C,1)
ans (a scalar) = 0
© National Instruments Corporation
6-31
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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 -0.209028
1.12431 -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. You can specify the intensity of the
noise with the optional input covariance parameter Ucov, which defaults to
identity.
For a continuous system, the covariance of the states is given by X, where
·
X is the differential Lyapunov solution (shown in Equation 6-14) with X
equal to zero for steady-state. Thus, for a system with output Y defined by:
Y = Cx + Du
the output covariance matrix (Ycov) is expressed as:
Ycov = CXC′ + DUcovD′
The output covariance for a discrete system follows analogously,
with X being the solution to Equation 6-12 in this case. Thus, a call to
lyapunov( )forms the core of rms( ).
The diagonal terms of the covariance matrix correspond to the expected
values of the squares of the power at each output. Taking the square root of
these diagonal terms, you obtain the rms (root mean square) power at each
output. For an example of rms( )responses, refer to Example 6-11.
Xmath Control Design Module
6-32
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
Example 6-11 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]. The goal in balancing a linear
system is to find a state transformation that resolves the trade-off between
controllability and observability, returning a transformed system whose
states are equally controllable and observable. This raises the issue of
quantifying a system’s controllability or observability. You can do this
by considering the system singular values associated with the mappings
between the inputs and states, and those associated with the state-output
mappings.
These singular values can be obtained from decompositions of two
quantities referred to as the controllability and observability grammians.
These quantities are represented by Wc and Wo respectively, and defined by
the following equation for a system with an asymptotically stable A matrix.
∞
Wc = etABB'etA'dt
∫
0
∞
Wo
=
etA'C'CetAdt
∫
(6-17)
0
For continuous systems, the controllability and observability grammians
satisfy the Lyapunov equations:
AWc + WcA' + BB' = 0
A'Wo + WoA + C'C = 0
(6-18)
© National Instruments Corporation
6-33
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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:
AWcA' + BB' = Wc
A'WoA + C'C = Wo
(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]). The
condition number of Wc reflects how well conditioned the system model is
with regard to pointwise state control.
The condition number of Wo reflects the condition of the model with regard
to zero-input state-observation.
A linear transformation T of the system {A,B,C} also results in a linear
transformation of the grammians. If the state vector is transformed as
ˆ
x = Tx, the system and grammian transformations in the following
equations:
–1
ˆ
A = T AT
–1
–1
ˆ
–1
Wc = T Wc(T')
ˆ
B = T B
ˆ
ˆ
Wo = T'WoT
C = CT
ˆ
D = D
Although the poles of the system (of the eigenvalues of A) do not change
under the transformation, the singular values (eigenvalues of the
grammians) do. However, the eigenvalues of the product of the grammians
are invariant under transformation, and these are the singular values of the
system input-to-state and state-to-output maps [LHPW87].
The system is defined as being internally balanced if for some
transformation T, Wc = Wo = Σ2
2
2
ˆ
ˆ
where
Σ2 = diagonal([σ12, σ22, …, σ22])
Xmath Control Design Module
6-34
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
For model reduction problems, consider the balanced model partition as:
·
x1
A11 A12 x1
A21 A22 x2
B1
B2
=
+
u
·
x2
x1
y = [C1C2]
+ Du
s2
Σ12
0 Σ22
0
Σ2
=
with
Σ21 = diagonal(σ12....σ2k )
Σ22 = diagonal(σk2 + 1....σn2)
The essence of a balanced model reduction is that if σ2k >> σ2
,
k + 1
the input/output behavior of the states in x2 is much less important than
that of the states in x1. Eliminating the part of the model corresponding to
x2 will result in a reduced-order model which retains the most important
input-output characteristics of the original system.
balance( )
[SysB,HSV,T] = balance(Sys)
system, returning the system transformed to a balanced form as SysB. HSV
contains the second-order modes of the balanced system, or the singular
values of H, where His as defined previously.
© National Instruments Corporation
6-35
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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. It then
performs a singular-value decomposition of both grammians:
[Sc,Uc,Vc]=SVD(Wc); [So,Uo,Vo]=SVD(Wo)
and constructs the H matrix from the square roots of the singular values:
H=diagonal(sqrt(So))*Uo'*Uc*diagonal(sqrt(Sc))
The singular-value decomposition of H returns the Hankel singular
values HSV. The transformation matrix T is obtained by backsolving
and retransforming. The algorithm is given in [Moo81]. When the
transformation has been found, the balanced system matrices then can be
obtained from the original system through Equation 6-18. For an example
of how to balance a system, refer to Example 6-12.
Example 6-12 Balancing a System
Taking a hypothetical system:
A=-[1,2,3,4;0,5,6,7;0,0,8,9;0,0,0,10];
B=[0;0;0;1];
C=[1,0,0,0];
D=0;
Sys = system(A,B,C,D);
Computing the controllability and observability grammians and noting
their rather high condition numbers:
Wc=lyapunov(A,B*B');
Wo=lyapunov(A',C'*C);
condition(Wc)
ans (a scalar) = 283.029
condition(Wo)
ans (a scalar) = 1112.66
You then balance the system,
[SysB,HSV,T]=balance(Sys);
[Ab,Bb,Cb]=ABCD(SysB);
Xmath Control Design Module
6-36
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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.
Modal Form of a System
The modes of a state-space system are defined as corresponding to the
eigenvalues of the system’s A matrix. The modes of a system are distinct
from the states of a system; because a given system can be arbitrarily
transformed, the states can be arbitrarily assigned. The modes, on the other
hand, do not change from realization to realization of a given system.
The modal decomposition of a system can be obtained mathematically
through a Laplace transform, partial fraction decomposition, and eigen
decomposition as shown in [Kai80]. The key advantage of a modal
decomposition is that it provides a means by which large systems can
be represented as a parallel combination of first-order systems. In addition,
the modal decomposition of a given system representation is often better
conditioned numerically.
The modal form is particularly useful with structured dynamic systems
whose poles primarily occur as complex pairs. When a system model has
been converted to modal form, it can be reduced to focus attention on the
particular modes whose dynamics are of interest.
modal( )
[SysMod, T] = modal(Sys)
The modal( )function uses eigenvalue decomposition to find the Jordan
form of the system matrix A (all eigenvalues on the diagonal). This
approach is appropriate for models without repeated eigenvalues; modal
decomposition of a system with repeated eigenvalues is numerically
unreliable. If a system with repeated or very closely spaced eigenvalues is
passed to modal( ), a warning appears noting that the results may not be
© National Instruments Corporation
6-37
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
accurate. Given a variable Sysbuilt from the matrices {A,B,C,D}, the
modal decomposition SysModis 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.
modal( )does not accept input systems in transfer-function form, as the
concept of modes applies only to a state-variable system representation and
modes and poles are not interchangeable terms. The poles of a transfer
function always correspond to the system modes (eigenvalues of the system
A matrix).
mreduce( )
SysRed = mreduce(Sys, keep)
The mreduce( )function computes a reduced-order form of a given
system by retaining the states indicated within the vector keep. States not
specified within this vector are eliminated to obtain a lower-order model
SysRed.
mreduce( )is implemented by partitioning the state vector x into two
subvectors, x1 (states to be retained in the reduction) and x2 (states to be
eliminated in the reduction), so that:
x1
x =
x2
Similarly, the A, B, and C matrices are partitioned according to this state
partition:
A11 A12
A21 A22
B1
B2
A =
B =
C =
C1 C2
The model reductions differ for the continuous and discrete-time cases
because the updates for the states being eliminated are handled differently
in the respective differential and difference equations. In both cases, the
eliminated states are taken to be constant over time. In the continuous case,
Xmath Control Design Module
6-38
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Chapter 6
State-Space Design
the derivative of x2 is set to zero, resulting in reduced-order state equations
of the form:
–1
–1
·
x1 = (A11 – A12A22A21)x1 + (B1 – A12A22B2)u
–1
22
–1
22
y = (C1–C2A A21)x1 + (D – C2A B2)u
In the discrete case, x2k + 1 is taken to be equal to x2k so that the state
equations become:
x1k + 1 = [A11 – A12(A22 – I)–1A21]x1k
+
[B1 – A12(A22 – I)–1B2]uk
y = [C1–C2(A22 – I)–1A21]x1k
+
[D – C2(A22 – I)–1B2]uk
When using mreduce( ), remember to remove states corresponding to
in a pair—will produce a meaningless system.
More complex model reduction algorithms, which are intended to model
complete system dynamics in the absence of one of more states, are
available with the Xmath Model Reduction Module, as shown in Figure 6-8
and in Example 6-13.
Example 6-13 Model Reduction Module
A= [0.37,0.26,0.22,0.67;
0,0.52,0.63,0.20;
0,0,0.76,0.39;
0,0,0.04,0.83]
B = [0,1.7e-5,0,0.0004]'
C = [1,0,1,0]
D = 0
Sys = system(A,B,C,D,{dt = 0.2});
[SysM, T] = modal(Sys)
SysM (a state space system) =
A
0.37 0 0 0
0 0.52 0 0
0 0 0.665289 0
© National Instruments Corporation
6-39
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
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.866186 -0.668844 -0.641745
0 0.499722 -0.71998 -0.653294
0 0 -0.17991 -0.370059
0 0 0.043691 -0.15629
eig(A)
ans (a column vector) =
0.37
0.52
0.665289
0.924711
T\A*T
ans (a square matrix) =
0.37 5.55112e-17 -6.245e-17 -2.77556e-16
Xmath Control Design Module
6-40
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
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:.2:10),{keep,
legend=["Original System";"Reduced System"]})
© National Instruments Corporation
6-41
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
A
Technical References
[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]
[ChB84]
[DeS74]
A.E. Bryson and Y.C. Ho, Applied Optimal Control, Hemisphere
Publishing Corporation, Washington, D.C., 1975.
R.V. Churchill and J.W. Brown, Complex Variables and Applications,
McGraw-Hill Book Company, New York, 1984.
C.A. Desoer and J.D. Schulman, “Zeros and Poles of Matrix Transfer
Functions and Their Dynamical Interpretation,” IEEE Transactions on
Circuits and Systems, CAS-21, pages 3–8, 1974.
[EmV82]
[FPE87]
[FPW90]
[GrD86]
A. Emami-Naeini and P. Van Dooren, “Computation of Zeros of Linear
Multivariable Systems,” 18, 4, pages 415–430, 1982.
G.F. Franklin, J.D. Powell, A. Emami-Naeini, Feedback Control of
Dynamic Systems, Addison-Wesley Publishing Company, New York, 1987.
G.F. Franklin, J.D. Powell, M.L. Workman, Digital Control of Dynamic
Systems, Addison-Wesley Publishing Company, New York, 1990.
R.M. Gray and L.D. Davisson, Random Processes, A Mathematical
Approach for Engineers, Prentice-Hall, Inc., Englewood Cliffs,
New Jersey, 1986.
[HW91]
[KaB61]
P. Hsu and J. Wendlandt, “The Wedge-A Controller Design Experiment,”
Preprints, IFAC Advances in Control Education, pages 169–174, 1991.
R.E. Kalman and R. Bucy, “New Results in Linear Filtering and
Prediction,” Transactions of the American Society of Mechanical
Engineers, 83D, page 95, 1961.
[Kal60]
R.E. Kalman, “A New Approach to Linear Filtering and Prediction
page 35, 1960.
© National Instruments Corporation
A-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Appendix A
[Kai80]
[Kai81]
[Lau79]
[Lau86]
Technical References
T. Kailath, Linear Systems, Prentice-Hall, Inc., Englewood Cliffs,
New Jersey, 1980.
T. Kailath, Lectures on Wiener and Kalman Filtering, Springer-Verlag,
New York, 1981.
A.J. Laub, “A Schur method for solving algebraic Riccati equations,”
IEEE Transactions on Automatic Control, AC-24, pages 913–621, 1979.
A.J. Laub, “Efficient Calculation of Frequency Response Matrices from
State Space Models,” ACM Transactions on Mathematical Software, 12, 1,
pages 26–33, 1986.
[Leo89]
A. Leon-Garcia, Probability and Random Processes for Electrical
Engineering, Addison-Wesley Publishing Company, New York, 1989.
[LHPW87]
A.J. Laub, M.T. Heath, C.C. Paige, and R.C. Ward, “Computation of
System Balancing Transformations and Other Applications of
Simultaneous Diagonalization Algorithms,” IEEE Transactions on
Automatic Control, AC-32, 2, pages 115–121, 1987.
[Oga70]
[PLS80]
K. Ogata, Modern Control Engineering, Prentice-Hall, Inc., Englewood
Cliffs, New Jersey, 1970.
T. Pappas, A.J. Laub, and N.R. Sandell, Jr., “On the Numerical Solution of
the Discrete-Time Algebraic Riccati Equation,” IEEE Transactions on
Automatic Control, AC-25, 4, pages 631–641, 1980.
[ShH92]
[Van79]
B. Shahian and M. Hassul, Control System Design Using MATRIXx,
Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1992.
P. Van Dooren, “The computation of Kronecker’s canonical form of a
singular pencil,” Linear Algebra and Applications, 27, pages 103–141,
1979.
[Van81]
P. Van Dooren, “The generalized eigenstructure problem in linear system
theory,” IEEE Transactions on Automatic Control, AC-26, 1, pages
111–129, 1981.
Xmath Control Design Module
A-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
B
Technical Support and
Professional Services
Visit the following sections of the National Instruments Web site at
ni.comfor technical support and professional services:
•
Support—Online technical support resources at ni.com/support
include the following:
–
Self-Help Resources—For answers and solutions, visit the
award-winning National Instruments Web site for software drivers
and updates, a searchable KnowledgeBase, product manuals,
step-by-step troubleshooting wizards, thousands of example
programs, tutorials, application notes, instrument drivers, and
so on.
–
Free Technical Support—All registered users receive free Basic
Service, which includes access to hundreds of Application
Engineers worldwide in the NI Discussion Forums at
ni.com/forums. National Instruments Application Engineers
make sure every question receives an answer.
For information about other technical support options in your
area, visit ni.com/servicesor contact your local office at
ni.com/contact.
•
Training and Certification—Visit ni.com/trainingfor
self-paced training, eLearning virtual classrooms, interactive CDs,
and Certification program information. You also can register for
instructor-led, hands-on courses at locations around the world.
System Integration—If you have time constraints, limited in-house
technical resources, or other project challenges, National Instruments
Alliance Partner members can help. To learn more, call your local NI office
or visit ni.com/alliance.
If you searched ni.comand could not find the answers you need, contact
your local office or NI corporate headquarters. Phone numbers for our
the Worldwide Offices section of ni.com/niglobalto access the branch
© National Instruments Corporation
B-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
office Web sites, which provide up-to-date contact information, support
phone numbers, email addresses, and current events.
Download from Www.Somanuals.com. All Manuals Search And Download.
Index
connection, 3-1
A
abcd, 1-16, 2-8
parallel, 3-1, 3-2
series, 3-2
adjoint system, 3-3
afeedback, 3-4
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
append, 3-6
appending dynamic systems, 3-6
autocorrelation function, 5-20
analysis, 2-7
B
balance, 6-35
continuous time Riccati equation, 6-13, 6-25
controllability, 6-4, 6-35
matrix, 6-2
controllable, 6-3
bilinear transform, 2-14
Bode
default frequency range
(deffreqrange), 5-10
format, 5-8
controllable partition of a state-space system, 6-3
conventions used in the manual, iv
coordinate transformation, 6-7
corner frequency, 5-12
frequency analysis, 5-7
plots, 5-7
bode, 5-10
cross-spectral density, 5-20
C
cancel, 2-2
cascaded systems, 5-21
Cauchy’s principle, 5-15
caution, A-1
decibel gain, 5-8, 5-16
deffreqrange, 5-10
deftimerange, 4-15
delay time, td, 4-18
check, 1-5, 4-6
convert keyword, 2-12
with system objects, 2-12
choosing a sample rate, 1-18
closed-loop system
eigenvalues, 6-14
simulate performance, 1-23
combinepf, 4-9
compensator
function, 6-17
impulse, 5-20
diagnostic tools (NI resources), B-1
discrete
Riccati equation, 6-13, 6-20
system, 2-4, 4-13
analysis, 2-7
direct, in predictor form, 6-22
LQG, 6-21
checking for, 2-12
© National Instruments Corporation
I-1
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Index
converting to continuous
(example), 1-18
estimator system, 1-17
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 exponential keyword, 2-16
documentation
constant gain, 3-8
loop, 3-11
filters, 5-20
freq, to find values of a transfer function at one
frequency, 1-10
conventions used in the manual, iv
NI resources, B-1
frequency response, 1-19, 5-5
open-loop, 5-14
domain, regular, 4-11
drivers (NI resources), B-1
duality, 6-6, 6-9, 6-10
principle, 6-14
dynamic system, 6-32
appending systems, 3-6
frequency response, 5-5
improper, 3-4
G
gain margin, 5-8, 5-9, 5-12
general interconnection around a system, 3-8
impulse response, 4-13
indexing, 2-12
E
eigenvalues, 4-2, 6-6
generalized solver, 6-26
encirclements, 5-15
equilibrium state, 6-30
error covariance, 6-19
optimal, 6-18
helicopter hover problem, 1-4
ad hoc approach, 1-4
discrete formulation, 1-4, 1-18
state feedback and observer design, 1-4,
1-13
help, technical support, B-1
Xmath Control Design Module
I-2
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
Index
I
impulse, 4-13
input, 4-13
magnitude, 5-5, 5-8
makecontinuous, 2-17
verifying discretization with, 2-18
makepoly, 2-3
continuous time, 4-13
discrete time, 4-13
response, 4-13, 6-33
margin, 5-12
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
inverse, 3-3
Markov parameters, 4-13
matrix
controllability, 6-2
inputs disturbance, 6-17
observability, 6-5
Riccati equation, 6-18
transformation, 6-9
MATRIXx help, 1-3
maximum overshoot, Mp, 4-18
measurement
inverted pendulum problem, 6-21
noise, 6-16, 6-17
update, 6-19
K
MIMO systems, 2-5
minimal, 1-11, 6-8
minimal realization of a system, 6-7
Kronecker delta function, 6-29
mreduce, 6-38
L
linear
quadratic Gaussian compensation, 6-21
quadratic regulator, 6-12
systems, defining, 2-1
transformation, 6-34
logarithmic plots, 5-8
LQG
names, 2-11
default for systems, extracting, 2-11
modifying state estimate, 1-15
services, B-1
neutral stability, 5-9
nichols, 5-14
compensator, 6-21
estimator, weighting matrix, 1-21
regulator, weighting matrix, 1-21
lqgcomp, 6-23
Lyapunov
noise
continuous equation, 6-28
discrete equation, 6-29
lyapunov, 6-30
intensity matrices, 6-17
measurement, 6-16, 6-17, 6-23
process, 6-16
© National Instruments Corporation
I-3
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Index
steady-state response to white noise, 6-29
white process, 6-23
nomenclature, 1-2
numden, 2-10
period, 2-10
phase, 5-5, 5-11
f, 5-8
margin, 5-8, 5-9, 5-12
numerical integration methods, 2-14
Nyquist
rate of change, 5-12
tracking, 5-6
contour, 5-19
poleplace, 6-10
plot, 5-17
stability criterion, 5-15
nyquist, 5-16
finding feedback gains, 1-14
finding observer gains, 1-15
poles, 1-1, 2-2, 4-3
open-loop, 5-15
placement, 6-1
O
problem, 6-10
unstable open-loop, 5-15
observability, 6-4, 6-35
matrix, 6-5
cancellation, 1-11
observable, 6-6
matching, 2-15
system, 6-6
power spectral density, 5-20, 5-21
observer, 6-5, 6-16
-based controller, 1-14
gains, finding with poleplace, 1-15
open-loop
programming examples (NI resources), B-1
Q
frequency response, 5-14
poles, 5-15
operators
linear system interconnection, 3-1
overloaded, 3-1
regular domain, 4-11
regulator
optimal
estimator, 6-18
linear quadratic, 6-12
optimal, 6-12
regulator, 6-12
output
residue, 1-2, 4-8
residues, definition, 4-5
continuous-time equation, 6-25, 6-26
discrete-time equation, 6-25, 6-27
riccati, 6-26
Riccati discrete-time equation, 6-27
rise time, tr, 4-18
P
parallel connection, 3-1, 3-2
partial fraction expansion, 4-9
PDM, 4-10, 5-6
peak time, tp, 4-18
rlocus, 1-7, 5-3
Xmath Control Design Module
I-4
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
rms, 6-32
root locus, plotting, 5-3
controllable partition of, 6-3
convert to transfer function form, 1-5
with numden, 2-10
S
sample period, extract with period, 2-10
sample rate, 2-4
choosing, 1-18
modes of, 6-37
scale system output, 1-16
Schur
observable partition of, 6-6
steady state, 6-30
system response, 6-28
step, 1-9, 4-18
decomposition, 6-30
solver, 6-26
selection matrix for adding or removing
inputs, 3-9
sensitivity of states, determining, 4-17
serial, 3-1
step response, 4-18
support, technical, B-1
system, 2-6
analysis, 2-7
series connection, 3-2
settling time, ts, 4-18
system, 1-23
cascaded, 5-21
connections, 3-1
continuous, 4-13
controllability, 6-4
controllable, 6-1
singular values, 6-35
singular-value decomposition, 6-36
SISO systems, 2-3
software (NI resources), B-1
square system, 3-5
stability
discrete, 4-13
general interconnection, 3-8
impulse input to, 4-13
initial values for states, 2-7
input names, 2-7
inverse, 3-3
checking for, 2-12
criteria, 5-8
keywords, 2-7
minimal, 6-1
criterion for Nyquist, 5-15
neutral, 5-9
stair, 6-9
objects, using check with, 2-12
observability, 6-4
observable, 6-1
output names, 2-7
reformat an existing system, 2-8
square, 3-5
staircase
algorithm, 6-3, 6-6
form, 6-9
state
states names, 2-7
covariance, 6-29
names, extracting, 2-11
sensitivity, determining, 4-17
transformation, 6-33
© National Instruments Corporation
I-5
Xmath Control Design Module
Download from Www.Somanuals.com. All Manuals Search And Download.
Index
troubleshooting (NI resources), B-1
T
technical support, B-1
time
domain simulation, general, 4-10
update, 6-19
uncontrollable modes, 6-7
checking for, 1-21
coefficients form, 2-3
converted to state space before
decomposition, 2-9
creating, 2-3
formed from partial fractions, 4-9
system, 2-10
unity-gain negative feedback, 1-9
unobservable modes
unobservable uncontrollable modes, 6-7
Web resources, B-1
wedge problem, 1-20
weighting, 6-13
system models, 2-2
variable, 2-3
matrices, 6-13, 6-20
white noise, 6-23
transform
bilinear, 2-14
Fourier, 5-20
trapezoidal, 2-14
transformation, 6-9
coordinate, 6-7
zeros, 1-2, 4-3
linear, 6-34
matrix, 6-9
state, 6-33
transmission, 4-4
zeros of the transfer function, 2-2
transmission zeros, 4-4
Xmath Control Design Module
I-6
ni.com
Download from Www.Somanuals.com. All Manuals Search And Download.
|