Usingcontrol Matlab
Usingcontrol Matlab
User's Guide
R2018a
How to Contact MathWorks
Phone: 508-647-7000
v
Using Model Objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-23
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-24
Model Creation
2
Transfer Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-3
Transfer Function Representations . . . . . . . . . . . . . . . . . 2-3
Commands for Creating Transfer Functions . . . . . . . . . . 2-4
Create Transfer Function Using Numerator and
Denominator Coefficients . . . . . . . . . . . . . . . . . . . . . . 2-4
Create Transfer Function Model Using Zeros, Poles, and
Gain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-5
vi Contents
Other Model Types in Discrete Time Representations . . 2-23
vii
Frequency Response Data (FRD) Model with Time
Delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-69
viii Contents
Use Model Arrays to Create Linear Parameter-Varying
Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-114
Approximate Nonlinear Systems using LPV Models . . . 2-115
Applications of Linear Parameter-Varying Models . . . . 2-116
Data Manipulation
3
Store and Retrieve Model Data . . . . . . . . . . . . . . . . . . . . . . 3-2
Model Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-2
Specify Model Properties at Model Creation . . . . . . . . . . 3-2
Examine and Change Properties of an Existing Model . . . 3-3
ix
Configure Display Format of Transfer Function in
Factorized Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-18
Model Interconnections
4
Why Interconnect Models? . . . . . . . . . . . . . . . . . . . . . . . . . 4-2
x Contents
Model Transformation
5
Conversion Between Model Types . . . . . . . . . . . . . . . . . . . 5-2
Explicit Conversion Between Model Types . . . . . . . . . . . 5-2
Automatic Conversion Between Model Types . . . . . . . . . 5-2
Recommended Working Representation . . . . . . . . . . . . . 5-3
xi
Model Simplification
6
Model Reduction Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-2
When to Reduce Model Order . . . . . . . . . . . . . . . . . . . . . 6-2
Choosing a Model Reduction Method . . . . . . . . . . . . . . . 6-4
xii Contents
Linear Analysis
xiii
Frequency Domain Analysis
8
Frequency-Domain Responses . . . . . . . . . . . . . . . . . . . . . . 8-2
Sensitivity Analysis
9
Model Array with Single Parameter Variation . . . . . . . . . . 9-2
xiv Contents
Passivity and Conic Sectors
10
About Passivity and Passivity Indices . . . . . . . . . . . . . . . 10-2
Control Design
xv
Classical Control Design
12
Choosing a Control Design Approach . . . . . . . . . . . . . . . 12-2
xvi Contents
Edit Compensator Dynamics . . . . . . . . . . . . . . . . . . . . . . 12-93
Compensator Editor . . . . . . . . . . . . . . . . . . . . . . . . . . 12-93
Graphical Compensator Editing . . . . . . . . . . . . . . . . . 12-96
Poles and Zeros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12-97
Lead and Lag Networks . . . . . . . . . . . . . . . . . . . . . . . 12-97
Notch Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12-98
xvii
State-Space Control Design
13
Extended and Unscented Kalman Filter Algorithms for
Online State Estimation . . . . . . . . . . . . . . . . . . . . . . . . 13-2
Extended Kalman Filter Algorithm . . . . . . . . . . . . . . . . 13-2
Unscented Kalman Filter Algorithm . . . . . . . . . . . . . . . 13-5
xviii Contents
Choosing an Automated Tuning Approach . . . . . . . . . . . 14-5
xix
Quick Loop Tuning of Feedback Loops in Control System
Tuner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-48
xx Contents
Tips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-84
Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-85
xxi
Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-114
Sensitivity Evaluation . . . . . . . . . . . . . . . . . . . . . . . . 14-115
Sensitivity Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-116
Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-116
Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-117
xxii Contents
Desired Loop Shape . . . . . . . . . . . . . . . . . . . . . . . . . 14-144
Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-144
Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-146
xxiii
Controller Poles Goal . . . . . . . . . . . . . . . . . . . . . . . . . . 14-176
Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-176
Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-176
Constrain Dynamics of Tuned Block . . . . . . . . . . . . . 14-177
Keep Poles Inside the Following Region . . . . . . . . . . 14-177
Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-178
xxiv Contents
Speed Up Tuning with Parallel Computing Toolbox
Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14-225
Loop-Shaping Design
15
Structure of Control System for Tuning With
looptune . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15-2
Gain-Scheduled Controllers
16
Gain Scheduling Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . 16-2
Gain Scheduling in Simulink . . . . . . . . . . . . . . . . . . . . . 16-2
Tune Gain Schedules . . . . . . . . . . . . . . . . . . . . . . . . . . 16-3
xxv
Tune Gain Schedules in Simulink . . . . . . . . . . . . . . . . . 16-15
Workflow for Tuning Gain Schedules . . . . . . . . . . . . . 16-15
xxvi Contents
Control System Tuning Examples
17
Tuning Multiloop Control Systems . . . . . . . . . . . . . . . . . . 17-3
xxvii
Multiloop Control of a Helicopter . . . . . . . . . . . . . . . . 17-216
Customization
Preliminaries
18
Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18-2
xxviii Contents
Setting Toolbox Preferences
19
Toolbox Preferences Editor . . . . . . . . . . . . . . . . . . . . . . . . 19-2
Overview of the Toolbox Preferences Editor . . . . . . . . . 19-2
Opening the Toolbox Preferences Editor . . . . . . . . . . . . 19-2
Units Pane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19-3
Style Pane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19-5
Options Pane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19-6
Control System Designer Pane . . . . . . . . . . . . . . . . . . . 19-7
xxix
Customizing Response Plots Using Plot Tools . . . . . . . 21-24
Properties You Can Customize Using Plot Tools . . . . . . 21-24
Opening and Working with Plot Tools . . . . . . . . . . . . . 21-24
Example of Changing Line Color Using Plot Tools . . . . 21-25
xxx Contents
Reliable Computations
23
Scaling State-Space Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 23-2
Why Scaling Is Important . . . . . . . . . . . . . . . . . . . . . . . . . . . 23-2
When to Scale Your Model . . . . . . . . . . . . . . . . . . . . . . . . . . 23-2
Manually Scale Your Model . . . . . . . . . . . . . . . . . . . . . . . . . . 23-3
xxxi
Linear System Modeling
33
1
Other model attributes stored as model data include time units, names for the model
inputs or outputs, and time delays. For more information about setting and retrieving
model attributes, see “Model Attributes”.
1-2
See Also
Note All model objects are MATLAB® objects, but working with them does not require a
background in object-oriented programming. To learn more about objects and object
syntax, see “Role of Classes in MATLAB” (MATLAB).
See Also
More About
• “Control System Modeling with Model Objects” on page 1-4
• “Types of Model Objects” on page 1-7
1-3
1 Linear System Model Objects
For example, the following control system contains a prefilter F, a plant G, and a
controller C, arranged in a single-loop configuration. The model also includes a
representation of sensor dynamics, S.
You can represent each of the components as a model object. You do not need to use the
same type of model object for each component. For example, represent the plant G as a
zero-pole-gain (zpk) model with a double pole at s = -1; C as a PID controller, and F and
S as transfer functions:
G = zpk([],[-1,-1],1);
C = pid(2,1.3,0.3,0.5);
S = tf(5,[1 4]);
F = tf(1,[1 1]);
You can then combine these elements build models that represent your control system or
the control system as a whole. For example, create the open-loop response SGC:
open_loop = S*G*C;
To build a model of the unfiltered closed-loop response, use the feedback command:
T = feedback(G*C,S);
To model the entire closed-loop system response from r to y, combine T with the filter
transfer function:
Try = T*F;
1-4
Control System Modeling with Model Objects
The results open_loop, T, and Try are also linear model objects. You can operate on
them with Control System Toolbox™ control design and analysis commands. For example,
plot the step response of the entire system:
stepplot(Try)
When you combine Numeric LTI models, the resulting Numeric LTI model represents the
aggregate system. The resulting model does not retain the original data from the
combined components. For example, T does not separately keep track of the dynamics of
the components G, C, and S that are combined to create T.
1-5
1 Linear System Model Objects
See Also
feedback
Related Examples
• “Numeric Model of SISO Feedback Loop” on page 4-6
• “Multi-Loop Control System” on page 4-10
• “MIMO Control System” on page 4-19
More About
• “Types of Model Objects” on page 1-7
1-6
Types of Model Objects
1-7
1 Linear System Model Objects
The diagram illustrates the following two overlapping broad classifications of model
object types:
1-8
See Also
• Dynamic System Models vs. Static Models — In general, Dynamic System Models
represent systems that have internal dynamics, while Static Models represent static
input/output relationships.
• Numeric Models vs. Generalized Models — Numeric Models are the basic numeric
representation of linear systems with fixed coefficients. Generalized Models represent
systems with tunable or uncertain components.
See Also
More About
• “What Are Model Objects?” on page 1-2
• “Dynamic System Models” on page 1-10
• “Static Models” on page 1-12
• “Numeric Models” on page 1-13
• “Generalized Models” on page 1-16
1-9
1 Linear System Model Objects
Most commands for analyzing linear systems, such as bode, margin, and
linearSystemAnalyzer, work on most Dynamic System Model objects. For Generalized
Models, analysis commands use the current value of tunable parameters and the nominal
value of uncertain parameters. Commands that generate response plots display random
samples of uncertain models.
1-10
See Also
See Also
More About
• “Numeric Linear Time Invariant (LTI) Models” on page 1-13
• “Identified LTI Models” on page 1-14
• “Identified Nonlinear Models” on page 1-14
• “Generalized and Uncertain LTI Models” on page 1-16
• “Control Design Blocks” on page 1-16
1-11
1 Linear System Model Objects
Static Models
Static Models represent static input/output relationships and generalize the notions of
matrix and numeric array to parametric or uncertain arrays. You can use static models to
create parametric or uncertain expressions, and to construct Generalized LTI models
whose coefficients are parametric or uncertain expressions. The Static Models family
includes:
For more information about using these objects to create parametric models, see “Models
with Tunable Coefficients” on page 1-19. For information about creating uncertain static
models, see “Uncertain Real Parameters” (Robust Control Toolbox) and “Uncertain
Matrices” (Robust Control Toolbox).
1-12
Numeric Models
Numeric Models
Numeric Linear Time Invariant (LTI) Models
Numeric LTI models are the basic numeric representation of linear systems or
components of linear systems. Use numeric LTI models for modeling dynamic
components, such as transfer functions or state-space models, whose coefficients are
fixed, numeric values. You can use numeric LTI models for linear analysis or control
design tasks.
The following table summarizes the available types of numeric LTI models.
You can use Numeric LTI models to represent block diagram components such as plant or
sensor dynamics. By connecting Numeric LTI models together, you can derive Numeric
LTI models of block diagrams. Use Numeric LTI models for most modeling, analysis, and
control design tasks, including:
1-13
1 Linear System Model Objects
• Analyzing linear system dynamics using analysis commands such as bode, step, or
impulse.
• Designing controllers for linear systems using the Control System Designer app or
the PID Tuner GUI.
• Designing controllers using control design commands such as pidtune, rlocus, or
lqr/lqg.
The following table summarizes the available types of identified LTI models.
The following table summarizes the available types of identified nonlinear models.
1-14
Numeric Models
1-15
1 Linear System Model Objects
Generalized Models
Uncertain LTI Models are a special type of Generalized LTI model that include uncertain
coefficients but not tunable coefficients. For more information about using uncertain
models, see “Uncertain State-Space Models” (Robust Control Toolbox) and “Create
Uncertain Frequency Response Data Models” (Robust Control Toolbox).
1-16
Generalized Models
Tunable Control Design Blocks include tunable parameter objects as well as tunable
linear models with predefined structure. For more information about using tunable
Control Design Blocks, see “Models with Tunable Coefficients” on page 1-19.
If you have Robust Control Toolbox software, you can use uncertain Control Design Blocks
to model uncertain parameters or uncertain system dynamics. For more information
about using uncertain blocks, see “Uncertain LTI Dynamics Elements” (Robust Control
Toolbox), “Uncertain Real Parameters” (Robust Control Toolbox), and “Uncertain Complex
Parameters and Matrices” (Robust Control Toolbox).
The following tables summarize the available types of Control Design Blocks.
1-17
1 Linear System Model Objects
Generalized Matrices
Generalized Matrices extend the notion of numeric matrices to matrices that include
tunable or uncertain values.
If you have Robust Control Toolbox software, you can create uncertain matrices by
building rational expressions involving uncertain parameters such as ureal or
ucomplex.
For more information about generalized matrices and their applications, see “Models with
Tunable Coefficients” on page 1-19.
1-18
Models with Tunable Coefficients
• Model a tunable (or parametric) component of a control system, such as a tunable low-
pass filter.
• Model a control system that contains both:
You can use tunable Generalized LTI models for parameter studies. For an example, see
“Study Parameter Variation by Sampling Tunable Model” on page 9-9. You can also use
tunable Generalized LTI models for tuning fixed control structures using tuning
commands such as systune or the Control System Tuner app. See “Multiloop,
Multiobjective Tuning”.
To create tunable components with a specific custom structure that is not covered by the
Control Design Blocks:
1 Use the tunable real parameter realp or the generalized matrix genmat to
represent the tunable coefficients of your component.
2 Use the resulting realp or genmat objects as inputs to tf or ss to model the
component. The result is a generalized state-space (genss) model of the component.
1-19
1 Linear System Model Objects
1 Model the nontunable components of your system using numeric LTI models on page
1-13.
2 Model each tunable component using Control Design Blocks or expressions involving
such blocks. See “Modeling Tunable Components” on page 1-19.
3 Use model interconnection commands such as series, parallel or connect, or
the arithmetic operators +, -, *, /, \, and ^, to combine all the components of your
system.
For an example of constructing a genss model of a control system with both fixed and
tunable components, see “Control System with Tunable Components” on page 2-83.
1-20
Models with Tunable Coefficients
w z
H
u y
B1 0 ... 0
.
0 B2 ..
.. ..
. . 0
0 ... 0 BN
B
w and z represent the inputs and outputs of the Generalized model.
H represents all portions of the Generalized model that have fixed (non-parametric)
coefficients. H is:
B represents the parametric components of the Generalized model, which are the Control
Design Blocks B1, . . . , BN. The Blocks property of the Generalized model stores a list of
the names of these blocks. If the Generalized model has blocks that occur multiple times
in B1, . . . , BN, these are only listed once in the Blocks property.
To access the internal representation of a Generalized model, including H and B, use the
getLFTModel command.
1-21
1 Linear System Model Objects
This Standard Form can represent any control structure. To understand why, consider the
control structure as an aggregation of fixed-coefficient elements interacting with the
parametric elements:
external
inputs H
w
Fixed system
components
(actuators, external
sensors, etc.) outputs
y1
uN z
B1 BN
yN
u1
B2
...
y2
u2
u := [u1 ,…, uN ]
y := [ y1 ,… , yN ] ,
and group the tunable control elements B1, . . . , BN into the block-diagonal configuration
C. P includes all the fixed components of the control architecture—actuators, sensors, and
other nontunable elements—and their interconnections.
1-22
Using Model Objects
• Attach additional information to the model using model attributes (properties). See
“Model Attributes”.
• Manipulate the model using arithmetic and model interconnection operations. See
“Model Interconnection”.
• Analyze the model response using commands such as bode and step. See “Linear
Analysis”.
• Perform parameter studies using model arrays. See “Model Arrays”.
• Design compensators. You can:
1-23
1 Linear System Model Objects
References
[1] Dorf, R.C. and R.H. Bishop, Modern Control Systems, Addison-Wesley, Menlo Park, CA,
1998.
1-24
2
Model Creation
2-2
Transfer Functions
Transfer Functions
Transfer Function Representations
Control System Toolbox software supports transfer functions that are continuous-time or
discrete-time, and SISO or MIMO. You can also have time delays in your transfer function
representation.
N ( s)
G ( s) = ,
D( s )
of polynomials N(s) and D(s), called the numerator and denominator polynomials,
respectively.
You can represent linear systems as transfer functions in polynomial or factorized (zero-
pole-gain) form. For example, the polynomial-form transfer function:
s2 - 3 s - 4
G ( s) =
s2 + 5 s + 6
( s + 1) (s - 4)
G ( s) = .
( s + 2)( s + 3)
The tf model object represents transfer functions in polynomial form. The zpk model
object represents transfer functions in factorized form.
MIMO transfer functions are arrays of SISO transfer functions. For example:
Ès - 3˘
Ís + 4˙
G ( s) = Í ˙
Í s + 1˙
ÍÎ s + 2 ˙˚
2-3
2 Model Creation
Command Description
tf Create tf objects representing continuous-time or discrete-
time transfer functions in polynomial form.
zpk Create zpk objects representing continuous-time or discrete-
time transfer functions in zero-pole-gain (factorized) form.
filt Create tf objects representing discrete-time transfer
functions using digital signal processing (DSP) convention.
s
Create the transfer function G ( s ) = :
2
s + 3s + 2
num = [1 0];
den = [1 3 2];
G = tf(num,den);
num and den are the numerator and denominator polynomial coefficients in descending
powers of s. For example, den = [1 3 2] represents the denominator polynomial
s2 + 3s + 2.
Tip Alternatively, you can specify the transfer function G(s) as an expression in s:
s = tf('s');
2 Specify G(s) as a ratio of polynomials in s.
2-4
See Also
s
Create the factored transfer function G ( s ) = 5 :
( s + 1 + i) ( s + 1 - i) ( s + 2)
Z = [0];
P = [-1-1i -1+1i -2];
K = 5;
G = zpk(Z,P,K);
Z and P are the zeros and poles (the roots of the numerator and denominator,
respectively). K is the gain of the factored form. For example, G(s) has a real pole at s = –
2 and a pair of complex poles at s = –1 ± i. The vector P = [-1-1i -1+1i -2] specifies
these pole locations.
G is a zpk model object, which is a data container for representing transfer functions in
zero-pole-gain (factorized) form.
See Also
filt | tf | zpk
Related Examples
• “MIMO Transfer Functions” on page 2-28
• “State-Space Models” on page 2-6
• “Discrete-Time Numeric Models” on page 2-23
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-5
2 Model Creation
State-Space Models
State-Space Model Representations
State-space models rely on linear differential equations or difference equations to
describe system dynamics. Control System Toolbox software supports SISO or MIMO
state-space models in continuous or discrete time. State-space models can include time
delays. You can represent state-space models in either explicit or descriptor (implicit)
form.
dx
= Ax + Bu
dt
y = Cx + Du
where x is the state vector. u is the input vector, and y is the output vector. A, B, C, and D
are the state-space matrices that express the system dynamics.
x [ n + 1] = Ax[ n] + Bu [ n]
y [ n] = Cx [ n ] + Du [ n]
where the vectors x[n], u[n], and y[n] are the state, input, and output vectors for the nth
sample.
2-6
State-Space Models
dx
E = Ax + Bu
dt
y = Cx + Du
where x is the state vector. u is the input vector, and y is the output vector. A, B, C, D, and
E are the state-space matrices.
Command Description
ss Create explicit state-space model.
dss Create descriptor (implicit) state-space model.
delayss Create state-space models with specified time delays.
dx
= Ax + Bu
dt
y = Cx + Du
where the state variables are the angular position θ and angular velocity dθ/dt:
Èq ˘
x = Í dq ˙ ,
Í ˙
ÎÍ dt ˙˚
2-7
2 Model Creation
u is the electric current, the output y is the angular velocity, and the state-space matrices
are:
È0 1˘ È0 ˘
A=Í ˙, B = Í ˙ , C = [ 0 1] , D = [ 0 ].
Î -5 -2 ˚ Î3 ˚
A = [0 1;-5 -2];
B = [0;3];
C = [0 1];
D = 0;
sys = ss(A,B,C,D);
sys is an ss model object, which is a data container for representing state-space models.
dx
E = Ax + Bu
dt
y = Cx + Du
use dss. This command creates a ss model with a nonempty E matrix, also called a
descriptor state-space model. See “MIMO Descriptor State-Space Models” on page 2-32
for an example.
See Also
delayss | dss | ss
Related Examples
• “MIMO State-Space Models” on page 2-31
• “Transfer Functions” on page 2-3
• “Discrete-Time Numeric Models” on page 2-23
2-8
See Also
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-9
2 Model Creation
For example, suppose you measure frequency response data for the SISO system you
want to model. You can measure such data by driving the system with a sine wave at a set
of frequencies ω1, ω2, ,...,ωn, as shown:
At steady state, the measured response yi(t) to the driving signal at each frequency ωi
takes the following form:
The measurement yields the complex frequency response G at each input frequency:
You can do most frequency-domain analysis tasks on frd models, but you cannot perform
time-domain simulations with them. For information on frequency response analysis of
linear systems, see Chapter 8 of [1].
Command Description
frd Create frd objects from frequency response data.
2-10
Frequency Response Data (FRD) Models
Command Description
frestimate Create frd objects by estimating the frequency response of a
Simulink® model. This approach requires Simulink Control
Design™ software. See “Frequency Response Estimation”
(Simulink Control Design) for more information.
load AnalyzerData
This command loads the data into the MATLAB workspace as the column vectors
freq and resp. The variables freq and resp contain 256 test frequencies and the
corresponding complex-valued frequency response points, respectively.
sys = frd(resp,freq);
sys is an frd model object, which is a data container for representing frequency
response data.
You can use frd models with many frequency-domain analysis commands. For example,
visualize the frequency response data using bode.
2-11
2 Model Creation
Tip By default, the frd command assumes that the frequencies are in radians/second. To
specify different frequency units, use the TimeUnit and FrequencyUnit properties of
the frd model object. For example:
sys = frd(resp,freq,'TimeUnit','min','FrequencyUnit','rad/TimeUnit')
See Also
frd | frestimate
Related Examples
• “MIMO Frequency Response Data Models” on page 2-37
• “Discrete-Time Numeric Models” on page 2-23
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-12
Proportional-Integral-Derivative (PID) Controllers
Form Formula
Parallel (pid object) Ki Kd s
C = Kp + + ,
s Tf s + 1
where:
• Kp = proportional gain
• Ki = integrator gain
• Kd = derivative gain
• Tf = derivative filter time
Standard (pidstd Ê ˆ
object) Á 1 Td s ˜
C = K p Á1 + + ˜,
Á Ti s Td
Á s + 1 ˜˜
Ë N ¯
where:
• Kp = proportional gain
• Ti = integrator time
• Td = derivative time
• N = derivative filter divisor
2-13
2 Model Creation
Use a controller form that is convenient for your application. For example, if you want to
express the integrator and derivative actions in terms of time constants, use standard
form.
26 .2 4.3 s
Create the following parallel-form PID controller: C = 29 .5 + - .
s 0 .06 s + 1
Kp = 29.5;
Ki = 26.2;
Kd = 4.3;
Tf = 0.06;
C = pid(Kp,Ki,Kd,Tf)
C is a pid model object, which is a data container for representing parallel-form PID
controllers. For more examples of how to create PID controllers, see the pid reference
page.
Ê ˆ
Á 1 0 .15 s ˜
Create the following standard-form PID controller: C = 29 .5 Á 1 + + ˜.
ÁÁ 1.13s 0 .15 s + 1 ˜
˜
Kp = 29.5; Ë 2.3 ¯
Ti = 1.13;
Td = 0.15;
N = 2.3;
C = pidstd(Kp,Ti,Td,N)
2-14
See Also
C is a pidstd model object, which is a data container for representing standard-form PID
controllers. For more examples of how to create standard-form PID controllers, see the
pidstd reference page.
See Also
pid | pidTuner | pidstd | pidtune
Related Examples
• “Transfer Functions” on page 2-3
• “Discrete-Time Proportional-Integral-Derivative (PID) Controllers” on page 2-24
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-15
2 Model Creation
You can represent PID controllers using the specialized model objects pid2 and
pidstd2. This topic describes the representation of 2-DOF PID controllers in MATLAB.
For information about automatic PID controller tuning, see “PID Controller Tuning”.
The relationship between the 2-DOF controller’s output (u) and its two inputs (r and y)
can be represented in either parallel or standard form. The two forms differ in the
parameters used to express the proportional, integral, and derivative actions of the
controller, as expressed in the following table.
2-16
Two-Degree-of-Freedom PID Controllers
Form Formula
Parallel (pid2 object) Ki K s
u = K p ( br - y ) + ( r - y) + d ( cr - y) .
s Tf s + 1
In this representation:
• Kp = proportional gain
• Ki = integrator gain
• Kd = derivative gain
• Tf = derivative filter time
• b = setpoint weight on proportional term
• c = setpoint weight on derivative term
Standard (pidstd2 object) È ˘
Í 1 Td s ˙
u = K p Í( br - y ) + (r - y) + T ( cr - y ) ˙ .
Í Ti s d ˙
s +1
ÍÎ N ˙˚
In this representation:
• Kp = proportional gain
• Ti = integrator time
• Td = derivative time
• N = derivative filter divisor
• b = setpoint weight on proportional term
• c = setpoint weight on derivative term
Use a controller form that is convenient for your application. For instance, if you want to
express the integrator and derivative actions in terms of time constants, use standard
form. For examples showing how to create parallel-form and standard-form controllers,
see the pid2 and pidstd2 reference pages, respectively.
2-17
2 Model Creation
Each of the components Cr(s) and Cy(s) is a PID controller, with different weights on the
proportional and derivative terms. For example, in continuous time, these components are
given by:
Ki cK d s
Cr ( s ) = bK p + + ,
s Tf s + 1
È K Kd s ˘
Cy ( s ) = - Í K p + i + ˙.
ÍÎ s Tf s + 1 ˙˚
You can access these components by converting the PID controller into a two-input, one-
output transfer function. For example, suppose that C2 is a 2-DOF PID controller, stored
as a pid2 object.
C2tf = tf(C2);
Cr = C2tf(1);
Cy = C2tf(2);
Cr(s) is the transfer function from the first input of C2 to the output. Similarly, Cy(s) is the
transfer function from the second input of C2 to the output.
2-18
Two-Degree-of-Freedom PID Controllers
Suppose that G is a dynamic system model, such as a zpk model, representing the plant.
Build the closed-loop transfer function from r to y. Note that the Cy(s) loop has positive
feedback, by the definition of Cy(s).
T = Cr*feedback(G,Cy,+1)
G.InputName = 'u';
G.OutputName = 'y';
C2.Inputname = {'r','y'};
C2.OutputName = 'u';
T = connect(G,C2,'r','y');
There are other configurations in which you can decompose a 2-DOF PID controller into
SISO components. For particular choices of C(s) and X(s), each of the following
configurations is equivalent to the 2-DOF architecture with C2(s). You can obtain C(s) and
X(s) for each of these configurations using the getComponents command.
Feedforward
For a continuous-time, parallel-form 2-DOF PID controller, the components are given by:
2-19
2 Model Creation
Ki Kd s
C ( s) = K p + + ,
s Tf s + 1
( c - 1 ) Kd s
X ( s) = ( b - 1 ) K p + .
Tf s + 1
[C,X] = getComponents(C2,'feedforward');
The following command constructs the closed-loop system from r to y for the feedforward
configuration.
T = G*(C+X)*feedback(1,G*C);
Feedback
For a continuous-time, parallel-form 2-DOF PID controller, the components are given by:
Ki cK d s
C ( s) = bK p + + ,
s Tf s + 1
(1 - c ) Kd s
X ( s) = (1 - b ) K p + .
Tf s + 1
2-20
Two-Degree-of-Freedom PID Controllers
[C,X] = getComponents(C2,'feedback');
The following command constructs the closed-loop system from r to y for the feedback
configuration.
T = G*C*feedback(1,G*(C+X));
Filter
In the filter configuration, the 2-DOF PID controller is decomposed into a conventional
SISO PID controller and a prefilter on the reference signal.
For a continuous-time, parallel-form 2-DOF PID controller, the components are given by:
Ki Kd s
C ( s) = K p + + ,
s Tf s + 1
X ( s) =
(bK pTf + cK d ) s2 + ( bK p + K iTf ) s + K i .
( K pTf + Kd ) s2 + ( K p + KiTf ) s + Ki
The filter X(s) can also be expressed as the ratio: –[Cr(s)/Cy(s)].
The following command constructs the closed-loop system from r to y for the filter
configuration.
T = X*feedback(G*C,1);
For an example illustrating the decomposition of a 2-DOF PID controller into these
configurations, see “Decompose a 2-DOF PID Controller into SISO Components” on page
5-8.
2-21
2 Model Creation
See Also
getComponents | pid2 | pidTuner | pidstd2 | pidtune
Related Examples
• “Discrete-Time Proportional-Integral-Derivative (PID) Controllers” on page 2-24
• “Proportional-Integral-Derivative (PID) Controllers” on page 2-13
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-22
Discrete-Time Numeric Models
z
Create the transfer function G ( z ) = with a sample time of 0.1 s.
2
z - 2z - 6
num = [1 0];
den = [1 -2 -6];
Ts = 0.1;
G = tf(num,den,Ts)
num and den are the numerator and denominator polynomial coefficients in descending
powers of z. G is a tf model object.
The sample time is stored in the Ts property of G. Access the sample time Ts, using dot
notation:
G.Ts
See Also
frd | ss | tf | zpk
More About
• “What Are Model Objects?” on page 1-2
2-23
2 Model Creation
Form Formula
Parallel (pid) Kd
C = K p + Ki IF ( z) + ,
Tf + DF ( z )
where:
• Kp = proportional gain
• Ki = integrator gain
• Kd = derivative gain
• Tf = derivative filter time
Standard (pidstd) Ê ˆ
Á 1 Td ˜
C = K p Á 1 + IF ( z) + ˜,
Á Ti Td
Á + DF ( ) ˜˜
z
Ë N ¯
where:
• Kp = proportional gain
• Ti = integrator time
• Td = derivative time
• N = derivative filter divisor
2-24
Discrete-Time Proportional-Integral-Derivative (PID) Controllers
Form Formula
2-DOF Parallel (pid2) The relationship between the 2-DOF controller’s output (u) and
its two inputs (r and y) is:
Kd
u = K p ( br - y ) + Ki IF ( z ) ( r - y ) + ( cr - y) .
Tf + DF ( z )
In this representation:
• Kp = proportional gain
• Ki = integrator gain
• Kd = derivative gain
• Tf = derivative filter time
• b = setpoint weight on proportional term
• c = setpoint weight on derivative term
2-DOF Standard È ˘
(pidstd2 object) Í 1 Td ˙
u = K p Í( br - y ) + IF ( z) ( r - y) + ( cr - y) ˙ .
Í Ti Td ˙
+ DF ( z )
ÍÎ N ˙˚
In this representation:
• Kp = proportional gain
• Ti = integrator time
• Td = derivative time
• N = derivative filter divisor
• b = setpoint weight on proportional term
• c = setpoint weight on derivative term
In all of these expressions, IF(z) and DF(z) are the discrete integrator formulas for the
integrator and derivative filter, respectively. Use the IFormula and DFormula properties
of the controller objects to set the IF(z) and DF(z) formulas. The next table shows
available formulas for IF(z) and DF(z). Ts is the sample time.
2-25
2 Model Creation
If you do not specify a value for IFormula, DFormula, or both when you create the
controller object, ForwardEuler is used by default. For more information about setting
and changing the discrete integrator formulas, see the reference pages for the controller
objects, pid, pidstd, pid2, and pidstd2.
Ts z + 1 Tz
This command creates a pidstd model with IF ( z) = and DF ( z ) = s .
2 z -1 z -1
You can set the discrete integrator formulas for a parallel-form controller in the same way,
using pid.
2-26
See Also
b = 0.5;
c = 0;
Ts = 0.1;
C2 = pidstd2(Kp,Ti,Td,N,b,c,Ts,'IFormula','Trapezoidal')
C2 =
1 Ts*(z+1)
u = Kp * [(b*r-y) + ---- * -------- * (r-y)]
Ti 2*(z-1)
Setting Td = 0 specifies a PI controller with no derivative term. As the display shows, the
values of N and c are not used in this controller. The display also shows that the
trapezoidal formula is used for the integrator.
See Also
pid | pid2 | pidstd | pidstd2
Related Examples
• “Proportional-Integral-Derivative (PID) Controllers” on page 2-13
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-27
2 Model Creation
È s-1 ˘
Í s+1 ˙
H ( s) = Í ˙.
Í s+ 2 ˙
ÍÎ s2 + 4 s + 5 ˙˚
You can specify H(s) by concatenation of its SISO entries. For instance,
or, equivalently,
s = tf('s')
h11 = (s-1)/(s+1);
h21 = (s+2)/(s^2+4*s+5);
H = [h11; h21]
This syntax mimics standard matrix concatenation and tends to be easier and more
readable for MIMO systems with many inputs and/or outputs.
Tip Use zpk instead of tf to create MIMO transfer functions in factorized form on page
2-5.
2-28
See Also
For example, for the rational transfer matrix H(s), the two cell arrays N and D should
contain the row-vector representations of the polynomial entries of
È s-1 ˘ È s +1 ˘
N ( s) = Í ˙, D( s) = Í ˙.
Îs + 2˚ Î s2 + 4 s + 5 ˚
s + 2
#2: -------------
s^2 + 4 s + 5
Notice that both N and D have the same dimensions as H. For a general MIMO transfer
matrix H(s), the cell array entries N{i,j} and D{i,j} should be row-vector
representations of the numerator and denominator of Hij(s), the ijth entry of the transfer
matrix H(s).
See Also
tf | zpk
Related Examples
• “Transfer Functions” on page 2-3
2-29
2 Model Creation
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-30
MIMO State-Space Models
# of columns =
# of inputs
A B
# of rows =
# of outputs C D
In this example, you create a state-space model for a rotating body with inertia tensor J,
damping force F, and three axes of rotation, related as:
dw
J + Fw = T
dt
y = w.
The system input T is the driving torque. The output y is the vector of angular velocities of
the rotating body.
2-31
2 Model Creation
dx
= Ax + Bu
dt
y = Cx + Du
rewrite it as:
dw
= - J -1 Fw + J -1T
dt
y = w.
A = - J -1 F, B = J -1 , C = I , D = 0.
These commands assume that J is the inertia tensor of a cube rotating about its corner,
and the damping force has magnitude 0.2.
sys_mimo is an ss model.
This example uses the same rotating-body system shown in “MIMO Explicit State-Space
Models” on page 2-31, where you inverted the inertia matrix J to obtain the value of the B
matrix. If J is poorly-conditioned for inversion, you can instead use a descriptor (implicit)
state-space model. A descriptor (implicit) state-space model is of the form:
2-32
MIMO State-Space Models
dx
E = Ax + Bu
dt
y = Cx + Du
Create a state-space model for a rotating body with inertia tensor J, damping force F, and
three axes of rotation, related as:
dw
J + Fw = T
dt
y = w.
The system input T is the driving torque. The output y is the vector of angular velocities of
the rotating body. You can write this system as a descriptor state-space model having the
following state-space matrices:
A = - F, B = I, C = I, D = 0, E = J.
These commands assume that J is the inertia tensor of a cube rotating about its corner,
and the damping force has magnitude 0.2.
The jet model during cruise flight at MACH = 0.8 and H = 40,000 ft. is
2-33
2 Model Creation
B = [ 0.0073 0
-0.4750 0.0077
0.1530 0.1430
0 0];
C = [0 1 0 0
0 0 0 1];
D = [0 0
0 0];
Use the following commands to specify this state-space model as an LTI object and attach
names to the states, inputs, and outputs.
sys_mimo = ss(A,B,C,D,'statename',states,...
'inputname',inputs,...
'outputname',outputs);
sys_mimo
a =
beta yaw roll phi
beta -0.0558 -0.9968 0.0802 0.0415
yaw 0.598 -0.115 -0.0318 0
roll -3.05 0.388 -0.465 0
phi 0 0.0805 1 0
b =
rudder aileron
beta 0.0073 0
yaw -0.475 0.0077
roll 0.153 0.143
2-34
MIMO State-Space Models
phi 0 0
c =
beta yaw roll phi
yaw rate 0 1 0 0
bank angle 0 0 0 1
d =
rudder aileron
yaw rate 0 0
bank angle 0 0
Continuous-time model.
The model has two inputs and two outputs. The units are radians for beta (sideslip angle)
and phi (bank angle) and radians/sec for yaw (yaw rate) and roll (roll rate). The rudder
and aileron deflections are in degrees.
tf(sys_mimo)
2-35
2 Model Creation
See Also
ss
Related Examples
• “State-Space Models” on page 2-6
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-36
MIMO Frequency Response Data Models
Frequency response data for a MIMO system includes a vector of complex response data
for each of the input/output (I/O) pair of the system. Thus, if you measure the frequency
response of each I/O pair of your system at a set of test frequencies, you can use the data
to create a frequency response model:
1 Load frequency response data in AnalyzerDataMIMO.mat.
load AnalyzerDataMIMO H11 H12 H21 H22 freq
This command loads the data into the MATLAB workspace as five column vectors
H11, H12, H21, H22, and freq. The vector freq contains 100 test frequencies. The
other four vectors contain the corresponding complex-valued frequency response of
each I/O pair of a two-input, two-output system.
The dimensions of Hresp are the number of outputs, number of inputs, and the
number of frequencies for which there is response data. Hresp(i,j,:) contains the
frequency response from input j to output i.
3 Create a frequency-response model.
H = frd(Hresp,freq);
H is an frd model object, which is a data container for representing frequency response
data.
You can use frd models with many frequency-domain analysis commands. For example,
visualize the response of this two-input, two-output system using bode.
2-37
2 Model Creation
Tip By default, the frd command assumes that the frequencies are in radians/second. To
specify different frequency units, use the TimeUnit and FrequencyUnit properties of
the frd model object. For example:
H = frd(Hresp,freq,'TimeUnit','min','FrequencyUnit','rad/TimeUnit')
See Also
frd
Related Examples
• “Frequency Response Data (FRD) Models” on page 2-10
More About
• “What Are Model Objects?” on page 1-2
• “Store and Retrieve Model Data” on page 3-2
2-38
Select Input/Output Pairs in MIMO Models
Note For more information about using cell arrays to create MIMO transfer
functions, see the tf reference page.
2 Select the response from the second input to the output of H.
H12 = H(1,2)
For any MIMO system H, the index notation H(i,j) selects the response from the jth
input to the ith output.
See Also
Related Examples
• “MIMO Transfer Functions” on page 2-28
• “MIMO State-Space Models” on page 2-31
More About
• “Store and Retrieve Model Data” on page 3-2
2-39
2 Model Creation
In discrete-time models, these properties are constrained to integer values that represent
delays expressed as integer multiples of the sample time. To approximate discrete-time
models with delays that are a fractional multiple of the sample time, use thiran.
To create the following first-order transfer function with a 2.1 s time delay:
1
G ( s ) = e-2.1s ,
s + 10
enter:
G = tf(1,[1 10],'InputDelay',2.1)
where InputDelay specifies the delay at the input of the transfer function.
Tip You can use InputDelay with zpk the same way as with tf:
G = zpk([],-10,1,'InputDelay',2.1)
For SISO transfer functions, a delay at the input is equivalent to a delay at the output.
Therefore, the following command creates the same transfer function:
G = tf(1,[1 10],'OutputDelay',2.1)
Use dot notation to examine or change the value of a time delay. For example, change the
time delay to 3.2 as follows:
2-40
Time Delays in Linear Systems
G.OutputDelay = 3.2;
ans =
3.2000
Tip An alternative way to create a model with a time delay is to specify the transfer
function with the delay as an expression in s:
1 Create a transfer function model for the variable s.
s = tf('s');
2 Specify G(s) as an expression in s.
G = exp(-2.1*s)/(s+10);
dx ( t )
= -2 x ( t ) + 3u ( t - 1 .5 )
dt
È x ( t - 0 .7 ) ˘
y (t ) = Í ˙.
Î - x (t ) ˚
This system has an input delay of 1.5. The first output has an output delay of 0.7, and the
second output is not delayed.
Note In contrast to SISO transfer functions, input delays are not equivalent to output
delays for state-space models. Shifting a delay from input to output in a state-space model
requires introducing a time shift in the model states. For example, in the model of this
2-41
2 Model Creation
example, defining T = t – 1.5 and X(T) = x(T + 1.5) results in the following equivalent
system:
dX ( T )
= -2 X ( T ) + 3u ( T )
dT
È X ( T - 2 .2 ) ˘
y (T ) = Í ˙.
Î - X ( T - 1 .5 ) ˚
All of the time delays are on the outputs, but the new state variable X is time-shifted
relative to the original state variable x. Therefore, if your states have physical meaning,
or if you have known state initial conditions, consider carefully before shifting time delays
between inputs and outputs.
A = -2;
B = 3;
C = [1;-1];
D = 0;
2 Create the model.
G = ss(A,B,C,D,'InputDelay',1.5,'OutputDelay',[0.7;0])
G is a ss model.
Tip Use delayss to create state-space models with more general combinations of input,
output, and state delays, of the form:
N
dx
= Ax (t ) + Bu( t ) + Â ( Ajx (t - tj ) + Bju (t - t j))
dt j =1
N
y (t ) = Cx (t ) + Du (t ) + Â ( Cjx( t - tj) + Dju (t - tj ))
j =1
2-42
Time Delays in Linear Systems
È - 0 .1 2 s +1 ˘
Íe e - 0 .3
s s + 10 ˙
H (s ) = Í ˙.
Í 10 -0.2 s - 1 ˙
e
ÍÎ s + 5 ˙˚
Time delays in MIMO systems can be specific to each I/O pair, as in this example. You
cannot use InputDelay and OutputDelay to model I/O-specific transport delays.
Instead, use ioDelay to specify the transport delay across each I/O pair.
s = tf('s');
2 Use the variable s to specify the transfer functions of H without the time delays.
H is a two-input, two-output tf model. Each I/O pair in H has the time delay specified by
the corresponding entry in tau.
2-43
2 Model Creation
H = tf(2,[1 -0.95],0.1,'InputDelay',25)
H =
2
z^(-25) * --------
z - 0.95
If system has a time delay that is not an integer multiple of the sampling time, you can
use the thiran command to approximate the fractional portion of the time delay with an
all-pass filter. See “Time-Delay Approximation” on page 2-48.
See Also
Related Examples
• “Closing Feedback Loops with Time Delays” on page 2-45
• “Convert Time Delay in Discrete-Time Model to Factors of 1/z” on page 2-64
More About
• “Time-Delay Approximation” on page 2-48
2-44
Closing Feedback Loops with Time Delays
2.3 e-2.1s
r 0.5 + s +10
y
- s
C G
G is the plant model, which has an input delay. C is a proportional-integral (PI) controller.
G = tf(1,[1 10],'InputDelay',2.1);
C = pid(0.5,2.3);
T = feedback(C*G,1);
The time delay in T is not an input delay as it is in G. Because the time delay is internal to
the closed-loop system, the software returns T as an ss model with an internal time delay
of 2.1 seconds.
T.InternalDelay
2-45
2 Model Creation
step(T)
Note Most analysis commands, such as step, bode and margin, support models with
internal delays.
The internal time delay is stored in the InternalDelay property of T. Use dot notation
to access InternalDelay. For example, to change the internal delay to 3.5 seconds,
enter:
T.InternalDelay = 3.5
2-46
See Also
You cannot modify the number of internal delays because they are structural properties of
the model.
See Also
Related Examples
• “Convert Time Delay in Discrete-Time Model to Factors of 1/z” on page 2-64
More About
• “Internal Delays” on page 2-73
2-47
2 Model Creation
Time-Delay Approximation
Many control design algorithms cannot handle time delays directly. For example,
techniques such as root locus, LQG, and pole placement do not work properly if time
delays are present. A common technique is to replace delays with all-pass filters that
approximate the delays.
To approximate time delays in continuous-time LTI models, use the pade command to
compute a Padé approximation. The Padé approximation is valid only at low frequencies,
and provides better frequency-domain approximation than time-domain approximation. It
is therefore important to compare the true and approximate responses to choose the right
approximation order and check the approximation validity.
Use the thiran command to approximate a time delay that is a fractional multiple of the
sample time as a Thiran all-pass filter.
For a time delay of tau and a sample time of Ts, the syntax thiran(tau,Ts) creates a
discrete-time transfer function that is the product of two terms:
• A term representing the integer portion of the time delay as a pure line delay, (1/z)N,
where N = ceil(tau/Ts).
• A term approximating the fractional portion of the time delay (tau - NTs) as a Thiran
all-pass filter.
Discretizing a Padé approximation does not guarantee good phase matching between the
continuous-time delay and its discrete approximation. Using thiran to generate a
discrete-time approximation of a continuous-time delay can yield much better phase
matching. For example, the following figure shows the phase delay of a 10.2-second time
delay discretized with a sample time of 1 s, approximated in three ways:
2-48
See Also
The Thiran filter yields the closest approximation of the 10.2-second delay.
See the thiran reference page for more information about Thiran filters.
See Also
absorbDelay | pade | thiran
Related Examples
• “Time-Delay Approximation in Continuous-Time Open-Loop Model” on page 2-50
• “Convert Time Delay in Discrete-Time Model to Factors of 1/z” on page 2-64
• “Approximate Different Delays with Different Approximation Orders” on page 2-60
2-49
2 Model Creation
Padé approximation is helpful when using analysis or design tools that do not support
time delays.
1
u e-2.6s y
s 2 + 0.9s + 1
s = tf('s');
P = exp(-2.6*s)/(s^2+0.9*s+1);
Pnd1 = pade(P,1)
Pnd1 =
-s + 0.7692
----------------------------------
s^3 + 1.669 s^2 + 1.692 s + 0.7692
h = bodeoptions;
h.PhaseMatching = 'on';
bodeplot(P,'-b',Pnd1,'-.r',{0.1,10},h)
legend('Exact delay','First-Order Pade','Location','SouthWest')
2-50
Time-Delay Approximation in Continuous-Time Open-Loop Model
The magnitude of P and Pnd1 match exactly. However, the phase of Pnd1 deviates
from the phase of P beyond approximately 1 rad/s.
4 Increase the Padé approximation order to extend the frequency band in which the
phase approximation is good.
Pnd3 = pade(P,3);
5 Compare the frequency response of P, Pnd1 and Pnd3.
bodeplot(P,'-b',Pnd3,'-.r',Pnd1,':k',{0.1 10},h)
legend('Exact delay','Third-Order Pade','First-Order Pade',...
'Location','SouthWest')
2-51
2 Model Creation
stepplot(P,'-b',Pnd3,'-.r',Pnd1,':k')
legend('Exact delay','Third-Order Pade','First-Order Pade',...
'Location','Southeast')
2-52
Time-Delay Approximation in Continuous-Time Open-Loop Model
Using the Padé approximation introduces a nonminimum phase artifact (“wrong way”
effect) in the initial transient response. The effect is quite pronounced in the first-
order approximation, which dips significantly below zero before changing direction.
The effect is reduced in the higher-order approximation, which far more closely
matches the exact system’s response.
Note Using too high an approximation order may result in numerical issues and
possibly unstable poles. Therefore, avoid Padé approximations with order N>10.
2-53
2 Model Creation
See Also
pade
Related Examples
• “Time-Delay Approximation in Continuous-Time Closed-Loop Model” on page 2-55
More About
• “Time-Delay Approximation” on page 2-48
• “Internal Delays” on page 2-73
2-54
Time-Delay Approximation in Continuous-Time Closed-Loop Model
Padé approximation is helpful when using analysis or design tools that do not support
time delays.
0.15 e-4.2s s +1
r 0.06 + + 0.006s s2 + 0.68s + 1 y
- s
C G
s = tf('s');
G = (s+1)/(s^2+.68*s+1)*exp(-4.2*s);
C = pid(0.06,0.15,0.006);
Tcl = feedback(G*C,1);
Tcl.InternalDelay
ans = 4.2000
2 Compute the first-order Padé approximation of Tcl.
Tnd1 = pade(Tcl,1);
h = bodeoptions;
h.PhaseMatching = 'on';
bodeplot(Tcl,'-b',Tnd1,'-.r',{.1,10},h);
legend('Exact delay','First-Order Pade','Location','SouthWest');
2-55
2 Model Creation
The magnitude and phase approximation errors are significant beyond 1 rad/s.
4 Compare the time domain response of Tcl and Tnd1 using stepplot.
stepplot(Tcl,'-b',Tnd1,'-.r');
legend('Exact delay','First-Order Pade','Location','SouthEast');
2-56
Time-Delay Approximation in Continuous-Time Closed-Loop Model
Using the Padé approximation introduces a nonminimum phase artifact (“wrong way”
effect) in the initial transient response.
5 Increase the Padé approximation order to see if this will extend the frequency with
good phase and magnitude approximation.
Tnd3 = pade(Tcl,3);
6 Observe the behavior of the third-order Padé approximation of Tcl. Compare the
frequency response of Tcl and Tnd3.
bodeplot(Tcl,'-b',Tnd3,'-.r',Tnd1,'--k',{.1,10},h);
legend('Exact delay','Third-Order Pade','First-Order Pade',...
'Location','SouthWest');
2-57
2 Model Creation
The magnitude and phase approximation errors are reduced when a third-order Padé
approximation is used.
Increasing the Padé approximation order extends the frequency band where the
approximation is good. However, too high an approximation order may result in numerical
issues and possibly unstable poles. Therefore, avoid Padé approximations with order
N>10.
See Also
pade
2-58
See Also
Related Examples
• “Approximate Different Delays with Different Approximation Orders” on page 2-60
More About
• “Time-Delay Approximation” on page 2-48
• “Internal Delays” on page 2-73
2-59
2 Model Creation
Load a sample continuous-time open-loop system that contains internal and output time
delays.
load(fullfile(matlabroot,'examples','control','PadeApproximation1.mat'),'sys')
sys
sys =
A =
x1 x2
x1 -1.5 -0.1
x2 1 0
B =
u1
x1 1
x2 0
C =
x1 x2
y1 0.5 0.1
D =
u1
y1 0
sys is a second-order continuous-time ss model with internal delay 3.4 s and output
delay 1.5 s.
2-60
Approximate Different Delays with Different Approximation Orders
Use the pade function to compute a third-order approximation of the internal delay and a
first-order approximation of the output delay.
P13 = pade(sys,inf,1,3);
size(P13)
The three input arguments following sys specify the approximation orders of any input,
output, and internal delays of sys, respectively. inf specifies that a delay is not to be
approximated. The approximation orders for the output and internal delays are one and
three respectively.
Approximating the time delays with pade absorbs delays into the dynamics, adding as
many states to the model as orders in the approximation. Thus, P13 is a sixth-order model
with no delays.
For comparison, approximate only the internal delay of sys, leaving the output delay
intact.
P3 = pade(sys,inf,inf,3);
size(P3)
P3.OutputDelay
ans = 1.5000
P3.InternalDelay
ans =
P3 retains the output delay, but the internal delay is approximated and absorbed into the
state-space matrices, resulting in a fifth-order model without internal delays.
Compare the frequency response of the exact and approximated systems sys, P13, P3.
h = bodeoptions;
h.PhaseMatching = 'on';
bode(sys,'b-',P13,'r-.',P3,'k--',h,{.01,10});
legend('sys','approximated output and internal delays','approximated internal delay onl
'location','SouthWest')
2-61
2 Model Creation
Notice that approximating the internal delay loses the gain ripple displayed in the exact
system.
See Also
pade
Related Examples
• “Time-Delay Approximation in Continuous-Time Open-Loop Model” on page 2-50
2-62
See Also
More About
• “Time-Delay Approximation” on page 2-48
• “Internal Delays” on page 2-73
2-63
2 Model Creation
Closing the feedback loop on a plant with input delays gives rise to internal delays in the
closed-loop system. Examine the order and internal delay of T.
order(T)
ans =
T.InternalDelay
ans =
2-64
Convert Time Delay in Discrete-Time Model to Factors of 1/z
Tnd = absorbDelay(T);
This command converts the internal delay to seven poles at z = 0. To confirm this,
examine the order and internal delay of Tnd.
order(Tnd)
ans =
Tnd.InternalDelay
ans =
Tnd has no internal delay, but it is a ninth-order model, due to the seven extra poles
introduced by absorbing the seven-unit delay into the model dynamics.
Despite this difference in representation, the responses of Tnd exactly match those of T.
stepplot(T,Tnd,'r--')
legend('T','Tnd')
2-65
2 Model Creation
bodeplot(T,Tnd,'r--')
legend('T','Tnd')
2-66
See Also
See Also
pade
Related Examples
• “Time-Delay Approximation in Continuous-Time Open-Loop Model” on page 2-50
More About
• “Time-Delay Approximation” on page 2-48
2-67
2 Model Creation
2-68
Frequency Response Data (FRD) Model with Time Delay
When you collect frequency response data for a system that includes time delays, you can
absorb the time delay into the frequency response as a phase shift. Alternatively, if you
are able to separate time delays from your measured frequency response, you can
represent the delays using the InputDelay, OutputDelay, or ioDelay properties of the
frd model object. The latter approach can give better numerical results, as this example
illustrates.
The frd model fsys includes a transport delay of 2 s. Load the model into the MATLAB®
workspace and inspect the time delay.
load(fullfile(matlabroot,'examples','control','frddelayexample.mat'),'fsys')
fsys.IODelay
ans =
A Bode plot of fsys shows the effect of the transport delay, causing the accumulation of
phase as frequency increases.
bodeplot(fsys)
2-69
2 Model Creation
The absorbDelay command absorbs all time delays directly into the frequency response,
resulting in an frd model with IODelay = 0.
fsys2 = absorbDelay(fsys);
fsys2.IODelay
ans =
Comparing the two ways of representing the delay shows that absorbing the delay into
the frequency response causes phase-wrapping.
2-70
See Also
bode(fsys,fsys2)
Phase wrapping can introduce numerical inaccuracy at high frequencies or where the
frequency grid is sparse. For that reason, if your system takes the form , you
might get better results by measuring frequency response data for G(s) and using
InputDelay, OutputDelay, or ioDelay to model the time delay .
See Also
absorbDelay
2-71
2 Model Creation
More About
• “Time-Delay Approximation” on page 2-48
2-72
Internal Delays
Internal Delays
Using the InputDelay, OutputDelay, and ioDelay properties, you can model simple
processes with transport delays. However, these properties cannot model more complex
situations, such as feedback loops with delays. In addition to the InputDelay and
OutputDelay properties, state-space (ss) models have an InternalDelay property.
This property lets you model the interconnection of systems with input, output, or
transport delays, including feedback loops with delays. You can use InternalDelay
property to accurately model and analyze arbitrary linear systems with delays. Internal
delays can arise from the following:
2-73
2 Model Creation
e-2s
s+2
e- 2 s
s + 2 + e- 2s
The delay term in the numerator can be represented as an output delay. However, the
delay term in the denominator cannot. In order to model the effect of the delay on the
feedback loop, the InternalDelay property is needed to keep track of internal coupling
between delays and ordinary dynamics.
Typically, you do not create state-space models with internal delays directly, by specifying
the A, B, C, and D matrices together with a set of internal delays. Rather, such models
arise when you interconnect models having delays. There is no limitation on how many
delays are involved and how the models are connected. For an example of creating an
internal delay by closing a feedback loop, see “Closing Feedback Loops with Time Delays”
on page 2-45.
• When a model interconnection gives rise to internal delays, the software returns an ss
model regardless of the interconnected model types. This occurs because only ss
supports internal delays.
• The software fully supports feedback loops. You can wrap a feedback loop around any
system with delays.
• When displaying the A, B, C, and D matrices, the software sets all delays to zero
(creating a zero-order Padé approximation). This approximation occurs for the display
only, and not for calculations using the model.
For some systems, setting delays to zero creates singular algebraic loops, which result
in either improper or ill-defined, zero-delay approximations. For these systems:
2-74
Internal Delays
• Entering sys returns only sizes for the matrices of a system named sys.
• Entering sys.A produces an error.
The limited display and the error do not imply a problem with the model sys itself.
2-75
2 Model Creation
You need not bother with this internal representation to use the tools. If, however, you
want to extract H or the matrices A, B1, B2, ... , you can use getDelayModel, For the
example:
P = 5*exp(-3.4*s)/(s+1);
C = 0.1 * (1 + 1/(5*s));
T = feedback(ss(P*C),1);
[H,tau] = getDelayModel(T,'lft');
size(H)
Note that H is a two-input, two-output model whereas T is SISO. The inverse operation
(combining H and tau to construct T) is performed by setDelayModel.
The following commands support internal delays for both continuous- and discrete-time
systems and have certain limitations:
2-76
See Also
To use these functions on a system with internal delays, use pade to approximate the
internal delays. See “Time-Delay Approximation” on page 2-48.
References
[1] P. Gahinet and L.F. Shampine, "Software for Modeling and Analysis of Linear Systems
with Delays," Proc. American Control Conf., Boston, 2004, pp. 5600-5605
See Also
Related Examples
• “Closing Feedback Loops with Time Delays” on page 2-45
2-77
2 Model Creation
You cannot use tunableTF to represent F, because the numerator and denominator
coefficients of a tunableTF block are independent. Instead, construct F using the
tunable real parameter object realp.
a = realp('a',10);
F = tf(a,[1 a]);
F is a genss object which has the tunable parameter a in its Blocks property. You can
connect F with other tunable or numeric models to create more complex control system
models. For example, see “Control System with Tunable Components” on page 2-83.
See Also
genss | realp | tunableTF
More About
• “Models with Tunable Coefficients” on page 1-19
• “Create Tunable Second-Order Filter” on page 2-79
• “Create State-Space Model with Both Fixed and Tunable Parameters” on page 2-81
• “Control System with Tunable Components” on page 2-83
2-78
Create Tunable Second-Order Filter
where the damping and the natural frequency are tunable parameters.
wn and zeta are realp parameter objects, with initial values 3 and 0.8, respectively.
The inputs to tf are the vectors of numerator and denominator coefficients expressed in
terms of wn and zeta.
F is a genss model. The property F.Blocks lists the two tunable parameters wn and
zeta.
F.Blocks
You can examine the number of tunable blocks in a generalized model using nblocks.
nblocks(F)
ans = 6
F has two tunable parameters, but the parameter wn appears five times - Twice in the
numerator and three times in the denominator.
2-79
2 Model Creation
nblocks(F)
ans = 4
In the new formulation, there are only three occurrences of the tunable parameter wn.
Reducing the number of occurrences of a block in a model can improve the performance
of calculations involving the model. However, the number of occurrences does not affect
the results of tuning the model or sampling it for parameter studies.
See Also
genss | nblocks | realp
More About
• “Models with Tunable Coefficients” on page 1-19
• “Create Tunable Low-Pass Filter” on page 2-78
• “Create State-Space Model with Both Fixed and Tunable Parameters” on page 2-81
• “Control System with Tunable Components” on page 2-83
2-80
Create State-Space Model with Both Fixed and Tunable Parameters
where a and b are tunable parameters, whose initial values are -1 and 3, respectively.
A is a generalized matrix whose Blocks property contains a and b. The initial value of A
is [1 2;0 -3], from the initial values of a and b.
sys =
Type "ss(sys)" to see the current value, "get(sys)" to see all properties, and "sys.Blo
sys is a generalized LTI model (genss) with tunable parameters a and b. Confirm that
the A property of sys is stored as a generalized matrix.
2-81
2 Model Creation
sys.A
ans =
Type "double(ans)" to see the current value, "get(ans)" to see all properties, and "ans
See Also
More About
• “Models with Tunable Coefficients” on page 1-19
• “Create Tunable Low-Pass Filter” on page 2-78
• “Create Tunable Second-Order Filter” on page 2-79
• “Control System with Tunable Components” on page 2-83
2-82
Control System with Tunable Components
Create models representing the plant and sensor dynamics. Since the plant and sensor
dynamics are fixed, represent them using numeric LTI models zpk and tf.
G = zpk([],[-1,-1],1);
S = tf(5,[1 4]);
a is a realp (real tunable parameter) object with initial value 10. Using a as a coefficient
in tf creates the tunable genss model object F.
Connect the models together to construct a model of the closed-loop response from to
.
2-83
2 Model Creation
T = feedback(G*C,S)*F
T =
Type "ss(T)" to see the current value, "get(T)" to see all properties, and "T.Blocks" t
T.Blocks
You can use tuning commands such as systune to tune the free parameters of T to meet
design requirements you specify.
See Also
Related Examples
• “Create Tunable Low-Pass Filter” on page 2-78
• “Create Tunable Second-Order Filter” on page 2-79
• “Create State-Space Model with Both Fixed and Tunable Parameters” on page 2-81
More About
• “Models with Tunable Coefficients” on page 1-19
2-84
Control System with Multichannel Analysis Points
The plant G has two inputs and two outputs. Therefore, the line marked y in the block
diagram represents two signals, y(1) and y(2). Similarly, r and e each represent two
signals.
Suppose you want to create tuning requirements or extract responses that require
injecting or measuring signals at the locations L and V. To do so, create an
AnalysisPoint block and include it in the closed-loop model of the control system as
shown in the following illustration.
To create a model of this system, first create the numeric LTI models and control design
blocks that represent the plant and controller elements. D is a tunable gain block, and
C_L and C_V are tunable PI controllers. Suppose the plant model is the following:
s = tf('s');
G = [87.8 -86.4 ; 108.2 -109.6]/(75*s+1);
2-85
2 Model Creation
D = tunableGain('Decoupler',eye(2));
C_L = tunablePID('C_L','pi');
C_V = tunablePID('C_V','pi');
AP_1 = AnalysisPoint('AP_1',2)
AP_1 =
Type "ss(AP_1)" to see the current value and "get(AP_1)" to see all properties.
AP_1.Location = {'L';'V'}
AP_1 =
Type "ss(AP_1)" to see the current value and "get(AP_1)" to see all properties.
The following diagram illustrates the input names, output names, and channel names
(locations) in the block AP_1.
The input and output names of the AnalysisPoint block are distinct from the channel
names. Use the channel names to refer to the analysis-point locations when extracting
responses or defining design goals for tuning. You can use the input and output names
AP_1.u and AP_1.y, for example, when interconnecting blocks using the connect
command.
2-86
Control System with Multichannel Analysis Points
You can now build the closed-loop model of the control system. First, join all the plant and
controller blocks along with the first AnalysisPoint block.
GC = G*AP_1*append(C_L,C_V)*D;
Then, close the feedback loop. Recall that GC has two inputs and outputs.
CL = feedback(GC,eye(2));
You can now use the analysis points for analysis or tuning. For example, extract the SISO
closed-loop transfer function from 'L' to the first output. Assign a name to the output so
you can reference it in analysis functions. The software automatically expands the
assigned name 'y' to the vector-valued output signals {y(1),y(2)}.
CL.OutputName = 'y';
TLy1 = getIOTransfer(CL,'L','y(1)');
bodeplot(TLy1);
2-87
2 Model Creation
See Also
AnalysisPoint
More About
• “Mark Signals of Interest for Control System Analysis and Design” on page 2-89
2-88
Mark Signals of Interest for Control System Analysis and Design
Analysis Points
Whether you model your control system in MATLAB or Simulink, use analysis points to
mark points of interest in the model. Analysis points allow you to access internal signals,
perform open-loop analysis, or specify requirements for controller tuning. In the block
diagram representation, an analysis point can be thought of as an access port to a signal
flowing from one block to another. In Simulink, analysis points are attached to the
outports of Simulink blocks. For example, in the following model, the reference signal, r,
and the control signal, u, are analysis points that originate from the outputs of the
setpoint and C blocks respectively.
Each analysis point can serve one or more of the following purposes:
• Input — The software injects an additive input signal at an analysis point, for
example, to model a disturbance at the plant input.
• Output — The software measures the signal value at a point, for example, to study the
impact of a disturbance on the plant output.
• Loop Opening — The software inserts a break in the signal flow at a point, for
example, to study the open-loop response at the plant input.
You can apply these purposes concurrently. For example, to compute the open-loop
response from u to y, you can treat u as both a loop opening and an input. When you use
an analysis point for more than one purpose, the software applies the purposes in this
sequence: output measurement, then loop opening, then input.
2-89
2 Model Creation
Using analysis points, you can extract open-loop and closed-loop responses from a control
system model. For example, suppose T represents the closed-loop system in the model
above, and u and y are marked as analysis points. T can be either a generalized state-
space model or an slLinearizer or slTuner interface to a Simulink model. You can
plot the closed-loop response to a step disturbance at the plant input with the following
commands:
Tuy = getIOTransfer(T,'u','y');
stepplot(Tuy)
Analysis points are also useful to specify design requirements when tuning control
systems with the systune command. For example, you can create a requirement that
attenuates disturbances at the plant input by a factor of 10 (20 dB) or more.
Req = TuningGoal.Rejection('u',10);
2-90
Mark Signals of Interest for Control System Analysis and Design
G = tf(10,[1 3 10]);
C = pid(0.2,1.5);
T = feedback(G*C,1);
With this model, you can obtain the closed-loop response from r to y. However, you
cannot analyze the open-loop response at the plant input or simulate the rejection of a
step disturbance at the plant input. To enable such analysis, mark the signal u as an
analysis point by inserting an AnalysisPoint block between the plant and controller.
AP = AnalysisPoint('u');
T = feedback(G*AP*C,1);
T.OutputName = 'y';
In creating the model T, you manually created the analysis point block AP and explicitly
included it in the feedback loop. When you combine models using the connect command,
you can instruct the software to insert analysis points automatically at the locations you
specify. For more information, see connect.
To mark an analysis point explicitly in the model, right-click a signal and, under Linear
Analysis Points, select an analysis point type.
2-91
2 Model Creation
You can select any of the following closed-loop analysis point types, which are equivalent
within an slLinearizer or slTuner interface; that is, they are treated the same way by
analysis functions, such as getIOTransfer, and tuning goals, such as
TuningGoal.StepTracking.
2-92
Mark Signals of Interest for Control System Analysis and Design
• Input Perturbation
• Output Measurement
• Sensitivity
• Complementary Sensitivity
If you want to introduce a permanent loop opening at a signal as well, select one of the
following open-loop analysis point types:
• Open-Loop Input
• Open-Loop Output
• Loop Transfer
• Loop Break
When you create an slLinearizer or slTuner interface for a model, any analysis
points defined in the model are automatically added to the interface. If you defined an
analysis point using:
To mark analysis points programmatically, use the addPoint command. For example,
consider the scdcascade model.
open_system('scdcascade')
2-93
2 Model Creation
ST = slTuner('scdcascade');
To add a signal as an analysis point, use the addPoint command, specifying the source
block and port number for the signal.
addPoint(ST,'scdcascade/C1',1);
If the source block has a single output port, you can omit the port number.
addPoint(ST,'scdcascade/G2');
For convenience, you can also mark analysis points using the:
addPoint(ST,'y2');
• Combined source block path and port number.
addPoint(ST,'scdcascade/C1/1')
• End of the full source block path when unambiguous.
addPoint(ST,'G1/1')
You can also add permanent openings to an slLinearizer or slTuner interface using
the addOpening command, and specifying signals in the same way as for addPoint. For
2-94
Mark Signals of Interest for Control System Analysis and Design
more information on how the software treats loop openings during linearization, see
“How the Software Treats Loop Openings” (Simulink Control Design).
addOpening(ST,'y1m');
You can also define analysis points by creating linearization I/O objects using the linio
command.
io(1) = linio('scdcascade/C1',1,'input');
io(2) = linio('scdcascade/G1',1,'output');
addPoint(ST,io);
As when you define analysis points directly in your model, if you specify a linearization I/O
object with:
When you specify response I/Os in a tool such as Linear Analysis Tool or Control System
Tuner, the software creates analysis points as needed.
You can also create tuning goals that constrain the system response at these points. The
tools to perform these operations operate in a similar manner for models created at the
command line and models created in Simulink.
To view the available analysis points, use the getPoints function. You can view the
analysis for models created:
2-95
2 Model Creation
For closed-loop models created at the command line, you can also use the model input
and output names when:
2-96
Mark Signals of Interest for Control System Analysis and Design
Use the same method to refer to analysis points for models created in Simulink. In
Simulink models, for convenience, you can use any unambiguous abbreviation of the
analysis point names returned by getPoints.
ioSys = getIOTransfer(ST,'u1','y1');
sensG2 = getSensitivity(ST,'G2');
R = TuningGoal.Margins('u1',10,60);
Finally, if some analysis points are vector-valued signals or multichannel locations, you
can use indices to select particular entries or channels. For example, suppose u is a two-
entry vector in a closed-loop MIMO model.
2-97
2 Model Creation
You can compute the open-loop response of the second channel and measure the impact
of a disturbance on the first channel.
L = getLoopTransfer(T,'u(2)',-1);
stepplot(getIOTransfer(T,'u(1)','y'))
When you create tuning goals in Control System Tuner, the software creates analysis
points as needed.
2-98
See Also
See Also
AnalysisPoint | getIOTransfer | getPoints
More About
• “Control System with Multichannel Analysis Points” on page 2-85
• “Mark Analysis Points in Closed-Loop Models” on page 4-13
2-99
2 Model Creation
Model Arrays
and so on. Model arrays are a convenient way to store and analyze such a collection.
Model arrays are collections of multiple linear models, stored as elements in a single
MATLAB array.
For all models collected in a single model array, the following attributes must be the
same:
Using model arrays, you can apply almost all of the basic model operations that work on
single model objects to entire sets of models at once. Functions operate on arrays model
by model, allowing you to manipulate an entire collection of models in a vectorized
fashion. You can also use analysis functions such as bode, nyquist, and step to model
2-100
Model Arrays
arrays to analyze multiple models simultaneously. You can access the individual models in
the collection through MATLAB array indexing.
Just as you might collect a set of two-by-two matrices in a multidimensional array, you can
collect this set of five transfer function models as a list in a model array under one
variable name, say, sys. Each element of the model array is a single model object.
2-101
2 Model Creation
The following illustration shows selection of models from the two-dimensional model
array m2d.
See Also
Related Examples
• “Query Array Size and Characteristics” on page 2-106
• “Select Models from Array” on page 2-103
• “Model Array with Variations in Two Parameters” on page 9-6
2-102
Select Models from Array
1 Load the transfer function array m2d into the MATLAB workspace.
step(m2d)
The step response shows that m2d contains six one-input, two-output models. The
step command plots all of the models in an array on a single plot.
3 (Optional) Examine the dimensions of m2d.
arraydim = size(m2d)
arraydim =
2 1 2 3
2-103
2 Model Creation
• The first entries of arraydim, 2 and 1, show that m2d is an array of two-output,
one-input transfer functions.
• The remaining entries in arraydim give the array dimensions of m2d, 2-by-3.
sys = m2d(:,:,2,1)
Tip You can also access models using single index referencing of the array
dimensions. For example,
sys = m2d(:,:,4)
m11 = m2d(1,1,:,:)
6 (Optional) Plot the step response of m11.
step(m11)
2-104
See Also
The step response shows that m11 is an array of six single-input, single-output (SISO)
models.
Note For frequency response data (FRD) models, the array indices can be followed
by the keyword 'frequency' and some expression selecting a subset of the
frequency points, as in:
sys(outputs,inputs,n1,...,nk,'frequency',SelectedFreqs)
See Also
More About
• “Model Arrays” on page 2-100
• “Query Array Size and Characteristics” on page 2-106
2-105
2 Model Creation
Array Size
Model arrays have two different sets of dimensions, the I/O dimensions and the array
dimensions. The I/O dimensions are the numbers of inputs and outputs of the models in
the array. (Each model in an array must have the same I/O dimensions.) The array
dimensions are the dimensions of the array itself. Load a saved model array and query its
dimensions.
load(fullfile(matlabroot,'examples','control','queryexample.mat'),'sysarr')
size(sysarr)
When you use the size command on a model array with no output argument, the display
shows the two sets of dimensions.
To obtain the array dimensions as a numeric array, use size with an output argument.
dims = size(sysarr)
dims = 1×4
3 1 2 4
The first two entries in dims are the I/O dimensions of the models in sysarr, which each
have three outputs and one input. The remaining entries in dims are the dimensions of
the array itself. Thus, sysarr is a 2-by-4 array of models.
To query the number of dimensions in the array, rather than the values of those
dimensions, use ndims.
dims = ndims(sysarr)
dims = 4
2-106
Query Array Size and Characteristics
In this case, sysarr has 4 = 2 + 2 dimensions: The I/O dimensions (outputs and inputs),
and the array dimensions. Query the I/O dimensions alone using the iosize command.
ios = iosize(sysarr)
ios = 1×2
3 1
N = nmodels(sysarr)
N = 8
Query commands such as isproper and isstable work on model arrays. For example,
query whether the models in sysarr are stable.
Bsiso = isstable(sysarr)
Bsiso = logical
1
By default, isstable returns 1 (true) if all of the models in the array are stable. The
commands returns 0 (false) if one or more of the models is not stable. To perform an
element-by-element query of a model array, use the 'elem' option.
Bsiso = isstable(sysarr,'elem')
1 1 1 1
1 1 1 1
Now isstable returns an array of Boolean values. The dimensions of this array match
the array dimensions of sysarr. Each entry in the array Bsiso indicates whether the
2-107
2 Model Creation
corresponding model of sysarr is stable. The 'elem' option works similarly for many
query commands.
See Also
More About
• “Model Arrays” on page 2-100
• “Select Models from Array” on page 2-103
2-108
Linear Parameter-Varying Models
dx ( t ) = A ( p) x ( t ) + B ( p) u ( t )
y ( t ) = C ( p) x ( t ) + D ( p) u ( t)
x ( 0 ) = x0
where
2-109
2 Model Creation
For example, the aerodynamic behavior of an aircraft is often scheduled over a grid of
incidence angle (α) and wind speed (V) values. For each scheduling parameter, a range of
values is chosen, such as α = 0:5:20 degrees, V = 700:100:1400 m/s. For each
combination of (α,V) values, a linear approximation of the aircraft behavior is obtained.
The local models are connected as shown in the following figure:
Each donut represents a local LTI model, and the connecting curves represent the
interpolation rules. The abscissa and ordinate of the surface are the scheduling
parameters (α, V).
This form is sometimes called the grid-based LPV representation. This is the form used by
the LPV System block. For meaningful interpolations of system matrices, all the local
models must use the same state basis.
2-110
Linear Parameter-Varying Models
The LPV system representation can be extended to allow offsets in dx, x, u and y
variables. This form is known as affine form of the LPV model. Mathematically, the
following represents an LPV system:
(
dx ( t ) = A ( p) x ( t ) + B ( p) u ( t ) + dx ( p) - A ( p) x ( p) - B( p)u( p) )
y ( t ) = C ( p) x ( t ) + D ( p) u ( t) + ( y ( p) - C ( p) x ( p) - D( p)u ( p) )
x ( 0 ) = x0
dx ( p) , x ( p) , u ( p) , y ( p) are the offsets in the values of dx(t), x(t), u(t) and y(t) at
a given parameter value p = p(t).
To obtain such representations of the linear system array, linearize a Simulink model over
a batch of operating points (see “Batch Linearization” (Simulink Control Design).) The
offsets correspond to the operating points at which you linearized the model.
You can obtain the offsets by returning additional linearization information when calling
functions such as linearize or getIOTransfer. You can then extract the offsets using
getOffsetsForLPV. For an example, see “LPV Approximation of a Boost Converter
Model” (Simulink Control Design).
In the affine representation, the linear model at a given point p = p* in the scheduling
space is:
( ) ( ) ( ) ( )
dDx(t, p* ) = A p* Dx t, p* + B p* Du t, p*
D y( t, p ) = C ( p* ) Dx ( t, p* ) + D ( p* ) D u ( t, p* )
*
The states of this linear model are related to the states of the overall LPV model
2-111
2 Model Creation
When parameters co-vary, that is, α and β increase together, an irregular grid is formed.
The system array parameters are available only along the diagonal in the parameter
plane.
2-112
Linear Parameter-Varying Models
If certain samples are missing from an otherwise regular grid, the grid is considered to be
irregular.
2-113
2 Model Creation
The system array size is equal to the grid size in scheduling space. In the aircraft
example, α takes 5 values in the 0–20 degrees range and V takes 8 values in the 700–1400
m/s range. If you define a linear model at every combination of (α,V) values (i.e., the grid
is regular), the grid size is 5-by-8. Therefore, the model array size must be 5-by-8.
The information about scheduling parameters is attached to the linear model array using
its SamplingGrid property. The value of the SamplingGrid property must be a
structure with as many fields as there are scheduling parameters. For each field, the
value must be set to all the values assumed by the corresponding variable in the
scheduling space.
For the aircraft example, you can define the SamplingGrid property as:
Alpha = 0:5:20;
V = 700:100:1400;
2-114
Linear Parameter-Varying Models
Note When obtaining linear models by linearization, do not reduce or alter the state
variables used by the models.
The operating region is usually of a high dimension because it consists of all the input and
state variables. Generating or interpolating local models in such high-dimensional spaces
is usually infeasible. A simpler approach is to use a small set of scheduling parameters as
a proxy for the operating space variables. The scheduling parameters are derived from
the inputs and state variables of the original system. You must choose the values carefully
so that for a fixed value of the scheduling parameters, the system behavior is
approximately linear. This approach is not always possible.
Suppose you use p( t) = x& 1 as a scheduling variable. At a given time instant t = t0, you
have:
x& 1 ª 2 x1 ( t0 ) x1 + 2 x2 ( t0 ) x2 - x&1 ( t0 )
x& 2 = -2 x1 - 3 x2 + 2u
y = x1 + 2
2-115
2 Model Creation
Thus, the dynamics are linear (affine) in the neighborhood of a given value of x& . The
approximation holds for all time spans and values of input u as long as of x& does not
deviate much from its nominal value at sampling point t0. Note that scheduling on input u
or states x1 or x2 does not help locally linearize the system. Therefore, they are not good
candidates for scheduling parameters.
For an example of this approach, see “Approximating Nonlinear Behavior Using an Array
of LTI Systems” (Simulink Control Design).
You can use LPV models to represent systems that exhibit multiple modes (regimes) of
operation. Examples of such systems include colliding bodies, systems controlled by
operator switches, and approximations of systems affected by dry friction and hysteresis
effects. For an example, see “Using LTI Arrays for Simulating Multi-Mode Dynamics” on
page 2-118.
This approach is useful for generating surrogate models that you can use in place of the
original system for enabling faster simulations, reducing memory footprint of target
hardware code, and hardware-in-loop (HIL) simulations. You can also use surrogate
models of this type for designing gain-scheduled controllers and for initializing the
parameter estimation tasks in Simulink. For an example of approximating a general
nonlinear system behavior by an LPV model, see “Approximating Nonlinear Behavior
Using an Array of LTI Systems” (Simulink Control Design).
LPV models can help speed up the simulation of physical component based systems, such
as those built using Simscape™ Multibody™ and Simscape Power Systems™ software.
For an example of this approach, see “LPV Approximation of a Boost Converter Model”
(Simulink Control Design).
See Also
LPV System | getOffsetsForLPV
2-116
See Also
More About
• “Using LTI Arrays for Simulating Multi-Mode Dynamics” on page 2-118
• “Approximating Nonlinear Behavior Using an Array of LTI Systems” (Simulink
Control Design)
• “LPV Approximation of a Boost Converter Model” (Simulink Control Design)
2-117
2 Model Creation
Introduction
We often encounter situations where an elastic body collides with, or presses against, a
possibly elastic surface. Examples of such situations are:
In these situations, the motion of the moving body exhibits different dynamics when it is
moving freely than when it is in contact with a surface. In the case of a bouncing ball, the
motion of the mass can be described by rigid body dynamics when it is falling freely.
When the ball collides and deforms while in contact with the surface, the dynamics have
to take into account the elastic properties of the ball and of the surface. A simple way of
modeling the impact dynamics is to use lumped mass spring-damper descriptions of the
colliding bodies. By adjusting the relative stiffness and damping coefficients of the two
bodies, we can model the various situations described above.
Figure 1 shows a mass-spring-damper model of the system. Mass 1 is falling freely under
the influence of gravity. Its elastic properties are described by stiffness constant and
damping coefficient . When this mass hits the fixed surface, the impact causes Mass 1
and Mass 2 to move downwards together. After a certain "residence time" during which
the Mass 1 deforms and recovers, it loses contact with Mass 2 completely to follow a
projectile motion. The overall dynamics are thus broken into two distinct modes - when
the masses are not in contact and when they are moving jointly.
2-118
Using LTI Arrays for Simulating Multi-Mode Dynamics
The unstretched (load-free) length of spring attached to Mass 1 is , while that of Mass 2
is . The variables and denote the positions of the two masses. When the
masses are not in contact ("Mode 1"), their motions are governed by the following
equations:
When Mass 1 touches Mass 2 ("Mode 2"), their displacements and velocities get
interlinked. The governing equations in this mode are:
LPV Representation
The governing equations are linear and time invariant. However, there are two distinct
behavioral modes corresponding to different equations of motion. Both modes are
governed by sets of second order equations. If we pick the positions and velocities of the
masses as state variables, we can represent each mode by a 4th order state-space
equation.
In the state-space view, it becomes possible to treat the two modes as a single system
whose coefficients change as a function of a certain condition which determines which
mode is active. The condition is, of course, whether the two masses are moving freely or
jointly. Such a representation, where the coefficients of a linear system are parameterized
by an external but measurable parameter is called a Linear Parameter Varying (LPV)
model. A common representation of an LPV model is by means of an array of linear state-
space models and a set of scheduling parameters that dictate the rules for choosing the
correct model under a given condition. The array of linear models must all be defined
using the same state variables.
2-119
2 Model Creation
For our example, we need two state-space models, one for each mode of operation. We
also need to define a scheduling variable to switch between them. We begin by writing the
above equations of motion in state-space form.
First mode: state-space representation of dynamics when the masses are not in contact.
A11 = [0 1; 0 0];
B11 = [0; -g];
C11 = [1 0];
D11 = 0;
A1 = blkdiag(A11, A12);
B1 = [B11; B12];
C1 = blkdiag(C11, C12);
D1 = [D11; D12];
sys1 = ss(A1,B1,C1,D1);
Second mode: state-space representation of dynamics when the masses are in contact.
A2 = [ 0 1, 0, 0; ...
-k1/m1, -c1/m1, k1/m1, c1/m1;...
2-120
Using LTI Arrays for Simulating Multi-Mode Dynamics
0, 0, 0, 1; ...
k1/m2, c1/m2, -(k1+k2)/m2, -(c1+c2)/m2];
sys2 = ss(A2,B2,C2,D2);
Now we stack the two models sys1 and sys2 together to create a state-space array.
sys = stack(1,sys1,sys2);
Use the information on whether the masses are moving freely or jointly for scheduling.
Let us call this parameter "FreeMove" which takes the value of 1 when masses are
moving freely and 0 when they are in contact and moving jointly. The scheduling
parameter information is incorporated into the state-space array object (sys) by using its
"SamplingGrid" property:
Whether the masses are in contact or not is decided by the relative positions of the two
masses; when , the masses are not in contact.
The state-space array sys has the necessary information to represent an LPV model. We
can simulate this model in Simulink using the "LPV System" block from the Control
System Toolbox™'s block library.
open_system('LPVBouncingMass')
open_system('LPVBouncingMass/Bouncing Mass Model','mask')
2-121
2 Model Creation
The block called "Bouncing Mass Model" is an LPV System block. Its parameters are
specified as follows:
• For "State-space array" field, specify the state-space model array sys that was created
above.
• For "Initial state" field, specify the initial positions and velocities of the two masses.
Note that the state vector is: . Specify its value as [h1 0 h2 0]'.
• Under the "Scheduling" tab, set the "Interpolation method" to "Nearest". This choice
causes only one of the two models in the array to be active at any time. In our
example, the behavior modes are mutually exclusive.
• Under the "Outputs" tab, uncheck all the checkboxes for optional output ports. We will
be observing only the positions of the two masses.
The constant block outputs a unit value. This serves as the input to the model and is
supplied from the first input port of the LPV block. The block has only one output port
which outputs the positions of the two masses as a 2-by-1 vector.
The second input port of the LPV block is for specifying the scheduling signal. As
discussed before, this signal represents the scheduling parameter "FreeMove" and takes
discrete values 0 (masses in contact) or 1 (masses not in contact). The value of this
parameter is computed as a function of the block's output signal. This computation is
performed by the blocks with cyan background color. We take the difference between the
two outputs (after demuxing) and compare the result to the unstretched length of spring
attached to Mass 1. The resulting Boolean result is converted into a double signal which
serves as the scheduling parameter value.
2-122
Using LTI Arrays for Simulating Multi-Mode Dynamics
open_system('LPVBouncingMass/Scope')
sim('LPVBouncingMass')
The yellow curve shows the position of Mass 1 while the magenta curve shows the
position of Mass 2. At the start of simulation, Mass 1 undergoes free fall until it hits Mass
2. The collision causes the Mass 2 to be displaced but it recoils quickly and bounces Mass
1 back. The two masses are in contact for the time duration where . When the
masses settle down, their equilibrium values are determined by the static settling due to
gravity. For example, the absolute location of Mass 1 is
2-123
2 Model Creation
Conclusions
This example shows how a Linear Parameter Varying model can be constructed by using
an array of state-space models and suitable scheduling variables. The example describes
the case of mutually exclusive modes, although a similar approach can be used in cases
where the dynamics behavior at a given value of scheduling parameters is influenced by
several linear models.
The LPV System block facilitates the simulation of parameter varying systems. The block
also supports code generation for various hardware targets.
2-124
Working with Linear Models
125
3
Data Manipulation
Model Properties
Model properties are the data fields that store all data about a dynamic system model.
Data stored in model properties includes model dynamics, such as transfer-function
coefficients, state-space matrices, and time delays. Model properties also let you specify
other model attributes such as sample time, channel names, and state names.
For information about the properties associated with each model type, see the
corresponding reference page, such as tf, pid, or ss.
You can specify other values for model properties at model creation using the
Name,Value pair syntax of the model-creation command. In this syntax, you specify the
name of the property you want to set, followed by the value. You can set multiple property
values in one command. For example, assign a transport delay and input and output
names to a new transfer function model.
H = tf(1,[1 10],'IODelay',6.5,'InputName','torque','OutputName','velocity')
H =
Some property values are reflected in the model display, such as the input and output
names. You can use Name,Value pair syntax when creating any type of model.
3-2
Store and Retrieve Model Data
sys =
A =
x1 x2
x1 -1.5 -0.1
x2 1 0
B =
u1
x1 1
x2 0
C =
x1 x2
y1 0.5 0.1
D =
u1
y1 0
The display shows that sys is a state-space model, and includes some of the property
values of sys. To see all properties of sys, use the get command.
get(sys)
A: [2x2 double]
B: [2x1 double]
C: [0.5000 0.1000]
D: 0
E: []
Scaled: 0
3-3
3 Data Manipulation
Use dot notation to access the values of particular properties. For example, display the A
matrix of sys.
Amat = sys.A
Amat = 2×2
-1.5000 -0.1000
1.0000 0
Dot notation also lets you change the value of individual model properties.
sys.InputDelay = 4.2;
sys.InputName = 'thrust';
sys.OutputName = 'velocity';
When you must change multiple property values at the same time to preserve the validity
of the model, such as changing the dimensions of the state-space matrices, you can use
the set command. For example, create a 1-state state-space model, and then replace the
matrices with new values representing a 2-state model.
sys2 = rss(1);
Anew = [-2, 1; 0.5 0];
Bnew = [1; -1];
Cnew = [0, -0.4];
3-4
See Also
set(sys2,'A',Anew,'B',Bnew,'C',Cnew)
sys2
sys2 =
A =
x1 x2
x1 -2 1
x2 0.5 0
B =
u1
x1 1
x2 -1
C =
x1 x2
y1 0 -0.4
D =
u1
y1 0.3426
See Also
Related Examples
• “Attach Metadata to Models” on page 3-9
• “Extract Model Coefficients” on page 3-6
3-5
3 Data Manipulation
Command Result
tfdata Extract transfer function coefficients
zpkdata Extract zero and pole locations and system gain
ssdata Extract state-space matrices
dssdata Extract descriptor state-space matrices
frdata Extract frequency response data from frd model
piddata Extract parallel-form PID data
pidstddata Extract standard-form PID data
get Access all model property values
s = tf('s');
H = exp(-2.5*s)/(s+12);
3-6
Extract Model Coefficients
The variables num and den are numerical arrays. Without the 'v' flag, tfdata
returns cell arrays.
Note For SISO transfer function models, you can also extract coefficients using:
num = H.Numerator{1};
den = H.Denominator{1};
In a SISO tf model, you can express a time delay as an input delay, an output
delay, or a transport delay (I/O delay).
get(H)
3-7
3 Data Manipulation
1 Create a transfer function that represents a PID controller with a first-order filter on
the derivative term.
Czpk = zpk([-6.6,-0.7],[0,-2],0.2)
2 Obtain the PID gains and filter constant.
[Kp,Ki,Kd,Tf] = piddata(Czpk)
This command returns the proportional gain Kp, integral gain Ki, derivative gain Kd,
and derivative filter time constant Tf. Because piddata automatically computes the
PID controller parameters, you can extract the PID coefficients without creating a
pid model.
See Also
Related Examples
• “Attach Metadata to Models” on page 3-9
More About
• “Store and Retrieve Model Data” on page 3-2
3-8
Attach Metadata to Models
The TimeUnit property of the tf model object specifies units of the time variable, time
delays (for continuous-time models), and the sample time Ts (for discrete-time models).
The default time units is seconds.
4s + 2
Create a SISO transfer function model sys = with time units in milliseconds:
2
s + 3 s + 10
num = [4 2];
den = [1 3 10];
sys = tf(num,den,'TimeUnit','milliseconds');
You can specify the time units of any dynamic system on page 1-10 in a similar way.
The system time units appear on the time- and frequency-domain plots. For multiple
systems with different time units, the units of the first system are used if the time and
frequency units in the “Toolbox Preferences Editor” on page 19-2 are auto.
Note Changing the TimeUnit property changes the system behavior. If you want to use
different time units without modifying system behavior, use chgTimeUnit.
1 Create two transfer function models with time units of milliseconds and seconds,
respectively.
3-9
3 Data Manipulation
The FrequencyUnit property specifies units of the frequency vector in the Frequency
property of the frd model object. The default frequency units are rad/TimeUnit, where
TimeUnit is the time unit specified in the TimeUnit property.
You can independently specify the units in which you measure the frequency points and
sample time in the FrequencyUnit and TimeUnit properties, respectively. You can also
specify the frequency units of a genfrd in a similar way.
The frequency units appear on the frequency-domain plots. For multiple systems with
different frequency units, the units of the first system are used if the frequency units in
the “Toolbox Preferences Editor” on page 19-2 is auto.
Note Changing the FrequencyUnit property changes the system behavior. If you want
to use different frequency units without modifying system behavior, use chgFreqUnit.
Extracting subsystems is useful when, for example, you want to analyze a portion of a
complex system.
3-10
Attach Metadata to Models
G1 = tf(3,[1 10]);
G2 = tf([1 2],[1 0]);
G = [G1,G2];
Gsub = G(:,1);
This command uses MATLAB indexing to specify a subsystem as G(out,in), where out
specifies the output indices and in specifies the input indices.
Using channel names, you can use MATLAB indexing to extract all the dynamics relating
to a particular channel. By using this approach, you can avoid having to keep track of
channel order in a complex MIMO model.
G.InputName = {'temperature';'pressure'};
Because G has two inputs, use a cell array to specify the two channel names.
Extract the subsystem of G that contains all dynamics from the 'temperature' input to
all outputs.
Gt = G(:,'temperature');
Note When you extract a subsystem from a state-space (ss) model, the resulting state-
space model may not be minimal. Use sminreal to eliminate unnecessary states in the
subsystem.
Input and output groups are useful for keeping track of inputs and outputs in complex
MIMO models.
3-11
3 Data Manipulation
H = rss(3,4,3);
2 Group the inputs as follows:
H.InputGroup.controls = [1 2];
H.OutputGroup.temperature = [1 3];
H.OutputGroup.measurements = [1 3 4];
InputGroup and OutputGroup are structures. The name of each field in the
structure is the name of the input or output group. The value of each field is a vector
that identifies the channels in that group.
3 Extract the subsystem corresponding to the controls inputs and the temperature
outputs.
Hc = H('temperature','controls')
You can see the relationship between H and the subsystem Hc in this illustration.
Hc
1 1
controls 2 temperature
2 H
3
3 4
3-12
See Also
See Also
Related Examples
• “Store and Retrieve Model Data” on page 3-2
• “Extract Model Coefficients” on page 3-6
• “Query Model Characteristics” on page 3-14
3-13
3 Data Manipulation
Bstab = logical
1
The isstable command returns 1 (true) if all system poles are in the open left-half
plane (for continuous-time models) or inside the open unit disk (for discrete-time models).
Otherwise, isstable command returns 0 (false). Here, the result shows that the model
is stable.
Bdel = logical
1
The returned value, 1, indicates that T has a time delay. For a state-space model, time
delay can be stored as input delay, output delay, internal delay, or a combination. Use
get(T) to determine which properties of T hold the time delay, and use dot notation to
access the delay values. The hasInternalDelay command tells you whether there is
any internal delay.
Bprop = logical
1
3-14
See Also
The returned value indicates that the system has relative degree less than or equal to 0.
This is true of a SISO system when it can be represented as a transfer function in which
the degree of the numerator does not exceed the degree of the denominator.
N = order(T)
N = 5
For a state-space model, order returns the number of states, which is 5 in this case. For
a tf or zpk model, the order is the number of states required for a state-space realization
of the system.
Bdisc = isdt(T)
Bdisc = logical
1
The returned value indicates that T is a discrete-time model. Similarly, use isct to query
whether T is a continuous-time model.
load(fullfile(matlabroot,'examples','control','queryexample.mat'),'Tmimo')
ios = iosize(Tmimo)
ios = 1×2
7 4
In the resulting array, the number of outputs is first. Therefore, Tmimo has 4 inputs and 7
outputs.
See Also
isproper | isstable | size
3-15
3 Data Manipulation
Related Examples
• “Select Models from Array” on page 2-103
More About
• “Store and Retrieve Model Data” on page 3-2
3-16
Customize Model Display
You can use the same steps to configure the display variable of transfer function models
in factorized form (zpk models).
By default, tf and zpk models are displayed in terms of s in continuous time and z in
discrete time. Use the Variable property change the display variable to 'p' (equivalent
to 's'), 'q' (equivalent to 'z'), 'z^-1', or 'q^-1'.
1
z -1
Create the discrete-time transfer function H ( z) =
2
z - 3z + 2
with a sample time of 1 s.
H =
z - 1
-------------
z^2 - 3 z + 2
H.Variable = 'q^-1'
H =
q^-1 - q^-2
-------------------
1 - 3 q^-1 + 2 q^-2
3-17
3 Data Manipulation
When you change the Variable property, the software computes new coefficients
and displays the transfer function in terms of the new variable. The num and den
properties are automatically updated with the new coefficients.
Tip Alternatively, you can directly create the same transfer function expressed in terms
of 'q^-1'.
For the inverse variables 'z^-1' and 'q^-1', tf interprets the numerator and
denominator arrays as coefficients of ascending powers of 'z^-1' or 'q^-1'.
You can configure the display of the factorized numerator and denominator polynomials to
highlight:
See the DisplayFormat property on the zpk reference page for more information about
these quantities.
1 Create a zpk model having a zero at s = 5, a pole at s = –10, and a pair of complex
poles at s = –3 ± 5i.
H = zpk(5,[-10,-3-5*i,-3+5*i],10)
H =
10 (s-5)
----------------------
3-18
Customize Model Display
The default display format, 'roots', displays the standard factorization of the
numerator and denominator polynomials.
2 Configure the display format to display the natural frequency of each polynomial root.
H.DisplayFormat = 'frequency'
H =
-0.14706 (1-s/5)
-------------------------------------------
(1+s/10) (1 + 1.029(s/5.831) + (s/5.831)^2)
You can read the natural frequencies and damping ratios for each pole and zero from
the display as follows:
H =
-0.14706 (1-0.2s)
-------------------------------------------
(1+0.1s) (1 + 1.029(0.1715s) + (0.1715s)^2)
You can read the time constants and damping ratios from the display as follows:
• Factors corresponding to real roots are displayed as (τs). The variable τ is the
time constant of the root. For example, the time constant of the zero of H is 0.2.
3-19
3 Data Manipulation
See Also
tf | zpk
Related Examples
• “Transfer Functions” on page 2-3
3-20
4
Model Interconnections
For example, you can interconnect dynamic system models of a plant G(s), a controller
C(s), sensor dynamics S(s), and a filter F(s) to construct a single model that represents
the entire closed-loop control system in the following illustration:
+
r F(s) C(s) G(s) y
-
S(s)
See Also
More About
• “Catalog of Model Interconnections” on page 4-3
4-2
Catalog of Model Interconnections
H2
N/A H1/H2 (division)
u H2-1 H1 y
4-3
4 Model Interconnections
z2 H2 w2
Arithmetic Operations
You can apply almost all arithmetic operations to dynamic system models, including those
shown below.
Operation Description
+ Addition
- Subtraction
* Multiplication
.* Element-by-element multiplication
/ Right matrix divide
\ Left matrix divide
inv Matrix inversion
' Conjugate system. For a system G, the transfer function of G' is:
4-4
See Also
Operation Description
^ Powers of a dynamic system model, as in the following syntax for
creating transfer functions:
s = tf('s');
G = 25/(s^2 + 10*s + 25);
In some cases, you might obtain better results using model interconnection commands,
such as feedback or connect, instead of model arithmetic. For example, the command T
= feedback(H1,H2) returns better results than the algebraic expression T = H1/
(1+H2*H1). The latter expression duplicates the poles of H1, which inflates the model
order and might lead to computational inaccuracy.
See Also
connect | feedback | parallel | series
Related Examples
• “Numeric Model of SISO Feedback Loop” on page 4-6
• “Multi-Loop Control System” on page 4-10
• “MIMO Control System” on page 4-19
More About
• “How the Software Determines Properties of Connected Models” on page 4-27
• “Recommended Model Type for Building Block Diagrams” on page 4-31
4-5
4 Model Interconnections
+
r F(s) C(s) G(s) y
-
S(s)
The feedback loop includes a plant G(s), a controller C(s), and a representation of sensor
dynamics, S(s). The system also includes a prefilter F(s).
1 Create model objects representing each of the components.
G = zpk([],[-1,-1],1);
C = pid(2,1.3,0.3,0.5);
S = tf(5,[1 4]);
F = tf(1,[1 1]);
The plant G is a zero-pole-gain (zpk) model with a double pole at s = –1. Model object
C is a PID controller. The models F and S are transfer functions.
2 Connect the controller and plant models.
H = G*C;
To combine models using the multiplication operator *, enter the models in reverse
order compared to the block diagram.
3
H
Construct the unfiltered closed-loop response T ( s ) = .
1 + HS
4-6
See Also
T = feedback(H,S);
T = H/(1+H*S)
This computation duplicates the poles of H, which inflates the model order and might
lead to computational inaccuracy.
4 Construct the entire closed-loop system response from r to y.
T_ry = T*F;
T_ry is a Numeric LTI Model representing the aggregate closed-loop system. T_ry does
not keep track of the coefficients of the components G, C, F, and S.
You can operate on T_ry with any Control System Toolbox control design or analysis
commands.
See Also
connect | feedback | parallel | series
Related Examples
• “Control System Model With Both Numeric and Tunable Components” on page 4-8
• “Multi-Loop Control System” on page 4-10
• “MIMO Control System” on page 4-19
More About
• “Catalog of Model Interconnections” on page 4-3
4-7
4 Model Interconnections
Suppose that the plant response is , and that the model of the sensor
dynamics is . The controller is a tunable PID controller, and the
prefilter is a low-pass filter with one tunable parameter, a.
Create models representing the plant and sensor dynamics. Because the plant and sensor
dynamics are fixed, represent them using numeric LTI models.
G = zpk([],[-1,-1],1);
S = tf(5,[1 4]);
To model the tunable components, use Control Design Blocks. Create a tunable
representation of the controller C.
C = tunablePID('C','PID');
a is a realp (real tunable parameter) object with initial value 10. Using a as a coefficient
in tf creates the tunable genss model object F.
4-8
See Also
Interconnect the models to construct a model of the complete closed-loop response from r
to y.
T = feedback(G*C,S)*F
T =
Type "ss(T)" to see the current value, "get(T)" to see all properties, and "T.Blocks" t
T.Blocks
When you create a genss model of a control system that has tunable components, you
can use tuning commands such as systune to tune the free parameters to meet design
requirements you specify.
See Also
feedback | tunablePID
More About
• “Control Design Blocks” on page 1-16
• “Dynamic System Models” on page 1-10
4-9
4 Model Interconnections
u
P y
sum1
+ e yp y1 +
ysp C Gp Dp sum3
- -
ym
+
sum2
+
dp dy
F
For more information about the Smith Predictor, see “Control of Processes with Long
Dead Time: The Smith Predictor”.
The connect command lets you construct the overall transfer function from ysp to y. To
use connect, specify the input and output channel names of the components of the block
diagram. connect automatically joins ports that have the same name, as shown in the
following figure.
error
B1 B1
connect(B1,B2)
error
B2 error B2
To build the closed loop model of the Smith Predictor system from ysp to y:
1 Create the components of the block diagram: the process model P, the predictor
model Gp, the delay model Dp, the filter F, and the PI controller C. Specify names for
the input and output channels of each model so that connect can automatically join
them to build the block diagram.
4-10
See Also
s = tf('s');
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u'; P.OutputName = 'y';
Gp = 5.6/(40.2*s+1);
Gp.InputName = 'u'; Gp.OutputName = 'yp';
Dp = exp(-93.9*s);
Dp.InputName = 'yp'; Dp.OutputName = 'y1';
F = 1/(20*s+1);
F.InputName = 'dy'; F.OutputName = 'dp';
C = pidstd(0.574,40.1);
C.Inputname = 'e'; C.OutputName = 'u';
2 Create the summing junctions needed to complete the block diagram.
The argument to sumblk is a formula that relates the input and output signals of the
summing junction. sumblk creates a summing junction with the input and output
signal names specified in the formula. For example, in sum1, the formula 'e = ysp
- ym' specifies an output signal named e, which is the difference between input
signals named ysp and ym.
3 Assemble the complete model from ysp to y.
T = connect(P,Gp,Dp,C,F,sum1,sum2,sum3,'ysp','y');
You can list the models and summing junctions in any order because connect
automatically interconnects them using their input and output channel names.
The last two arguments specify the input and output signals of the multi-loop control
structure. Thus, T is a ss model with input ysp and output y.
See Also
connect | sumblk
4-11
4 Model Interconnections
Related Examples
• “Control System Model With Both Numeric and Tunable Components” on page 4-8
• “MIMO Control System” on page 4-19
• “Mark Analysis Points in Closed-Loop Models” on page 4-13
More About
• “How the Software Determines Properties of Connected Models” on page 4-27
4-12
Mark Analysis Points in Closed-Loop Models
For this example, create a model of a Smith predictor, the SISO multiloop control system
shown in the following block diagram.
Points marked by x are analysis points to mark for this example. For instance, if you want
to calculate the step response of the closed-loop system to a disturbance at the plant
input, you can use an analysis point at u. If you want to calculate the response of the
system with one or both of the control loops open, you can use the analysis points at yp or
dp.
To build this system, first create the dynamic components of the block diagram. Specify
names for the input and output channels of each model so that connect can
automatically join them to build the block diagram.
s = tf('s');
% Process model
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u';
P.OutputName = 'y';
% Predictor model
Gp = 5.6/(40.2*s+1);
4-13
4 Model Interconnections
Gp.InputName = 'u';
Gp.OutputName = 'yp';
% Delay model
Dp = exp(-93.9*s);
Dp.InputName = 'yp';
Dp.OutputName = 'y1';
% Filter
F = 1/(20*s+1);
F.InputName = 'dy';
F.OutputName = 'dp';
% PI controller
C = pidstd(0.574,40.1);
C.Inputname = 'e';
C.OutputName = 'u';
Create the summing junctions needed to complete the block diagram. (For more
information about creating summing junctions, see sumblk).
sum1 = sumblk('e = ysp - ym');
sum2 = sumblk('ym = yp + dp');
sum3 = sumblk('dy = y - y1');
T =
Type "ss(T)" to see the current value, "get(T)" to see all properties, and "T.Blocks" t
When you use the APs input argument, the connect command automatically inserts an
AnalysisPoint block into the generalized state-space (genss) model, T. The
automatically generated block is named AnalysisPoints_. The three channels of
AnalysisPoints_ correspond to the three locations specified in the APs argument to
4-14
Mark Analysis Points in Closed-Loop Models
the connect command. Use getPoints to see a list of the available analysis points in
the model.
getPoints(T)
ans =
{'dp'}
{'u' }
{'yp'}
Use these locations as inputs, outputs, or loop openings when you extract responses from
the model. For example, extract and plot the response at the system output to a step
disturbance at the plant input, u.
Tp = getIOTransfer(T,'u','y');
stepplot(Tp)
4-15
4 Model Interconnections
Similarly, calculate the open-loop response of the plant and controller by opening both
feedback loops.
openings = {'dp','yp'};
L = getIOTransfer(T,'ysp','y',openings);
bodeplot(L)
4-16
See Also
When you create a control system model, you can create an AnalysisPoint block
explicitly and assign input and output names to it. You can then include it in the input
arguments to connect as one of the blocks to combine. However, using the APs
argument to connect as illustrated in this example is a simpler way to mark points of
interest when building control system models.
See Also
AnalysisPoint | connect | sumblk
4-17
4 Model Interconnections
Related Examples
• “Control System with Multichannel Analysis Points” on page 2-85
More About
• “Mark Signals of Interest for Control System Analysis and Design” on page 2-89
4-18
MIMO Control System
pL L
+ e CL
r D G y
pV CV
-
V
The plant G is a distillation column with two inputs and two outputs. The two inputs are
the reflux L and boilup V. The two outputs are the concentrations of two chemicals,
represented by the vector signal y = [y1,y2]. You can represent this plant model as:
1 È 87 .8 -86 .4 ˘
G ( s) = Í ˙.
75s + 1 Î108.2 -109 .6 ˚
The vector setpoint signal r = [r1,r2] specifies the desired concentrations of the two
chemicals. The vector error signal e represents the input to D, a static 2-by-2 decoupling
matrix. CL and CV represent independent PI controllers that control the two inputs of G.
When you construct the closed-loop model, connect uses the input and output
names to form connections between the block diagram components. Therefore, you
must assign names to the inputs and outputs of the transfer function G in either of the
following ways: .
4-19
4 Model Interconnections
• You can assign input and output names to individual signals by specifying signal
names in a cell array, as in G.InputName = {'L','V'}
• Alternatively, you can use vector signal naming, which the software automatically
expands. For example, the command G.OutputName = 'y' assigns the names
'y(1)' and 'y(2)' to the outputs of G.
2 Create tunable Control Design Blocks representing the decoupling matrix D and the
PI controllers CL and CV.
D = tunableGain('Decoupler',eye(2));
D.u = 'e';
D.y = {'pL','pV'};
Note u and y are shorthand notations for the InputName and OutputName
properties, respectively. Thus, for example, entering:
D.u = 'e';
D.y = {'pL','pV'};
is equivalent to entering:
D.InputName = 'e';
D.OutputName = {'pL','pV'};
The summing junction produces the error signals e by taking the difference between
r and y.
Sum represents the transfer function for the summing junction described by the
formula 'e = r - y'. The second argument to sumblk specifies that the inputs and
outputs of Sum are each vector signals of length 2. The software therefore
automatically assigns the signal names {'r(1)','r(2)','y(1)','y(2)'} to
Sum.InputName and {'e(1)','e(2)'} to Sum.OutputName.
4-20
See Also
CLry = connect(G,D,C_L,C_V,Sum,'r','y');
The arguments to the connect function include all the components of the closed-loop
system, in any order. connect automatically combines the components using the
input and output names to join signals.
The last two arguments to connect specify the output and input signals of the
closed-loop model, respectively. The resulting genss model CLry has two-inputs and
two outputs.
See Also
connect | sumblk
Related Examples
• “Control System Model With Both Numeric and Tunable Components” on page 4-8
• “Multi-Loop Control System” on page 4-10
• “MIMO Control System” on page 4-19
More About
• “Catalog of Model Interconnections” on page 4-3
4-21
4 Model Interconnections
In this example, you obtain the response from Azref to Az of the MIMO feedback loop of
the following block diagram.
You can compute the closed-loop response using one of the following three approaches:
You can use whichever of these approaches is most convenient for your application.
Load the plant Aerodyn and the controller Autopilot into the MATLAB® workspace.
These models are stored in the datafile MIMOfeedback.mat.
load(fullfile(matlabroot,'examples','control','MIMOfeedback.mat'))
4-22
MIMO Feedback Loop
T1 = connect(Autopilot,Aerodyn,'Azref','Az');
The connect function combines the models by joining the inputs and outputs that have
matching names. The last two arguments to connect specify the input and output signals
of the resulting model. Therefore, T1 is a state-space model with input Azref and output
Az. The connect function ignores the other inputs and outputs in Autopilot and
Aerodyn.
When you use the feedback function, think of the closed-loop system as a feedback
interconnection between an open-loop plant-controller combination L and a diagonal
unity-gain feedback element K. The following block diagram shows this interconnection.
L = series(Autopilot,Aerodyn,'Fin');
FeedbackChannels = {'Alpha','Mach','Az','q'};
K = ss(eye(4),'InputName',FeedbackChannels,...
'OutputName',FeedbackChannels);
4-23
4 Model Interconnections
T2 = feedback(L,K,'name',+1);
Compute the closed-loop response from Azref to Az using feedback, using indices to
specify the interconnections between Aerodyn and Autopilot.
L = series(Autopilot,Aerodyn,1,4);
K = ss(eye(4));
T3 = feedback(L,K,[1 2 3 4],[4 3 6 5],+1);
Compare the step response from Azref to Az to confirm that the three approaches yield
the same results.
step(T1,T2('Az','Azref'),T3(6,5),2)
4-24
See Also
See Also
connect | feedback
Related Examples
• “Multi-Loop Control System” on page 4-10
• “MIMO Control System” on page 4-19
4-25
4 Model Interconnections
More About
• “How the Software Determines Properties of Connected Models” on page 4-27
4-26
How the Software Determines Properties of Connected Models
4-27
4 Model Interconnections
See Also
More About
• “Rules That Determine Model Type” on page 4-29
4-28
Rules That Determine Model Type
P = ss([-0.8,0.4;0.4,-1.0],[-3.0;1.4],[0.3,0],0);
C = pid(-0.13,-0.61);
CL = feedback(P*C,1)
The ss model has the highest precedence among Numeric LTI models. Therefore,
combining P and C with any model interconnection command returns an ss model.
Combining Numeric LTI models with Generalized LTI models on page 1-16 or with Control
Design Blocks on page 1-16 results in Generalized LTI models.
F = tunableTF('F',0,1);
CLF = F*CL
Note The software automatically converts all models to the resulting model type before
performing the connection operation.
See Also
connect | feedback | parallel | series
Related Examples
• “Numeric Model of SISO Feedback Loop” on page 4-6
4-29
4 Model Interconnections
More About
• “How the Software Determines Properties of Connected Models” on page 4-27
• “Recommended Model Type for Building Block Diagrams” on page 4-31
4-30
Recommended Model Type for Building Block Diagrams
You can represent block diagram components with any model type. However, certain
connection operations yield better numerical accuracy for models in ss form.
For example, interconnect two models in series using different model types to see how
different representations introduce numerical inaccuracies.
Load the models Pd and Cd. These models are ninth-order and second-order discrete-time
transfer functions, respectively.
load numdemo Pd Cd
Compute the open-loop transfer function L = Pd*Cd using the tf, zpk, ss, and frd
representations.
Ltf = Pd*Cd;
Lzp = zpk(Pd)*Cd;
Lss = ss(Pd)*Cd;
w = logspace(-1,3,100);
Lfrd = frd(Pd,w)*Cd;
Plot the magnitude of the frequency response to compare the four representations.
bodemag(Ltf,Lzp,Lss,Lfrd)
legend('tf','zpk','ss','frd')
4-31
4 Model Interconnections
See Also
More About
• “Rules That Determine Model Type” on page 4-29
4-32
5
Model Transformation
In general, you can convert from any model type to any other. However, there are a few
limitations. For example, you cannot convert:
• frd models to analytic model types such as ss, tf, or zpk (unless you perform system
identification with System Identification Toolbox software).
• ss models with internal delays to tf or zpk.
You can convert between Numeric LTI models and Generalized LTI models.
• Converting a Generalized LTI model to a Numeric LTI model evaluates any Control
Design Blocks at their current (nominal) value.
• Converting a Numeric LTI model to a Generalized LTI model creates a Generalized LTI
model with an empty Blocks property.
sys = ss(0,1,1,0)
[num,den] = tfdata(sys)
tfdata automatically converts the state-space model sys to transfer function form to
return numerator and denominator data.
Conversions to state-space form are not uniquely defined. For this reason, automatic
conversions to state space do not occur when the result depends on the choice of state
5-2
See Also
coordinates. For example, the initial and kalman commands require state-space
models.
• The accuracy of computations using high-order transfer functions (tf or zpk models)
is sometimes poor, particularly for MIMO or high-order systems. Conversions to a
transfer function representation can incur a loss of accuracy.
• When you convert tf or zpk models to state space using ss, the software
automatically performs balancing and scaling operations. Balancing and scaling
improves the numeric accuracy of computations involving the model. For more
information about balancing and scaling state-space models, see “Scaling State-Space
Models” on page 23-2.
In addition, converting back and forth between model types can introduce additional
states or orders, or introduce numeric inaccuracies. For example, conversions to state
space are not uniquely defined, and are not guaranteed to produce a minimal realization
for MIMO models. For a given state-space model sys,
ss(tf(sys))
can return a model with different state-space matrices, or even a different number of
states in the MIMO case.
See Also
frd | pid | ss | tf | zpk
Related Examples
• “Convert From One Model Type to Another” on page 5-4
5-3
5 Model Transformation
In general, you can convert a model from one type to another type using the model-
creation command for the target type. For example, you can use the tf command to
convert an ss model to transfer function form, or use the ss command to convert a zpk
model to state-space form.
pid_sys = pid(1,1.5,3)
pid_sys =
1
Kp + Ki * --- + Kd * s
s
with Kp = 1, Ki = 1.5, Kd = 3
C = tf(pid_sys)
C =
3 s^2 + s + 1.5
---------------
s
C is a tf representation of pid_sys. C has the same dynamics as pid_sys, but stores the
dynamic parameters as transfer-function numerator and denominator coefficients instead
of proportional, integral, and derivative gains.
You can similarly convert transfer function models to pid models, provided the tf model
object represents a parallel-form PID controller with .
5-4
See Also
In general, you can use the technique of this example to convert any type of model to
another type of model. For more specific information about converting to a particular
model type, see the reference page for that model type.
See Also
frd | pid | ss | tf | zpk
More About
• “Conversion Between Model Types” on page 5-2
5-5
5 Model Transformation
a = realp('a',10);
F = tf(a,[1 a]);
Typically, once of you have a generalized model, you tune the parameters of the model
using a tuning command such as systune. For this example, instead of tuning the model,
manually change the value of the tunable component of F.
F.Blocks.a.Value = 5;
Get the current value of the generalized model by converting it to a numeric model.
F_cur_val = tf(F)
F_cur_val =
5
-----
s + 5
5-6
See Also
See Also
realp | showBlockValue | tf
More About
• “Models with Tunable Coefficients” on page 1-19
• “Conversion Between Model Types” on page 5-2
• “Convert From One Model Type to Another” on page 5-4
5-7
5 Model Transformation
Obtain a 2-DOF PID controller. For this example, create a plant model, and tune a 2-DOF
PID controller for it.
C2 is a pid2 controller object. The control architecture for C2 is as shown in the following
illustration.
This control system can be equivalently represented in several other architectures that
use only SISO components. In the feedforward configuration, the 2-DOF controller is
represented as a SISO PID controller and a feedforward compensator.
[Cff,Xff] = getComponents(C2,'feedforward')
5-8
Decompose a 2-DOF PID Controller into SISO Components
Cff =
1 s
Kp + Ki * --- + Kd * --------
s Tf*s+1
Xff =
-10.898 (s+0.2838)
------------------
(s+8.181)
This command returns the SISO PID controller Cff as a pid object. The feedforward
compensator X is returned as a zpk object.
Tff = G*(Cff+Xff)*feedback(1,G*Cff);
Decompose C2 using the feedback configuration and construct that closed-loop system.
5-9
5 Model Transformation
In the filter configuration, the 2-DOF controller is represented as a SISO PID controller
and prefilter on the reference signal.
Decompose C2 using the filter configuration. Construct that closed-loop system as well.
Construct the closed-loop system for the original 2-DOF controller, C2. To do so, convert
C2 to a two-input, one-output transfer function, and use array indexing to access the
channels.
Ctf = tf(C2);
Cr = Ctf(1);
Cy = Ctf(2);
T = Cr*feedback(G,Cy,+1);
stepplot(T,Tff,Tfb,Tfr)
legend('2-DOF','feedforward','feedback','filter','Location','Southeast')
5-10
See Also
The plots coincide, demonstrating that all the systems are equivalent.
Using a 2-DOF PID controller can yield improved performance compared to a 1-DOF
controller. For more information, see “Tune 2-DOF PID Controller (Command Line)” on
page 11-15.
See Also
getComponents | pid2 | pidstd2
5-11
5 Model Transformation
Related Examples
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
5-12
Discretize a Compensator
Discretize a Compensator
This example shows how to convert a compensator from continuous to discrete time using
several discretization methods, to identify a method that yields a good match in the
frequency domain.
You might design a compensator in continuous time, and then need to convert it to
discrete time for a digital implementation. When you do so, you want the discretization to
preserve frequency-domain characteristics that are essential to your performance and
stability requirements.
One valid controller for this system includes a notch filter in series with an integrator.
Create a model of this controller.
notch = tf([1,0.5,9],[1,5,9]);
integ = pid(0,0.34);
C = integ*notch;
bodeplot(C)
5-13
5 Model Transformation
The notch filter centered at 3 rad/s counteracts the effect of the resonance in G. This
configuration allows higher loop gain for a faster overall response.
Cdz = c2d(C,0.5);
The c2d command supports several different discretization methods. Since this command
does not specify a method, c2d uses the default method, Zero-Order Hold (ZOH). In the
ZOH method, the time-domain response of the discretized compensator matches the
continuous-time response at each time step.
5-14
Discretize a Compensator
The discretized controller Cdz has a sample time of 0.5 s. In practice, the sample time you
choose might be constrained by the system in which you implement your controller, or by
the bandwidth of your control system.
bodeplot(C,Cdz)
legend('C','Cdz');
The vertical line marks the Nyquist frequency, , where is the sample time. Near
the Nyquist frequency, the response of the discretized compensator is distorted relative to
the continuous-time response. As a result, the discretized notched filter may not properly
counteract the plant resonance.
5-15
5 Model Transformation
To fix this, try discretizing the compensator using the Tustin method and compare to the
ZOH result. The Tustin discretization method often yields a better match in the frequency
domain than the ZOH method.
Cdt = c2d(C,0.5,'tustin');
plotopts = bodeoptions;
plotopts.Ylim = {[-60,40],[-225,0]};
bodeplot(C,Cdz,Cdt,plotopts)
legend('C','Cdz','Cdt')
The Tustin method preserves the depth of the notch. However, the method introduces a
frequency shift that is unacceptable for many applications. You can remedy the frequency
shift by specifying the notch frequency as the prewarping frequency in the Tustin
transform.
5-16
Discretize a Compensator
Discretize the compensator using the Tustin method with frequency prewarping, and
compare the results.
discopts = c2dOptions('Method','tustin','PrewarpFrequency',3.0);
Cdtp = c2d(C,0.5,discopts);
bodeplot(C,Cdt,Cdtp,plotopts)
legend('C','Cdt','Cdtp')
5-17
5 Model Transformation
Using the Tustin method with frequency prewarping yields a better-matching frequency
response than Tustin without prewarping.
See Also
c2d | c2dOptions
More About
• “Continuous-Discrete Conversion Methods” on page 5-26
• “Improve Accuracy of Discretized System with Time Delay” on page 5-19
5-18
Improve Accuracy of Discretized System with Time Delay
For systems with time delays that are not integer multiples of the sample time, the
Tustin and Matched methods by default round the time delays to the nearest multiple of
the sample time. To improve the accuracy of these methods for such systems, c2d can
optionally approximate the fractional portion of the time delay by a discrete-time all-pass
filter (a Thiran filter). In this example, discretize the system both without and with an
approximation of the fractional portion of the delay and compare the results.
G = tf(1,[1,0.2,4],'ioDelay',2.5);
discopts = c2dOptions('Method','tustin','PrewarpFrequency',2);
Gt = c2d(G,1,discopts)
Warning: Rounding delays to the nearest multiple of the sampling period. For more accur
Gt =
The software warns you that it rounds the fractional time delay to the nearest multiple of
the sample time. In this example, the time delay of 2.5 times the sample time (2.5 s)
converts to an additional factor of z^(-3) in Gt.
5-19
5 Model Transformation
plotopts = bodeoptions;
plotopts.Ylim = {[-100,20],[-1080,0]};
bodeplot(G,Gt,plotopts);
legend('G','Gt')
There is a phase lag between the discretized system Gt and the continuous-time system G,
which grows as the frequency approaches the Nyquist frequency. This phase lag is largely
due to the rounding of the fractional time delay. In this example, the fractional time delay
is half the sample time.
5-20
Improve Accuracy of Discretized System with Time Delay
discopts.FractDelayApproxOrder = 3;
Gtf = c2d(G,1,discopts);
The FractDelayApproxOrder option specifies the order of the Thiran filter that
approximates the fractional portion of the delay. The other options in discopts are
unchanged. Thus Gtf is a Tustin discretization of G with prewarp at 2 rad/s.
plotopts.PhaseMatching = 'on';
bodeplot(G,Gt,Gtf,plotopts);
legend('G','Gt','Gtf','Location','SouthWest')
5-21
5 Model Transformation
The magnitudes of Gt and Gtf are identical. However, the phase of Gtf provides a better
match to the phase of the continuous-time system through the resonance. As the
frequency approaches the Nyquist frequency, this phase match deteriorates. A higher-
order approximation of the fractional delay would improve the phase matching closer to
the Nyquist frequencies. However, each additional order of approximation adds an
additional order (or state) to the discretized system.
If your application requires accurate frequency-matching near the Nyquist frequency, use
c2dOptions to make c2d approximate the fractional portion of the time delay as a
Thiran filter.
See Also
c2d | c2dOptions | thiran
More About
• “Continuous-Discrete Conversion Methods” on page 5-26
• “Discretize a Compensator” on page 5-13
5-22
Convert Discrete-Time System to Continuous Time
Convert the following second-order discrete-time system to continuous time using the
zero-order hold (ZOH) method:
G = zpk(-0.5,[-2,5],1,0.1);
Gcz = d2c(G)
Warning: The model order was increased to handle real negative poles.
Gcz =
When you call d2c without specifying a method, the function uses ZOH by default. The
ZOH interpolation method increases the model order for systems that have real negative
poles. This order increase occurs because the interpolation algorithm maps real negative
poles in the domain to pairs of complex conjugate poles in the domain.
Gct = d2c(G,'tustin')
Gct =
5-23
5 Model Transformation
bode(G,Gcz,Gct)
legend('G','Gcz','Gct')
In this case, the Tustin method provides a better frequency-domain match between the
discrete system and the interpolation. However, the Tustin interpolation method is
5-24
See Also
undefined for systems with poles at z = -1 (integrators), and is ill-conditioned for systems
with poles near z = 1.
See Also
d2c | d2cOptions
More About
• “Continuous-Discrete Conversion Methods” on page 5-26
• “Discretize a Compensator” on page 5-13
5-25
5 Model Transformation
5-26
Continuous-Discrete Conversion Methods
Zero-Order Hold
The Zero-Order Hold (ZOH) method provides an exact match between the continuous-
and discrete-time systems in the time domain for staircase inputs.
The ZOH block generates the continuous-time input signal u(t) by holding each sample
value u(k) constant over one sample period:
u ( t ) = u [ k] , kTs £ t £ ( k + 1 ) Ts
The signal u(t) is the input to the continuous system H(s). The output y[k] results from
sampling y(t) every Ts seconds.
Conversely, given a discrete system Hd(z), d2c produces a continuous system H(s). The
ZOH discretization of H(s) coincides with Hd(z).
5-27
5 Model Transformation
You can use the ZOH method to discretize SISO or MIMO continuous-time models with
time delays. The ZOH method yields an exact discretization for systems with input delays,
output delays, or transfer delays.
For systems with internal delays (delays in feedback loops), the ZOH method results in
approximate discretizations. The following figure illustrates a system with an internal
delay.
H(s)
e-ts
For such systems, c2d performs the following actions to compute an approximate ZOH
discretization:
1
Decomposes the delay τ as t = kTs + r with 0 £ r < Ts .
2
Absorbs the fractional delay r into H(s).
3 Discretizes H(s) to H(z).
4 Represents the integer portion of the delay kTs as an internal discrete-time delay z–k.
The final discretized model appears in the following figure:
H(z)
H(s) e-sr
z-k
First-Order Hold
The First-Order Hold (FOH) method provides an exact match between the continuous-
and discrete-time systems in the time domain for piecewise linear inputs.
5-28
Continuous-Discrete Conversion Methods
FOH differs from ZOH by the underlying hold mechanism. To turn the input samples u[k]
into a continuous input u(t), FOH uses linear interpolation between samples:
t - kTs
u ( t ) = u [ k] +
Ts
(u [k + 1] - u [k]) , kTs £ t £ (k + 1) Ts
In general, this method is more accurate than ZOH for systems driven by smooth inputs.
This FOH method differs from standard causal FOH and is more appropriately called
triangle approximation (see [2], p. 228). The method is also known as ramp-invariant
approximation.
You can use the FOH method to discretize SISO or MIMO continuous-time models with
time delays. The FOH method handles time delays in the same way as the ZOH method.
See “ZOH Method for Systems with Time Delays” on page 5-28.
Impulse-Invariant Mapping
The impulse-invariant mapping produces a discrete-time model with the same impulse
response as the continuous time system. For example, compare the impulse response of a
first-order continuous system with the impulse-invariant discretization:
G = tf(1,[1,1]);
Gd1 = c2d(G,0.01,'impulse');
impulse(G,Gd1)
5-29
5 Model Transformation
The impulse response plot shows that the impulse responses of the continuous and
discretized systems match.
Tustin Approximation
The Tustin or bilinear approximation yields the best frequency-domain match between the
continuous-time and discretized systems. This method relates the s-domain and z-domain
transfer functions using the approximation:
5-30
Continuous-Discrete Conversion Methods
sTs 1 + sTs / 2
z=e ª .
1 - sTs / 2
In c2d conversions, the discretization Hd(z) of a continuous transfer function H(s) is:
2 z-1
H d ( z ) = H ( s¢ ) , s¢ =
Ts z + 1
1 + sTs / 2
H ( s ) = Hd ( z¢ ) , z¢ =
1 - sTs / 2
When you convert a state-space model using the Tustin method, the states are not
preserved. The state transformation depends upon the state-space matrices and whether
the system has time delays. For example, for an explicit (E = I) continuous-time model
with no time delays, the state vector w[k] of the discretized model is related to the
continuous-time state vector x(t) by:
Ê T ˆ T T
w [ kTs ] = Á I - A s ˜ x ( kTs ) - s Bu ( kTs ) = x ( kTs ) - s ( Ax ( kTs ) + Bu ( kTs ) ) .
Ë 2 ¯ 2 2
Ts is the sample time of the discrete-time model. A and B are state-space matrices of the
continuous-time model.
If your system has important dynamics at a particular frequency that you want the
transformation to preserve, you can use the Tustin method with frequency prewarping.
This method ensures a match between the continuous- and discrete-time responses at the
prewarp frequency.
The Tustin approximation with frequency prewarping uses the following transformation of
variables:
w z-1
H d ( z ) = H ( s¢ ) , s¢ =
tan (w Ts / 2 ) z + 1
5-31
5 Model Transformation
This change of variable ensures the matching of the continuous- and discrete-time
frequency responses at the prewarp frequency ω, because of the following
correspondence:
(
H ( jw ) = H d e jw Ts )
Tustin Approximation for Systems with Time Delays
You can use the Tustin approximation to discretize SISO or MIMO continuous-time models
with time delays.
By default, the Tustin method rounds any time delay to the nearest multiple of the sample
time. Therefore, for any time delay tau, the integer portion of the delay, k*Ts, maps to a
delay of k sampling periods in the discretized model. This approach ignores the residual
fractional delay, tau - k*Ts.
You can to approximate the fractional portion of the delay by a discrete all-pass filter
(Thiran filter) of specified order. To do so, use the FractDelayApproxOrder option of
c2dOptions. See “Improve Accuracy of Discretized System with Time Delay” on page 5-
19 for an example.
To understand how the Tustin method handles systems with time delays, consider the
following SISO state-space model G(s). The model has input delay τi, output delay τo, and
internal delay τ.
G(s)
e-tis e-tos
H(s)
e-ts
The following figure shows the general result of discretizing G(s) using the Tustin
method.
5-32
Continuous-Discrete Conversion Methods
Gd(z)
z-m F(z)
By default, c2d converts the time delays to pure integer time delays. The c2d command
computes the integer delays by rounding each time delay to the nearest multiple of the
sample time Ts. Thus, in the default case, mi = round(τi/Ts), mo = round(τo/Ts), and m =
round(τ/Ts).. Also in this case, Fi(z) = Fo(z) = F(z) = 1.
The Thiran filters add additional states to the model. The maximum number of additional
states for each delay is FractDelayApproxOrder.
For example, for the input delay τi, the order of the Thiran filter Fi(z) is:
If ceil(τi/Ts) < FractDelayApproxOrder, the Thiran filter Fi(z) approximates the entire
input delay τi. If ceil(τi/Ts) > FractDelayApproxOrder, the Thiran filter only
approximates a portion of the input delay. In that case, c2d represents the remainder of
the input delay as a chain of unit delays z–mi, where
mi = ceil(τi/Ts) – FractDelayApproxOrder
When you discretizetf and zpk models using the Tustin method, c2d first aggregates all
input, output, and transfer delays into a single transfer delay τTOT for each channel. c2d
then approximates τTOT as a Thiran filter and a chain of unit delays in the same way as
described for each of the time delays in ss models.
5-33
5 Model Transformation
For more information about Thiran filters, see the thiran reference page and [4].
zi = esi Ts
where:
You can use zero-pole matching to discretize SISO continuous-time models with time
delay, except that the method does not support ss models with internal delays. The zero-
pole matching method handles time delays in the same way as the Tustin approximation.
See “Tustin Approximation for Systems with Time Delays” on page 5-32.
Least Squares
The least squares method minimizes the error between the frequency responses of the
continuous-time and discrete-time systems up to the Nyquist frequency using a vector-
fitting optimization approach. This method is useful when you want to capture fast system
dynamics but must use a larger sample time, for example, when computational resources
are limited.
This method is supported only by the c2d function and only for SISO systems.
As with Tustin approximation and zero-pole matching, the least squares method provides
a good match between the frequency responses of the original continuous-time system
and the converted discrete-time system. However, when using the least squares method
with:
5-34
See Also
• The same sample time as Tustin approximation or zero-pole matching, you get a
smaller difference between the continuous-time and discrete-time frequency
responses.
• A lower sample time than what you would use with Tustin approximation or zero-pole
matching, you can still get a result that meets your requirements. Doing so is useful if
computational resources are limited, since the slower sample time means that the
processor must do less work.
References
[1] Åström, K.J. and B. Wittenmark, Computer-Controlled Systems: Theory and Design,
Prentice-Hall, 1990, pp. 48-52.
[2] Franklin, G.F., Powell, D.J., and Workman, M.L., Digital Control of Dynamic Systems
(3rd Edition), Prentice Hall, 1997.
[3] Smith, J.O. III, "Impulse Invariant Method", Physical Audio Signal Processing, August
2007. http://www.dsprelated.com/dspbooks/pasp/
Impulse_Invariant_Method.html.
[4] T. Laakso, V. Valimaki, "Splitting the Unit Delay", IEEE Signal Processing Magazine,
Vol. 13, No. 1, p.30-60, 1996.
See Also
c2d | c2dOptions | d2c | d2cOptions | d2d | d2dOptions | thiran
Related Examples
• “Discretize a Compensator” on page 5-13
• “Improve Accuracy of Discretized System with Time Delay” on page 5-19
• “Convert Discrete-Time System to Continuous Time” on page 5-23
5-35
5 Model Transformation
Upsampling a system can be useful, for example, when you need to implement a digital
controller at a faster rate than you originally designed it for.
G = tf([1,0.4],[1,-0.7],0.3);
G_d2d = d2d(G,0.1)
G_d2d =
z - 0.4769
----------
z - 0.8879
By default, d2d uses the zero-order-hold (ZOH) method to resample the system. The
resampled system has the same order as G.
G_up = upsample(G,3)
G_up =
z^3 + 0.4
---------
z^3 - 0.7
5-36
Upsample Discrete-Time System
The second input, 3, tells upsample to resample G at a sample time three times faster
than the sample time of G. This input to upsample must be an integer.
Compare the step responses of the original model G with the resampled models G_d2d
and G_up.
step(G,'-r',G_d2d,':g',G_up,'--b')
legend('G','d2d','upsample','Location','SouthEast')
5-37
5 Model Transformation
The step response of the upsampled model G_up matches exactly the step response of the
original model G. The response of the resampled model G_d2d matches only at every third
sample.
Compare the frequency response of the original model with the resampled models.
bode(G,'-r',G_d2d,':g',G_up,'--b')
legend('G','d2d','upsample','Location','SouthWest')
In the frequency domain as well, the model G_up created with the upsample command
matches the original model exactly up to the Nyquist frequency of the original model.
Using upsample provides a better match than d2d in both the time and frequency
domains. However, upsample increases the model order, which can be undesirable.
5-38
See Also
Additionally, upsample is only available where the original sample time is an integer
multiple of the new sample time.
See Also
d2d | d2dOptions | upsample
More About
• “Choosing a Resampling Command” on page 5-40
5-39
5 Model Transformation
See Also
d2d | d2dOptions | upsample
Related Examples
• “Upsample Discrete-Time System” on page 5-36
5-40
6
Model Simplification
• You are working with a relatively high-order model obtained from linearizing a
Simulink model, performing a finite-element calculation, interconnecting model
elements, or other source.
• You want to improve the simulation speed of a Simulink model at a certain operating
point. In that case, you can linearize a portion of the model at that operating point and
compute a reduced-order simplification or approximation of the linearized model. You
can then replace the portion of the model with an LTI Block containing the reduced-
order model.
• You design a high-order controller that you want to implement as a lower-order
controller, such as a PID controller. For example, controller design using Linear-
Quadratic-Gaussian methods or H∞ synthesis techniques can yield a high-order result.
In this case, you can try reducing the plant order before synthesis, reducing the
controller order after synthesis, or both.
• You want to simplify a model obtained by identification with System Identification
Toolbox software.
The following diagram illustrates the relationship between model reduction and control
design.
6-2
Model Reduction Basics
Plant reduction
G GR
Controller design
Controller design
Controller reduction
C CR
Higher Lower
Order Order
• Discarding states that do not contribute to the system dynamics, such as structurally
disconnected states or canceling pole-zero pairs.
• Discarding low-energy states that contribute relatively little to system dynamics.
• Focusing on a particular frequency region and discarding dynamics outside that
region. For example, if your control bandwidth is limited by actuator dynamics,
discard higher-frequency dynamics.
In any case, when you reduce model order, you want to preserve model characteristics
that are important for your application. Whenever you compute a reduced-order model,
verify that the reduced model preserves time-domain or frequency-domain behavior that
you care about. For example, for control design, it is useful to verify that the reduced
closed-loop system is stable. It is also useful to check that the reduced open-loop transfer
function CRGR adequately matches the original models where the open-loop gain GC is
close to 1 (in the gain crossover region).
6-3
6 Model Simplification
Sometimes, approximation can yield better results, even if the model looks like a good
candidate for simplification. For example, models with near pole-zero cancellations are
sometimes better reduced by approximation than simplification. Similarly, using balred
to reduce state-space models can yield more accurate results than minreal.
When you use a reduced-order model, always verify that the simplification or
approximation preserves model characteristics that are important for your application.
For example, compare the frequency responses of the original and reduced models using
6-4
See Also
bodeplot or sigmaplot. Or, compare the open-loop responses for the original and
reduced plant and controller models.
See Also
Apps
Model Reducer
Functions
balred | freqsep | minreal | sminreal
Related Examples
• “Balanced Truncation Model Reduction” on page 6-17
• “Mode-Selection Model Reduction” on page 6-57
• “Pole-Zero Simplification” on page 6-47
6-5
6 Model Simplification
This example uses a model of the Los Angeles University Hospital building. The building
has eight floors, each with three degrees of freedom: two displacements and one rotation.
The input-output relationship for any one of these displacements is represented as a 48-
state model, where each state represents a displacement or its rate of change (velocity).
Load the building model and open Model Reducer with that model.
load build.mat
modelReducer(G)
Select the model in the Data Browser to display some information about the model in the
Preview section. Double-click the model to see more detailed information.
6-6
Reduce Model Order Using the Model Reducer App
Model Reducer has three model reduction methods: Balanced Truncation, Mode
Selection, and Pole/Zero Simplification. For this example, click Balanced Truncation.
6-7
6 Model Simplification
Model Reducer opens the Balanced Truncation tab and automatically generates a
reduced-order model. The top plot compares the original and reduced model in the
frequency domain. The bottom plot shows the energy contribution of each state, where
the states are sorted from high energy to low energy. The order of the reduced model, 14,
is highlighted in the bar chart. In the reduced model, all states with lower energy
contribution than this one are discarded.
Suppose that you want to preserve the first, second, and third peaks of the model
response, around 5.2 rad/s, 13 rad/s, and 25 rad/s. Try other model orders to see whether
6-8
Reduce Model Order Using the Model Reducer App
you can achieve this goal with a lower model order. Compute a 5th-order and a 10th-order
approximation in one of the following ways:
Model Reducer computes two new reduced-order models and displays them on the
response plot with the original model G. To examine the three peaks more closely, Zoom in
on the relevant frequency range. The 10th-order model captures the three peaks
successfully, while the 5th-order model only approximates the first two peaks. (For
information about zooming and other interactions with the analysis plots, see “Visualize
Reduced-Order Models in the Model Reducer App” on page 6-67.)
6-9
6 Model Simplification
In addition to the frequency response plot of all three models, Model Reducer lets you
examine the absolute and relative error between the original and reduced models. Select
Absolute error plot to see the difference between the building and reduced models.
6-10
Reduce Model Order Using the Model Reducer App
The 5th-order reduced model has at most -60dB error in the frequency region of the first
two peaks, below about 30 rad/s. The error increases at higher frequencies. The 10th-
order reduced model has smaller error over all frequencies.
Store the reduced models in the Data Browser by clicking Create Reduced Model. The
5th-order and 10th-order reduced models appear in the Data Browser with names
GReduced5 and Greduced10.
You can continue to change the model-reduction parameters and generate additional
reduced models. As you do so, GReduced5 and Greduced10 remain unchanged in the
Data Browser.
6-11
6 Model Simplification
You can also focus the balanced truncation on the model dynamics in a particular
frequency interval. For example, approximate only the second peak of the building model
around 13 rad/s. First, select the Model response plot to see the Bode plots of models.
Then check Select frequency range checkbox. Model Reducer analyzes state
contributions in the highlighted frequency interval only.
You can drag the boundaries to change the frequency range interactively. As you change
the frequency interval, the Hankel Singular Value plot reflects the changes in the energy
contributions of the states.
Enter the frequency limits [10 22] into the text box next to Select frequency range.
The 5th-order reduced model captures the essential dynamics. The 10th-order model has
almost the same dynamics as the original building model within this frequency range.
6-12
Reduce Model Order Using the Model Reducer App
Optionally, store these additional models in the Data Browser by clicking Create
Reduced Model.
You can compare time-domain responses of the stored reduced models and the original in
the Plots tab. In the Data Browser, control-click to select the models you want to
compare, G, GReduced5, and GReduced10. Then, click Step. Model Reducer creates a
step plot with all three models.
6-13
6 Model Simplification
Zooming on the transient behavior of this plot shows that GReduced10 captures the time
domain behavior of the original model well. However, the response of GReduced5
deviates from the original model after about 3 seconds.
Comparison of the reduced and original models in the time and frequency domains shows
that GReduced10 adequately captures the dynamics of interest. Export that model to the
MATLAB® workspace for further analysis and design. In the Model Reducer tab, click
Export Model. Clear the check boxes for G and Greduced5, and click Export to export
Greduced10.
6-14
See Also
See Also
Model Reducer
Related Examples
• “Model Reduction Basics” on page 6-2
• “Balanced Truncation Model Reduction” on page 6-17
• “Pole-Zero Simplification” on page 6-47
6-15
6 Model Simplification
6-16
Balanced Truncation Model Reduction
For more general information about model reduction, see “Model Reduction Basics” on
page 6-2.
1 Open the app, and import an LTI model to reduce. For instance, suppose that there is
a model named build in the MATLAB workspace. The following command opens
Model Reducer and imports the model.
modelReducer(build)
2
In the Data Browser, select the model to reduce. Click Balanced Truncation.
6-17
6 Model Simplification
In the Balanced Truncation tab, Model Reducer displays a plot of the frequency
response of the original model and a reduced version of the model. The frequency
response is a Bode plot for SISO models, and a singular-value plot for MIMO models.
The app also displays a Hankel singular-value plot of the original model.
6-18
Balanced Truncation Model Reduction
The Hankel singular-value plot shows the relative energy contributions of each state
in the system. Model Reducer computes an initial reduced-order model based on
these values. The highlighted bar is the lowest-energy state in the initial reduced-
order model. Model Reducer discards states that have lower Hankel singular values
than the highlighted bar.
6-19
6 Model Simplification
3 Try different reduced-model orders to find the lowest-order model that preserves the
dynamics that are important for your application. To specify different orders, either:
• Enter model orders in the Reduced model orders field. You can enter a single
integer or an array of integers, such as 10:14 or [8,11,12].
• Click a bar on the Hankel singular-value plot to specify the lowest-energy state of
the reduced-order model. Ctrl-click to specify multiple values.
When you change the specified reduced model order, Model Reducer automatically
computes a new reduced-order model. If you specify multiple model orders, Model
Reducer computes multiple reduced-order models and displays their responses on
the plot.
6-20
Balanced Truncation Model Reduction
4 Optionally, examine the absolute or relative error between the original and reduced-
order model, in addition to the frequency response. Select the error-plot type using
the buttons on the Balanced Truncation tab.
6-21
6 Model Simplification
For more information about using the analysis plots, see “Visualize Reduced-Order
Models in the Model Reducer App” on page 6-67.
5 If low-frequency dynamics are not important to your application, you can clear the
Preserve DC Gain checkbox. Doing so sometimes yields a better match at higher
frequencies between the original and reduced-order models.
6-22
Balanced Truncation Model Reduction
When you check or clear the Preserve DC Gain checkbox, Model Reducer
automatically computes new reduced-order models. For more information about this
option, see “Compare Truncated and DC Matched Low-Order Model Approximations”
on page 6-30.
6 Optionally, limit the Hankel singular-value computation to a specific frequency range.
Such a limit is useful when the model has modes outside the region of interest to your
particular application. When you apply a frequency limit, Model Reducer
determines which states to truncate based on their energy contribution within the
specified frequency range only. Neglecting energy contributions outside that range
can yield an even lower-order approximation that is still adequate for your
application.
6-23
6 Model Simplification
• In the text box, entering a vector of the form [fmin,fmax]. Units are rad/
TimeUnit, where TimeUnit is the TimeUnit property of the model you are
reducing.
• On the response plot or error plot, dragging the boundaries of the shaded region
or the shaded region itself. Model Reducer analyzes the state contributions
within the shaded region only.
When you check or clear the Select frequency range checkbox or change the
selected range, Model Reducer automatically computes new reduced-order models.
further, click . The new models appear in the Data Browser. If you have specified
multiple orders, each reduced model appears separately. Model names reflect the
reduced model order.
6-24
Balanced Truncation Model Reduction
After creating reduced models in the Data Browser, you can continue changing the
reduction parameters and create reduced models with different orders for analysis
and comparison.
You can now perform further analysis with the reduced model. For example:
• Examine other responses of the reduced system, such as the step response or Nichols
plot. To do so, use the tools on the Plots tab. See “Visualize Reduced-Order Models in
the Model Reducer App” on page 6-67 for more information.
• Export reduced models to the MATLAB workspace for further analysis or control
To create a MATLAB script you can use for further model-reduction tasks at the command
line, click Create Reduced Model, and select Generate MATLAB Script.
6-25
6 Model Simplification
Model Reducer creates a script that uses the balred command to perform model
reduction with the parameters and options you have set on the Balanced Truncation
tab. The script opens in the MATLAB editor.
To do so, first examine the contribution of the various states to the overall model
behavior. Choose the approximation order based on the number of states that make a
significant contribution to the overall model behavior.
For this example, load a high-order model. hplant is a 23rd-order SISO model.
6-26
Balanced Truncation Model Reduction
ans = 23
Examine the relative amount of energy per state in hplant using a Hankel singular-value
(HSV) plot.
hsvplot(hplant)
Small Hankel singular values indicate that the associated states contribute little to the
behavior of the system. The plot shows that two states account for most of the energy in
the system. Therefore, try simplifying the model to just first or second order.
6-27
6 Model Simplification
opts = balredOptions('StateElimMethod','Truncate');
hplant1 = balred(hplant,1,opts);
hplant2 = balred(hplant,2,opts);
The second argument to balred specifies the target approximation order, so that
hplant1 is a first-order approximation and hplant2 is a second-order approximation of
hplant. By default, balred discards the states with the smallest Hankel singular values,
and alters the remaining states to preserve the DC gain of the system. Setting the
StateElimMethod option to Truncate causes balred to discard low-energy states
without altering the remaining states.
When working with reduced-order models, it is important to verify that the approximation
does not introduce inaccuracies at frequencies that are important for your application.
Therefore, compare the frequency responses of the original and approximated systems.
For MIMO systems, use the sigmaplot command. For this SISO system, examine a Bode
plot.
bodeplot(hplant,hplant2,hplant1)
legend('Original','2nd order','1st order')
6-28
Balanced Truncation Model Reduction
The second-order approximation hplant2 matches the original 23rd-order system very
well, especially at lower frequencies. The first-order system does not match as well.
In general, as you decrease the order of the approximated model, the frequency response
of the approximated model begins to differ from the original model. Choose an
approximation that is sufficiently accurate in the bands that are important to you. For
example, in a control system you might want good accuracy inside the control bandwidth.
Accuracy at frequencies far above the control bandwidth, where the gain rapidly rolls off,
might be less important.
You can also validate the approximation in the time domain. For instance, examine the
step responses of the original and reduced-order systems.
6-29
6 Model Simplification
stepplot(hplant,hplant2,'r--',hplant1,'g--')
legend('Original','2nd order','1st order','Location','SouthEast')
This result confirms that the second-order approximation is a good match to the original
23rd-order system.
6-30
Balanced Truncation Model Reduction
• Discard the states that make the smallest contribution to system behavior, altering the
remaining states to preserve the DC gain of the system.
• Discard the low-energy states without altering the remaining states.
Which method you choose depends on what dynamics are most important to your
application. In general, preserving DC gain comes at the expense of accuracy in higher-
frequency dynamics. Conversely, state truncation can yield more accuracy in fast
transients, at the expense of low-frequency accuracy.
6-31
6 Model Simplification
T =
2 s (s+2)
--------------------------------
(s+0.004277) (s+1.588) (s+4.418)
Compute two second-order approximations to T, one that preserves the DC gain and one
that truncates the lowest-energy state without changing the other states. Use
balredOptions to specify the approximation methods, MatchDC and Truncate,
respectively.
matchopt = balredOptions('StateElimMethod','MatchDC');
truncopt = balredOptions('StateElimMethod','Truncate');
Tmatch = balred(T,2,matchopt);
Ttrunc = balred(T,2,truncopt);
bodeplot(T,Tmatch,Ttrunc)
legend('Original','DC Match','Truncate')
6-32
Balanced Truncation Model Reduction
The truncated model Ttrunc matches the original model well at high frequencies, but
differs considerably at low frequency. Conversely, Tmatch yields a good match at low
frequencies as expected, at the expense of high-frequency accuracy.
You can also see the differences between the two methods by examining the time-domain
response in different regimes. Compare the slow dynamics by looking at the step
response of all three models with a long time horizon.
stepplot(T,Tmatch,'r--',Ttrunc,1500)
legend('Original','DC Match','Truncate')
6-33
6 Model Simplification
As expected, on long time scales the DC-matched approximation Tmatch has a very
similar response to the original model.
stepplot(T,Tmatch,'r',Ttrunc,'g--',0.5)
legend('Original','DC Match','Truncate')
6-34
Balanced Truncation Model Reduction
On short time scales, the truncated approximation Ttrunc provides a better match to the
original model. Which approximation method you should use depends on which regime is
most important for your application.
When computing a reduced-order approximation, the balred command (or the Model
Reducer app) does not eliminate unstable poles because doing so would fundamentally
change the system dynamics. Instead, the software decomposes the model into stable and
unstable parts and reduces the stable part of the model.
6-35
6 Model Simplification
If your model has near-unstable poles, you might want to ensure that the reduced-order
approximation preserves these dynamics. This example shows how to use the Offset
option of balred to preserve poles that are close to the stable-unstable boundary. You
can achieve the same result in the Model Reducer app, on the Balanced Truncation
tab, under Options, using the Offset field, as shown:
load('reduce.mat','gasf35unst')
gasf35unst is a 25-state SISO model with two unstable poles (Re(s) > 0). Examine the
system poles to find the near-unstable poles.
pzplot(gasf35unst)
axis([-0.0015 0.0015 -0.0005 0.0005])
6-36
Balanced Truncation Model Reduction
The pole-zero plot shows several poles (marked by x) that fall in the left half-plane, but
relatively close to the imaginary axis. These are the near-unstable poles. Two of these fall
within 0.0005 of instability. Three more fall within 0.001 of instability.
hsvplot(gasf35unst)
6-37
6 Model Simplification
The plot shows the two unstable modes, but you cannot easily determine the energy
contribution of the near-unstable poles. In your application, you might want to reduce the
model without discarding those poles nearest to instability, even if they are of relatively
low energy. Use the Offset option of balred to calculate a reduced-order system that
preserves the two stable poles that are closest to the imaginary axis. The Offset option
sets the boundary between poles that balred can discard, and poles that balred must
preserve (treat as unstable).
opts = balredOptions('Offset',0.0005);
gasf_arr = balred(gasf35unst,[10 15],opts);
Providing balred an array of target approximation orders [10 15] causes balred to
return an array of approximated models. The array gasf_arr contains two models, a
6-38
Balanced Truncation Model Reduction
bodeplot(gasf35unst,gasf_arr,'r--')
The 15th order approximation is a good frequency-domain match to the original model.
However, the 10th-order approximation shows changes in high-frequency dynamics,
6-39
6 Model Simplification
which might be too large to be acceptable. The 15th-order approximation is likely a better
choice.
load(fullfile(matlabroot,'examples','control','build.mat'),'G')
bodeplot(G)
6-40
Balanced Truncation Model Reduction
G is a 48th-order model with several large peak regions around 5.2 rad/s, 13.5 rad/s, and
24.5 rad/s, and smaller peaks scattered across many frequencies. Examine the Hankel
singular-value plot to see the energy contributions of the model's 48 states.
hsvd(G)
6-41
6 Model Simplification
The singular-value plot suggests that you can discard at least 20 states without significant
impact on the overall system response. Suppose that for your application you are only
interested in the dynamics near the second large peak, between 10 rad/s and 22 rad/s.
Try a few reduced-model orders based on the Hankel singular value plot. Compare their
frequency responses to the original model, especially in the region of that peak.
G18 = balred(G,18);
G10 = balred(G,10);
bodeplot(G,G18,G10,logspace(0.5,1.5,100));
legend('Original','Order 18','Order 10');
6-42
Balanced Truncation Model Reduction
The 18th-order model is a good match to the dynamics in the region of interest. In the
10th order model, however, there is some degradation of the match.
Focus the model reduction on the region of interest to obtain a good match with a lower-
order approximation. First, examine the state energy contributions in that frequency
region only. Use hsvdOptions to specify the frequency interval for hsvd.
hopt = hsvdOptions('FreqIntervals',[10,22]);
hsvd(G,hopt)
6-43
6 Model Simplification
Comparing this plot to the previous Hankel singular-value plot shows that in this
frequency region, many fewer states contribute significantly to the dynamics than
contribute to the overall dynamics.
Try the same reduced-model orders again, this time choosing states to eliminate based
only on their contribution to the frequency interval. Use balredOptions to specify the
frequency interval for balred.
bopt = balredOptions('StateElimMethod','Truncate','FreqIntervals',[10,22]);
GLim18 = balred(G,18,bopt);
GLim10 = balred(G,10,bopt);
bodeplot(G,GLim18,GLim10,logspace(0.5,1.5,100));
legend('Original','Order 18','Order 10');
6-44
See Also
See Also
Apps
Model Reducer
Functions
balred | hsvplot
6-45
6 Model Simplification
Related Examples
• “Mode-Selection Model Reduction” on page 6-57
• “Pole-Zero Simplification” on page 6-47
• “Model Reduction Basics” on page 6-2
6-46
Pole-Zero Simplification
Pole-Zero Simplification
Pole-zero simplification reduces the order of your model exactly by canceling pole-zero
pairs or eliminating states that have no effect on the overall model response. Pole-zero
pairs can be introduced, for example, when you construct closed-loop architectures.
Normal small errors associated with numerical computation can convert such canceling
pairs to near-canceling pairs. Removing these states preserves the model response
characteristics while simplifying analysis and control design. Types of pole-zero
simplification include:
• Structural elimination — Eliminate states that are structurally disconnected from the
inputs or outputs. Eliminating structurally disconnected states is a good first step in
model reduction because the process does not involve any numerical computation. It
also preserves the state structure of the remaining states. At the command line,
perform structural elimination with sminreal.
• Pole-zero cancellation or minimal realization — Eliminate canceling or near-canceling
pole-zero pairs from transfer functions. Eliminate unobservable or uncontrollable
states from state-space models. At the command line, perform this kind of
simplification with minreal.
1 Open the app and import a model to reduce. For instance, suppose that there is a
model named build in the MATLAB workspace. The following command opens
Model Reducer and imports the LTI model build.
modelReducer(build)
2
In the Data Browser, select the model to reduce. Click Pole-Zero
Simplification.
6-47
6 Model Simplification
6-48
Pole-Zero Simplification
The pole-zero map marks pole locations with x and zero locations with o.
Note The frequency response is a Bode plot for SISO models, and a singular-value
plot for MIMO models.
3 Optionally, change the tolerance with which Model Reducer identifies canceling
pole-zero pairs. Model Reducer cancels pole-zero pairs that fall within the tolerance
specified by the Simplification of pole-zero pairs value. In this case, no pole-zero
6-49
6 Model Simplification
pairs are close enough together for Model Reducer to cancel them at the default
tolerance of 1e-05. To cancel pairs that are a little further apart, move the slider to
the right or enter a larger value in the text box.
The blue x and o marks on the pole-zero map show the near-canceling pole-zero pairs
in the original model that are eliminated from the simplified model. Poles and zeros
remaining in the simplified model are marked with red x and o.
4 Try different simplification tolerances while observing the frequency response of the
original and simplified model. Remove as many poles and zeros as you can while
preserving the system behavior in the frequency region that is important for your
6-50
Pole-Zero Simplification
application. Optionally, examine absolute or relative error between the original and
simplified model. Select the error-plot type using the buttons on the Pole-Zero
Simplification tab.
For more information about using the analysis plots, see “Visualize Reduced-Order
Models in the Model Reducer App” on page 6-67.
5 When you have a simplified model that you want to store and analyze further, click
. The new model appears in the Data Browser with a name that reflects the
reduced model order.
6-51
6 Model Simplification
After creating a reduced model in the Data Browser, you can continue changing the
simplification parameters and create reduced models with different orders for
analysis and comparison.
You can now perform further analysis with the reduced model. For example:
• Examine other responses of the reduced system, such as the step response or Nichols
plot. To do so, use the tools on the Plots tab. See “Visualize Reduced-Order Models in
the Model Reducer App” on page 6-67 for more information.
• Export reduced models to the MATLAB workspace for further analysis or control
To create a MATLAB script you can use for further model-reduction tasks at the command
line, click Create Reduced Model, and select Generate MATLAB Script.
6-52
Pole-Zero Simplification
Model Reducer creates a script that uses the minreal command to perform model
reduction with the parameters you have set on the Pole-Zero Simplification tab. The
script opens in the MATLAB editor.
Create a model of the following system, where C is a PI controller, and G has a zero at
rad/s. Such a low-frequency zero can arise from derivative action somewhere in
the plant dynamics. For example, the plant may include a component that computes speed
from position measurements.
6-53
6 Model Simplification
G = zpk(3e-8,[-1,-3],1);
C = pid(1,0.3);
T = feedback(G*C,1)
T =
(s+0.3) (s-3e-08)
----------------------
s (s+4.218) (s+0.7824)
In the closed-loop model T, the integrator from C very nearly cancels the low-
frequency zero of G.
Force a cancelation of the integrator with the zero near the origin.
Tred = minreal(T,1e-7)
Tred =
(s+0.3)
--------------------
(s+4.218) (s+0.7824)
By default, minreal reduces transfer function order by canceling exact pole-zero pairs or
near pole-zero pairs within sqrt(eps). Specifying 1e-7 as the second input causes
minreal to eliminate pole-zero pairs within rad/s of each other.
The reduced model Tred includes all the dynamics of the original closed-loop model T,
except for the near-canceling zero-pole pair.
6-54
Pole-Zero Simplification
bode(T,Tred,'r--')
legend('T','Tred')
Because the canceled pole and zero do not match exactly, some extreme low-frequency
dynamics evident in the original model are missing from Tred. In many applications, you
can neglect such extreme low-frequency dynamics. When you increase the matching
tolerance of minreal, make sure that you do not eliminate dynamic features that are
relevant to your application.
6-55
6 Model Simplification
See Also
Apps
Model Reducer
Functions
minreal | sminreal
Related Examples
• “Balanced Truncation Model Reduction” on page 6-17
• “Mode-Selection Model Reduction” on page 6-57
• “Model Reduction Basics” on page 6-2
6-56
Mode-Selection Model Reduction
For more general information about model reduction, see “Model Reduction Basics” on
page 6-2.
1 Open the app and import an LTI model to reduce. For instance, suppose that there is
a model named Gms in the MATLAB workspace. The following command opens Model
Reducer and imports the model.
modelReducer(Gms)
2
In the Data Browser, select the model to reduce. Click Mode Selection.
6-57
6 Model Simplification
In the Mode Selection tab, Model Reducer displays a plot of the frequency
response of the original model and a reduced version of the model. The app also
displays a pole-zero map of both models.
6-58
Mode-Selection Model Reduction
The pole-zero map marks pole locations with x and zero locations with o.
Note The frequency response is a Bode plot for SISO models, and a singular-value
plot for MIMO models.
6-59
6 Model Simplification
3 Model Reducer eliminates poles that lie outside the shaded region. Change the
shaded region to capture only the dynamics you want to preserve in the reduced
model. There are two ways to do so.
• On either the response plot or the pole-zero map, drag the boundaries of the
shaded region or the shaded region itself.
• On the Mode Selection tab, enter lower and upper cutoff frequencies.
When you change the shaded regions or cutoff frequencies, Model Reducer
automatically computes a new reduced-order model. All poles retained in the reduced
6-60
Mode-Selection Model Reduction
model fall within the shaded region on the pole-zero map. The reduced model might
contain zeros that fall outside the shaded region.
4 Optionally, examine absolute or relative error between the original and simplified
model. Select the error-plot type using the buttons on the Mode Selection tab.
For more information about using the analysis plots, see “Visualize Reduced-Order
Models in the Model Reducer App” on page 6-67.
5 When you have one or more reduced models that you want to store and analyze
6-61
6 Model Simplification
After creating a reduced model in the Data Browser, you can continue adjusting the
mode-selection region to create reduced models with different orders for analysis and
comparison.
You can now perform further analysis with the reduced model. For example:
• Examine other responses of the reduced system, such as the step response or Nichols
plot. To do so, use the tools on the Plots tab. See “Visualize Reduced-Order Models in
the Model Reducer App” on page 6-67 for more information.
• Export reduced models to the MATLAB workspace for further analysis or control
To create a MATLAB script you can use for further model-reduction tasks at the command
line, click Create Reduced Model, and select Generate MATLAB Script.
6-62
Mode-Selection Model Reduction
Model Reducer creates a script that uses the freqsep command to perform model
reduction with the parameters you have set on the Mode Selection tab. The script opens
in the MATLAB editor.
For this example, load the model Gms and examine its frequency response.
6-63
6 Model Simplification
Gms has two sets of resonances, one at relatively low frequency and the other at relatively
high frequency. Suppose that you want to tune a controller for Gms, but the actuator in
your system is limited to a bandwidth of about 3 rad/s, in between the two groups of
resonances. To simplify calculation and tuning using Gms, you can use mode selection to
eliminate the high-frequency dynamics.
[Gms_s,Gms_f] = freqsep(Gms,30);
freqsep decomposes Gms into slow and fast components such that Gms = Gms_s +
Gms_f. All modes (poles) with natural frequency less than 30 are in Gms_s, and the
higher-frequency poles are in Gms_f.
bodeplot(Gms,Gms_s,Gms_f)
legend('original','slow','fast')
6-64
Mode-Selection Model Reduction
The slow component, Gms_s, contains only the lower-frequency resonances and matches
the DC gain of the original model. Examine the orders of both models.
order(Gms)
ans = 18
order(Gms_s)
ans = 10
When the high-frequency dynamics are unimportant for your application, you can use the
10th-order Gms_s instead of the original 18th-order model. If neglecting low-frequency
dynamics is appropriate for your application, you can use Gms_f. To select modes that fall
between a low-frequency and a high-frequency cutoff, use additional calls to freqsep.
6-65
6 Model Simplification
See Also
Model Reducer | freqsep
Related Examples
• “Balanced Truncation Model Reduction” on page 6-17
• “Pole-Zero Simplification” on page 6-47
• “Model Reduction Basics” on page 6-2
6-66
Visualize Reduced-Order Models in the Model Reducer App
For more general information about model reduction, see “Model Reduction Basics” on
page 6-2.
Error Plots
By default, for any model reduction method, Model Reducer shows a frequency-response
plot of both the original and reduced models. This plot is a Bode plot for SISO models,
and a singular-value plot for MIMO models.
To more closely examine the differences between an original model and a reduced model,
you can use absolute error or relative error plots. On any model reduction tab, click
Absolute error plot or Relative error plot to view these plots.
6-67
6 Model Simplification
• Absolute error plot — Shows the singular values of G-Gr, where G is the original
model and Gr is the current reduced model.
• Relative error plot — Shows the singular values of (G-Gr)/G. This plot is useful
when the model has very high or very low gain in the region that is important to your
application. In such regions, absolute error can be misleading.
For SISO models, the singular-value plot is the magnitude of the frequency response.
Response Plots
After you click to add one or more reduced models to the Data Browser, compare
additional responses of the original and reduced models using the Plots tab.
In the Data Browser, select one or more models to plot. (Ctrl-click to select multiple
models.) Then, on the Plots tab, click the type of plot you want to create.
6-68
Visualize Reduced-Order Models in the Model Reducer App
In the Data Browser, select the model to add. Then, on the Plots tab, click the icon
corresponding to the plot you want to update. Plots you have created appear on the left
side of the plot gallery.
6-69
6 Model Simplification
Plot Characteristics
On any plot in Model Reducer:
• To see response information and data values, click a line on the plot.
6-70
Visualize Reduced-Order Models in the Model Reducer App
6-71
6 Model Simplification
Plot Tools
Mouse over any plot to access plot tools at the upper right corner of the plot.
6-72
Visualize Reduced-Order Models in the Model Reducer App
•
and — Zoom in and zoom out. Click to activate, and drag the cursor over the
region to zoom. The zoom icon turns dark when zoom is active. Right-click while zoom
is active to access additional zoom options.Click the icon again to deactivate.
•
— Pan. Click to activate, and drag the cursor across the plot area to pan. The pan
icon turns dark when pan is active. Right-click while pan is active to access additional
pan options. Click the icon again to deactivate.
•
— Legend. By default, the plot legend is active. To toggle the legend off and on,
click this icon. To move the legend, drag it to a new location on the plot.
To change the way plots are tiled or sorted, use the options on the View tab.
6-73
6 Model Simplification
See Also
Model Reducer
Related Examples
• “Balanced Truncation Model Reduction” on page 6-17
• “Mode-Selection Model Reduction” on page 6-57
• “Pole-Zero Simplification” on page 6-47
6-74
Linear Analysis
75
7
Time Responses
sys =
8 s^2 + 18 s + 32
-----------------------
s^3 + 6 s^2 + 14 s + 24
You can plot the step and impulse responses of this system using the step and impulse
commands:
subplot(2,1,1)
step(sys)
subplot(2,1,2)
impulse(sys)
7-2
Plotting System Responses
You can also simulate the response to an arbitrary signal, for example, a sine wave, using
the lsim command. The input signal appears in gray and the system's response in blue.
clf
t = 0:0.01:4;
u = sin(10*t);
lsim(sys,u,t) % u,t define the input signal
7-3
7 Time Domain Analysis
You can use the plotting commands with continuous or discrete TF, SS, or ZPK models.
For state-space models, you can also plot the unforced response from some given initial
state, for example:
7-4
Plotting System Responses
Frequency Responses
sys =
8 s^2 + 18 s + 32
-----------------------
7-5
7 Time Domain Analysis
s^3 + 6 s^2 + 14 s + 24
bode(sys)
grid
nyquist(sys)
grid
7-6
Plotting System Responses
nichols(sys)
grid
7-7
7 Time Domain Analysis
The poles and zeros of a system contain valuable information about its dynamics, stability,
and limits of performance. For example, consider the feedback loop in Figure 1 where
7-8
Plotting System Responses
For the gain value k = 0.7, you can plot the closed-loop poles and zeros using pzmap:
s = tf('s');
G = -(2*s+1)/(s^2+3*s+2);
k = 0.7;
T = feedback(G*k,1);
pzmap(T)
grid, axis([-2 0 -1 1])
7-9
7 Time Domain Analysis
The closed-loop poles (marked by blue x's) lie in the left half-plane so the feedback loop is
stable for this choice of gain k. You can read the damping ratio of the closed-loop poles
from this chart (see labels on the radial lines). Here the damping ratio is about 0.7,
suggesting a well-damped closed-loop response as confirmed by:
clf
step(T)
7-10
Plotting System Responses
To further understand how the loop gain k affects closed-loop stability, you can plot the
locus of the closed-loop poles as a function of k:
rlocus(G)
grid
7-11
7 Time Domain Analysis
Clicking where the locus intersects the y axis reveals that the closed-loop poles become
unstable for k = 1.51. So the loop gain should remain smaller than 1.5 for closed-loop
stability.
7-12
Plotting System Responses
Response Characteristics
Using the example from the previous section, plot the closed-loop step response:
step(T)
7-13
7 Time Domain Analysis
Now, right-click on the plot to display the Peak Response and Settling Time
Characteristics, and click on the blue dots to read the corresponding overshoot and
settling time values:
7-14
Plotting System Responses
All commands mentioned so far fully support multi-input multi-output (MIMO) systems. In
the MIMO case, these commands produce arrays of plots. For example, the step response
of the two-input, two-output system
sys = rss(3,2,2);
sys.A = [-0.5 -0.3 -0.2 ; 0 -1.3 -1.7; 0.4 1.7 -1.3];
is a 2-by-2 array of plots where each column shows the step response of a particular input
channel:
step(sys)
7-15
7 Time Domain Analysis
If desired, you can group all four responses on a single plot by right-clicking on the plot
and selecting the I/O Grouping -> All submenu. The resulting plot is shown below.
7-16
Plotting System Responses
The following additional plots are useful for analyzing MIMO systems:
• Singular value plot (sigma), which shows the principal gains of the frequency
response
• Pole/zero map for each I/O pair (iopzplot)
sigma(sys)
grid
7-17
7 Time Domain Analysis
Comparing Systems
You can plot multiple systems at once using any of the response plot commands. You can
assign a specific color, marker, or line style to each system for easy comparison. Using the
feedback example above, plot the closed-loop step response for three values of the loop
gain k in three different colors:
k1 = 0.4;
T1 = feedback(G*k1,1);
k2 = 1;
T2 = feedback(G*k2,1);
step(T,'b',T1,'r',T2,'g')
legend('k = 0.7','k = 0.4','k = 1')
7-18
Plotting System Responses
7-19
7 Time Domain Analysis
Time-Domain Responses
When you perform time-domain analysis of a dynamic system model, you may want one or
more of the following:
Control System Toolbox time-domain analysis commands can obtain these results for any
kind of dynamic system model (for example, continuous or discrete, SISO or MIMO, or
arrays of models) except for frequency response data models.
See Also
Related Examples
• “Time-Domain Response Data and Plots” on page 7-21
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
7-20
Time-Domain Response Data and Plots
Create a transfer function model and plot its response to a step input at = 0.
When call step without output arguments, it plots the step response on the screen.
Unless you specify a time range to plot, step automatically chooses a time range that
illustrates the system dynamics.
7-21
7 Time Domain Analysis
Calculate the step response data from = 0 (application of the step input) to = 8 s.
[y,t] = step(H,8);
When you call step with output arguments, the command returns the step response data
y. The vector t contains corresponding time values.
Plot the response of H to an impulse input applied at = 0. Plot the response with a grid.
opts = timeoptions;
opts.Grid = 'on';
impulseplot(H,opts)
Use the timeoptions command to define options sets for customizing time-domain plots
with commands like impulseplot and stepplot.
7-22
See Also
Calculate 200 points of impulse response data from = 1 (one second after application of
the impulse input) to = 3s.
[y,t] = impulse(H,linspace(1,3,200));
As for step, you can omit the time vector to allow impulse to automatically select a time
range.
See Also
impulse | impulseplot | step | stepplot | timeoptions
Related Examples
• “Time-Domain Characteristics on Response Plots” on page 7-24
• “Time-Domain Responses of Multiple Models” on page 7-36
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
More About
• “Time-Domain Responses” on page 7-20
7-23
7 Time Domain Analysis
You can use similar procedures to display system characteristics on impulse response
plots or initial value response plots, such as peak response or settling time.
Create a transfer function model and plot its response to a step input at t = 0.
H = tf([8 18 32],[1 6 14 24]);
stepplot(H)
Right-click anywhere in the figure and select Characteristics > Peak Response from
the menu.
7-24
Time-Domain Characteristics on Response Plots
A marker appears on the plot indicating the peak response. Horizontal and vertical dotted
lines indicate the time and amplitude of that response.
7-25
7 Time Domain Analysis
Click the marker to view the value of the peak response and the overshoot in a datatip.
7-26
See Also
You can use a similar procedure to select other characteristics such as settling time and
rise time from the Characteristics menu and view the values.
See Also
impulse | lsiminfo | step | stepinfo
Related Examples
• “Numeric Values of Time-Domain System Characteristics” on page 7-29
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
7-27
7 Time Domain Analysis
More About
• “Time-Domain Responses” on page 7-20
7-28
Numeric Values of Time-Domain System Characteristics
Create a dynamic system model and get numeric values of the system’s step response
characteristics.
The output is a structure that contains values for several step response characteristics. To
access these values or refer to them in other calculations, use dot notation. For example,
data.Overshoot is the overshoot value.
Calculate the time it takes the step response of H to settle within 0.5% of its final value.
data = stepinfo(H,'SettlingTimeThreshold',0.005);
t05 = data.SettlingTime
t05 = 4.8896
By default, stepinfo defines the settling time as the time it takes for the output to settle
within 0.02 (2%) of its final value. Specifying a more stringent
'SettlingTimeThreshold' of 0.005 results in a longer settling time.
For more information about the options and the characteristics, see the stepinfo
reference page.
7-29
7 Time Domain Analysis
See Also
lsiminfo | stepinfo
Related Examples
• “Time-Domain Characteristics on Response Plots” on page 7-24
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
More About
• “Time-Domain Responses” on page 7-20
7-30
Time-Domain Responses of Discrete-Time Model
You can use the techniques of this example with commands such as impulse, initial,
impulseplot, and initialpot to obtain time-domain responses of discrete-time
models.
Create a discrete-time transfer function model and plot its response to a step input at =
0.
H = tf([-0.06,0.4],[1,-1.6,0.78],0.1);
step(H)
7-31
7 Time Domain Analysis
For discrete-time models, step plots the response at multiples of the sample time,
assuming a hold between samples.
[y,t] = step(H,0.5:0.1:2.5);
When you specify a time vector for the response of a discrete-time model, the time step
must match the sample time Ts of the discrete-time model. The vector t contains the time
points between 0.5 and 2.5 seconds, at multiples of the sample time of H, 0.1 s. The vector
y contains the corresponding step response values.
7-32
See Also
See Also
impulse | impulseplot | initial | initialplot | step | stepplot
Related Examples
• “Time-Domain Responses of MIMO Model” on page 7-34
• “Time-Domain Responses of Multiple Models” on page 7-36
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
More About
• “Time-Domain Responses” on page 7-20
7-33
7 Time Domain Analysis
Create a MIMO model and plot its response to a t = 0 impulse at all inputs.
H = rss(2,2,2);
H.InputName = 'Control';
H.OutputName = 'Temperature';
impulse(H)
7-34
See Also
impulse plots the response of each output to an impulse applied at each input. (Because
rss generates a random state-space model, you might see different responses from those
pictured.) The first column of plots shows the response of each output to an impulse
applied at the first input, Control(1). The second column shows the response of each
output to an impulse applied at the second input, Control(2).
Calculate the impulse responses of all channels of H, and examine the size of the output.
[y,t] = impulse(H);
size(y)
ans = 1×3
207 2 2
The first dimension of the data array y is the number of samples in the time vector t. The
impulse command determines this number automatically if you do not supply a time
vector. The remaining dimensions of y are the numbers of outputs and inputs in H. Thus,
y(:,i,j) is the response at the i th output of H to an impulse applied at the j th input.
See Also
impulse | impulseplot | initial | initialplot | step | stepplot
Related Examples
• “Time-Domain Responses of Multiple Models” on page 7-36
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
More About
• “Time-Domain Responses” on page 7-20
7-35
7 Time Domain Analysis
For this example, obtain two models whose time responses you want to compare, and plot
them on a single step plot. For instance, you can compare a third-order plant G, and the
closed-loop response of G with a controller C1 having integral action.
G = zpk([],[-5 -5 -10],100);
C1 = pid(0,4.4);
CL1 = feedback(G*C1,1);
step(G,CL1);
7-36
Time-Domain Responses of Multiple Models
When you provide multiple models to step as input arguments, the command displays the
responses of both models on the same plot. If you do not specify a time range to plot,
step attempts to choose a time range that illustrates the dynamics of all the models.
Compare the step response of the closed-loop model with another controller. Specify plot
colors and styles for each response.
C2 = pid(2.9,7.1);
CL2 = feedback(G*C2,1);
step(G,'b--',CL1,'g-',CL2,'r-')
7-37
7 Time Domain Analysis
You can specify custom plot color and style for each response in the plot. For example,
'g-' specifies a solid green line for response CL2. For additional plot customization
options, use stepplot.
See Also
Linear System Analyzer | impulse | impulseplot | initial | initialplot | step |
stepplot
7-38
See Also
Related Examples
• “Time-Domain Responses of MIMO Model” on page 7-34
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
More About
• “Time-Domain Responses” on page 7-20
7-39
7 Time Domain Analysis
For example, compare a third-order plant G, and the closed-loop responses of G with two
different controllers, C1 and C2.
G = zpk([],[-5 -5 -10],100);
C1 = pid(0,4.4);
T1 = feedback(G*C1,1);
C2 = pid(2.9,7.1);
T2 = feedback(G*C2,1);
Open the Linear System Analyzer tool to examine the responses of the plant and the
closed-loop systems.
linearSystemAnalyzer(G,T1,T2)
7-40
Joint Time-Domain and Frequency-Domain Analysis
By default, the Linear System Analyzer launches with a plot of the step response of the
three systems. Click to add a legend to the plot.
Add plots of the impulse responses to the Linear System Analyzer display.
In the Linear System Analyzer, select Edit > Plot Configurations to open the Plot
Configurations dialog box.
7-41
7 Time Domain Analysis
Select the two-plot configuration. In the Response Type area, select Bode Magnitude for
the second plot type.
Click OK to add the Bode plots to the Linear System Analyzer display.
Right-click anywhere in the Bode Magnitude plot and select Characteristics > Peak
Response from the menu.
7-42
Joint Time-Domain and Frequency-Domain Analysis
Markers appear on the plot indicating the peak response values. Horizontal and vertical
dotted lines indicate the frequency and amplitude of those responses. Click on a marker
to view the value of the peak response in a datatip.
You can use a similar procedure to select other characteristics such as settling time and
rise time from the Characteristics menu and view the values.
You can also change the type of plot displayed in the Linear System Analyzer. For
example, to change the first plot type to a plot of the impulse response, right-click
anywhere in the plot. Select Plot Types > Impulse
7-43
7 Time Domain Analysis
The displayed plot changes to show the impulse of the three systems.
See Also
Linear System Analyzer | impulse | impulseplot | initial | initialplot | step |
stepplot
Related Examples
• “Time-Domain Responses of Multiple Models” on page 7-36
More About
• “Time-Domain Responses” on page 7-20
7-44
Response from Initial Conditions
sys_dc =
A =
Current w
Current -4 -0.03
w 0.75 -10
B =
Volts
Current 2
w 0
C =
Current w
w 0 1
D =
Volts
w 0
This example uses the SISO, 2-state model sys_dc. This model represents a DC motor.
The input is an applied voltage, and the output is the angular rate of the motor ω. The
states of the model are the induced current (x1), and ω (x2). The model display in the
command window shows the labeled input, output, and states.
Plot the undriven evolution of the motor's angular rate from an initial state in which the
induced current is 1.0 amp and the initial rotation rate is 5.0 rad/s.
7-45
7 Time Domain Analysis
x0 = [1.0 5.0];
initial(sys_dc,x0)
initial plots the time evolution from the specified initial condition on the screen.
Unless you specify a time range to plot, initial automatically chooses a time range that
illustrates the system dynamics.
Calculate the time evolution of the output and the states of sys_dc from =0
(application of the step input) to = 1 s.
t = 0:0.01:1;
[y,t,x] = initial(sys_dc,x0,t);
7-46
See Also
The vector y contains the output at each time step in t. The array x contains the state
values at each time step. Therefore, in this example x is a 2-by-101 array. Each row of x
contains the values of the two states of sys_dc at the corresponding time step.
See Also
impulse | initial | initialplot | step
Related Examples
• “Time-Domain Response Data and Plots” on page 7-21
• “Numeric Values of Time-Domain System Characteristics” on page 7-29
More About
• “Time-Domain Responses” on page 7-20
7-47
7 Time Domain Analysis
In the block parameters, set the LTI system variable parameter to the LTI model to
import. For state-space models, set the Initial states parameter to a vector to specify
non-zero initial states.
To specify a model for the LTI System block, set the LTI system variable block
parameter to either:
• The variable name of an LTI model in the MATLAB® workspace or model workspace,
such as sys.
• A MATLAB expression that evaluates to an LTI model, such as tf(1,[1 1]).
For example, you can specify a state-space (ss), zero-pole-gain (zpk), or transfer function
(tf) model. You can simulate SISO models or MIMO models, and continuous-time or
discrete-time models.
7-48
Import LTI Model Objects into Simulink
7-49
7 Time Domain Analysis
This example simulates the system response to a step input at t = 2 s. Use the LTI System
block to import an LTI model object anywhere in your Simulink model to simulate the
linear system response to any input.
The LTI System block has one input and one output, even when you specify a MIMO
model for the block. In that case, the block input and output become vector signals. For
instance, the model LTISystemBlockMIMO uses an LTI system block to represent a
MIMO plant in a control system.
In this model, the LTI System specified in the block is Gm, a 2-output, 2-input transfer
function model stored in the model workspace. A Mux block combines the two controller
outputs into a vector signal for the LTI System block input. Similarly, a Demux block
separates the vector output of the LTI System block into two scalar signals.
7-50
Import LTI Model Objects into Simulink
7-51
7 Time Domain Analysis
This example simulates a closed-loop system response to a t = 50 s step at the first input
and a t = 150 s step at the second input. You can use the LTI system block anywhere you
want to insert an LTI system into a Simulink model.
See Also
LTI System
7-52
Analysis of Systems with Time Delays
For example, consider the following control loop, where the plant is modeled as first-order
plus dead time:
You can model the closed-loop system from r to y with the following commands:
s = tf('s');
P = 5*exp(-3.4*s)/(s+1);
C = 0.1 * (1 + 1/(5*s));
T = feedback(P*C,1);
T is a state-space model with an internal delay. For more information about models with
internal delays, see “Closing Feedback Loops with Time Delays” on page 2-45.
stepplot(T)
7-53
7 Time Domain Analysis
For more complicated interconnections, you can name the input and output signals of
each block and use connect to automatically take care of the wiring. Suppose, for
example, that you want to add feedforward to the control loop of the previous model.
7-54
Analysis of Systems with Time Delays
F = 0.3/(s+4);
P.InputName = 'u';
P.OutputName = 'y';
C.InputName = 'e';
C.OutputName = 'uc';
F.InputName = 'r';
F.OutputName = 'uf';
Sum1 = sumblk('e','r','y','+-'); % e = r-y
Sum2 = sumblk('u','uf','uc','++'); % u = uf+uc
Tff = connect(P,C,F,Sum1,Sum2,'r','y');
stepplot(T,Tff)
legend('No feedforward','Feedforward')
7-55
7 Time Domain Analysis
The state-space representation keeps track of the internal delays in both models.
Gain ripple:
7-56
Analysis of Systems with Time Delays
s = tf('s');
G = exp(-5*s)/(s+1);
T = feedback(G,.5);
bodemag(T)
Gain oscillations:
G = 1 + 0.5 * exp(-3*s);
bodemag(G)
7-57
7 Time Domain Analysis
G = exp(-s) * (0.8*s^2+s+2)/(s^2+s);
T = feedback(G,1);
stepplot(T)
7-58
Analysis of Systems with Time Delays
Chaotic response:
G = 1/(s+1) + exp(-4*s);
T = feedback(1,G);
stepplot(T,150)
7-59
7 Time Domain Analysis
You can use Control System Toolbox tools to model and analyze these and other strange-
appearing artifacts of internal delays.
See Also
Related Examples
• “Closing Feedback Loops with Time Delays” on page 2-45
7-60
See Also
More About
• “Time Delays in Linear Systems” on page 2-40
• “Internal Delays” on page 2-73
7-61
8
Frequency-Domain Responses
When you perform frequency-domain analysis of a dynamic system model, you may want
one or more of the following:
• A plot of the system response as a function of frequency, or plots of pole and zero
locations.
• Numerical values of the system response in a data array.
• Numerical values of characteristics of the system response such as stability margins,
peak gains, or singular values.
Control System Toolbox frequency-domain analysis commands can obtain these results for
any kind of dynamic system model (for example, continuous or discrete, SISO or MIMO,
or arrays of models).
• pzplot, iopzplot — Plot pole and zero locations in the complex plane.
If you have a generalized state-space (genss) model of a control system, you can extract
various transfer functions from it for analysis using frequency-domain and time-domain
analysis commands. Extract responses from such models using getIOTransfer,
getLoopTransfer, getSensitivity, and getCompSensitivity.
8-2
See Also
See Also
Related Examples
• “Frequency Response of a SISO System” on page 8-4
• “Frequency Response of a MIMO System” on page 8-6
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
8-3
8 Frequency Domain Analysis
H = tf([10,21],[1,1.4,26]);
bode(H)
When you can bode without output arguments, it plots the frequency response on the
screen. Unless you specify a frequency range to plot, bode automatically chooses a
frequency range based on the system dynamics.
8-4
See Also
[mag,phase,w] = bode(H,{1,13});
When you call bode with output arguments, the command returns vectors mag and phase
containing the magnitude and phase of the frequency response. The cell array input
{1,13} tells bode to calculate the response at a grid of frequencies between 1 and 13
rad/s. bode returns the frequency points in the vector w.
See Also
bode | bodeoptions | bodeplot
Related Examples
• “Frequency Response of a MIMO System” on page 8-6
• “Numeric Values of Frequency-Domain Characteristics of SISO Model” on page 8-
13
More About
• “Frequency-Domain Responses” on page 8-2
8-5
8 Frequency Domain Analysis
Calculate the frequency response of a MIMO model and examine the size of the output.
H = rss(2,2,2);
H.InputName = 'Control';
H.OutputName = 'Temperature';
[mag,phase,w] = bode(H);
size(mag)
ans = 1×3
2 2 70
The first and second dimension of the data array mag are the number of outputs and
inputs of H. The third dimension is the number of points in the frequency vector w. (The
bode command determines this number automatically if you do not supply a frequency
vector.) Thus, mag(i,j,:) is the frequency response from the j th input of H to the i th
output, in absolute units. The phase data array phase takes the same form as mag.
bode(H)
8-6
Frequency Response of a MIMO System
bode plots the magnitude and the phase of the frequency response of each input/output
pair in H. (Because rss generates a random state-space model, you might see different
responses from those pictured.) The first column of plots shows the response from the
first input, Control(1), to each output. The second column shows the response from the
second input, Control(2), to each output.
sigma(H)
8-7
8 Frequency Domain Analysis
sigma plots the singular values of the MIMO system H as a function of frequency. The
maximum singular value at a particular frequency is the maximum gain of the system over
all linear combinations of inputs at that frequency. Singular values can provide a better
indication of the overall response, stability, and conditioning of a MIMO system than a
channel-by-channel Bode plot.
When you call sigma with output arguments, the command returns the singular values in
the data array sv. The cell array input {0.1,10} tells sigma to calculate the singular
values at a grid of frequencies between 0.1 and 10 rad/s. sigma returns these frequencies
in the vector w. Each row of sv contains the singular values of H at the frequencies of w.
8-8
See Also
See Also
bode | bodeplot | sigma | sigmaplot
Related Examples
• “Numeric Values of Frequency-Domain Characteristics of SISO Model” on page 8-
13
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
8-9
8 Frequency Domain Analysis
You can use similar procedures to display system characteristics on other types of
response plots.
Right-click anywhere in the figure and select Characteristics > Peak Response from
the menu.
8-10
Frequency-Domain Characteristics on Response Plots
A marker appears on the plot indicating the peak response. Horizontal and vertical dotted
lines indicate the frequency and magnitude of that response. The other menu options add
other system characteristics to the plot.
8-11
8 Frequency Domain Analysis
Click the marker to view the magnitude and frequency of the peak response in a datatip.
See Also
Related Examples
• “Numeric Values of Frequency-Domain Characteristics of SISO Model” on page 8-
13
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
• “Pole and Zero Locations” on page 8-16
8-12
Numeric Values of Frequency-Domain Characteristics of SISO Model
H = tf([10,21],[1,1.4,26]);
bodeplot(H)
8-13
8 Frequency Domain Analysis
[gpeak,fpeak] = getPeakGain(H);
gpeak_dB = mag2db(gpeak)
gpeak_dB = 17.7579
getPeakGain returns both the peak location fpeak and the peak gain gpeak in absolute
units. Using mag2db to convert gpeak to decibels shows that the gain peaks at almost 18
dB.
Find the band within which the system gain exceeds 0 dB, or 1 in absolute units.
wc = getGainCrossover(H,1)
wc = 2×1
1.2582
12.1843
The Bode response plot shows that the gain of H tends toward a finite value as the
frequency approaches zero. The dcgain command finds this value in absolute units.
k = dcgain(H);
Find the frequency at which the response of H rolls off to –10 dB relative to its dc value.
fb = bandwidth(H,-10);
bandwidth returns the first frequency at which the system response drops below the dc
gain by the specified value in dB.
8-14
See Also
See Also
bandwidth | getGainCrossover | getPeakGain
Related Examples
• “Pole and Zero Locations” on page 8-16
More About
• “Frequency-Domain Responses” on page 8-2
8-15
8 Frequency Domain Analysis
Examining the pole and zero locations can be useful for tasks such as stability analysis or
identifying near-canceling pole-zero pairs for model simplification. This example
compares two closed-loop systems that have the same plant and different controllers.
G = zpk([],[-5 -5 -10],100);
C1 = pid(2.9,7.1);
CL1 = feedback(G*C1,1);
C2 = pid(29,7.1);
CL2 = feedback(G*C2,1);
The controller C2 has a much higher proportional gain. Otherwise, the two closed-loop
systems CL1 and CL2 are the same.
Graphically examine the pole and zero locations of CL1 and CL2.
pzplot(CL1,CL2)
grid
8-16
Pole and Zero Locations
pzplot plots pole and zero locations on the complex plane as x and o marks,
respectively. When you provide multiple models, pzplot plots the poles and zeros of each
model in a different color. Here, there poles and zeros of CL1 are blue, and those of CL2
are green.
The plot shows that all poles of CL1 are in the left half-plane, and therefore CL1 is stable.
From the radial grid markings on the plot, you can read that the damping of the
oscillating (complex) poles is approximately 0.45. The plot also shows that CL2 contains
poles in the right half-plane and is therefore unstable.
8-17
8 Frequency Domain Analysis
zero and pole return column vectors containing the zero and pole locations of the
system.
See Also
pole | pzplot | zero
Related Examples
• “Numeric Values of Frequency-Domain Characteristics of SISO Model” on page 8-13
More About
• “Frequency-Domain Responses” on page 8-2
8-18
Assessing Gain and Phase Margins
Stability generally means that all internal signals remain bounded. This is a standard
requirement for control systems to avoid loss of control and damage to equipment. For
linear feedback systems, stability can be assessed by looking at the poles of the closed-
loop transfer function. Consider for example the SISO feedback loop:
For a unit loop gain k, you can compute the closed-loop transfer function T using:
G = tf([.5 1.3],[1 1.2 1.6 0]);
T = feedback(G,1);
ans =
-0.2305 + 1.3062i
-0.2305 - 1.3062i
-0.7389 + 0.0000i
The feedback loop for k=1 is stable since all poles have negative real parts.
8-19
8 Frequency Domain Analysis
rlocus(G)
Clicking on the point where the locus intersects the y axis reveals that this feedback loop
is stable for
8-20
Assessing Gain and Phase Margins
This range shows that with k=1, the loop gain can increase 270% before you lose stability.
Changes in the loop gain are only one aspect of robust stability. In general, imperfect
plant modeling means that both gain and phase are not known exactly. Because modeling
errors are most damaging near the gain crossover frequency (frequency where open-loop
gain is 0dB), it also matters how much phase variation can be tolerated at this frequency.
The phase margin measures how much phase variation is needed at the gain crossover
frequency to lose stability. Similarly, the gain margin measures what relative gain
variation is needed at the gain crossover frequency to lose stability. Together, these two
numbers give an estimate of the "safety margin" for closed-loop stability. The smaller the
stability margins, the more fragile stability is.
You can display the gain and phase margins on a Bode plot as follows. First create the
plot:
bode(G), grid
8-21
8 Frequency Domain Analysis
Then, right-click on the plot and select the Characteristics -> Minimum Stability
Margins submenu. Finally, click on the blue dot markers. The resulting plot is shown
below:
8-22
Assessing Gain and Phase Margins
This indicates a gain margin of about 9 dB and a phase margin of about 45 degrees. The
corresponding closed-loop step response exhibits about 20% overshoot and some
oscillations.
8-23
8 Frequency Domain Analysis
[Gm,Pm] = margin(2*G);
GmdB = 20*log10(Gm) % gain margin in dB
Pm % phase margin in degrees
GmdB =
2.7471
Pm =
8-24
Assessing Gain and Phase Margins
8.6328
and the closed-loop response has poorly damped oscillations, a sign of near instability.
Some systems have multiple gain crossover or phase crossover frequencies, which leads
to multiple gain or phase margin values. For example, consider the feedback loop
8-25
8 Frequency Domain Analysis
G = tf(20,[1 7]) * tf([1 3.2 7.2],[1 -1.2 0.8]) * tf([1 -8 400],[1 33 700]);
T = feedback(G,1);
step(T), title('Closed-loop response for k=1')
8-26
Assessing Gain and Phase Margins
To assess how robustly stable this loop is, plot its Bode response:
bode(G), grid
8-27
8 Frequency Domain Analysis
Then, right-click on the plot and select the Characteristics -> All Stability Margins
submenu to show all the crossover frequencies and associated stability margins. The
resulting plot is shown below.
8-28
Assessing Gain and Phase Margins
Note that there are two 180 deg phase crossings with corresponding gain margins of
-9.35dB and +10.6dB. Negative gain margins indicate that stability is lost by decreasing
the gain, while positive gain margins indicate that stability is lost by increasing the gain.
This is confirmed by plotting the closed-loop step response for a plus/minus 6dB gain
variation about k=1:
k1 = 2; T1 = feedback(G*k1,1);
k2 = 1/2; T2 = feedback(G*k2,1);
step(T,'b',T1,'r',T2,'g',12),
legend('k = 1','k = 2','k = 0.5')
8-29
8 Frequency Domain Analysis
The plot shows increased oscillations for both smaller and larger gain values.
You can use the command allmargin to compute all stability margins. Note that gain
margins are expressed as gain ratios, not dB. Use mag2db to convert the values to dB.
m = allmargin(G)
GainMargins_dB = mag2db(m.GainMargin)
m =
8-30
Assessing Gain and Phase Margins
GainMargins_dB =
-9.3510 10.6091
Interactive GUI
To gain additional insight into the connection between stability margins and closed-loop
responses, click on the link below to launch an interactive GUI for tuning the loop gain k
and seeing the effect on margins and closed-loop responses.
margin_gui
8-31
8 Frequency Domain Analysis
See Also
margin | pole
8-32
See Also
Related Examples
• “Pole and Zero Locations” on page 8-16
8-33
8 Frequency Domain Analysis
Many processes involve dead times, also referred to as transport delays or time lags.
Controlling such processes is challenging because delays cause linear phase shifts that
limit the control bandwidth and affect closed-loop stability.
Using the state-space representation, you can create accurate open- or closed-loop
models of control systems with delays and analyze their stability and performance without
approximation. The state-space (SS) object automatically keeps track of "internal" delays
when combining models, see the "Specifying Time Delays" tutorial for more details.
where the process model P has a 2.6 second dead time and the compensator C is a PI
controller:
To analyze the closed-loop response, construct a model T of the closed-loop transfer from
ysp to y. Because there is a delay in this feedback loop, you must convert P and C to state
space and use the state-space representation for analysis:
8-34
Analyzing Control Systems with Delays
T = feedback(P*C,1)
T =
A =
x1 x2 x3
x1 -0.36 -1.24 -0.18
x2 1 0 0
x3 0 1 0
B =
u1
x1 0.5
x2 0
x3 0
C =
x1 x2 x3
y1 0.12 0.48 0.36
D =
u1
y1 0
The result is a third-order model with an internal delay of 2.6 seconds. Internally, the
state-space object T tracks how the delay is coupled with the remaining dynamics. This
structural information is not visible to users, and the display above only gives the A,B,C,D
values when the delay is set to zero.
Use the STEP command to plot the closed-loop step response from ysp to y:
step(T)
8-35
8 Frequency Domain Analysis
The closed-loop oscillations are due to a weak gain margin as seen from the open-loop
response P*C:
margin(P*C)
8-36
Analyzing Control Systems with Delays
bode(T)
grid, title('Closed-loop frequency response')
8-37
8 Frequency Domain Analysis
To improve the design, you can try to notch out the resonance near 1 rad/s:
step(Tnotch), grid
8-38
Analyzing Control Systems with Delays
Many control design algorithms cannot handle time delays directly. A common
workaround consists of replacing delays by their Pade approximations (all-pass filters).
Because this approximation is only valid at low frequencies, it is important to compare the
true and approximate responses to choose the right approximation order and check the
approximation validity.
Use the PADE command to compute Pade approximations of LTI models with delays. For
the PI control example above, you can compare the exact closed-loop response T with the
response obtained for a first-order Pade approximation of the delay:
8-39
8 Frequency Domain Analysis
T1 = pade(T,1);
step(T,'b',T1,'r',100)
grid, legend('Exact','First-Order Pade')
The approximation error is fairly large. To get a better approximation, try a second-order
Pade approximation of the delay:
T2 = pade(T,2);
step(T,'b',T2,'r',100)
grid, legend('Exact','Second-Order Pade')
8-40
Analyzing Control Systems with Delays
The responses now match closely except for the non-minimum phase artifact introduced
by the Pade approximation.
Sensitivity Analysis
Delays are rarely known accurately, so it is often important to understand how sensitive a
control system is to the delay value. Such sensitivity analysis is easily performed using
LTI arrays and the InternalDelay property.
For example, to analyze the sensitivity of the notched PI control above, create 5 models
with delay values ranging from 2.0 to 3.0:
tau = linspace(2,3,5); % 5 delay values
Tsens = repsys(Tnotch,[1 1 5]); % 5 copies of Tnotch
8-41
8 Frequency Domain Analysis
for j=1:5
Tsens(:,:,j).InternalDelay = tau(j); % jth delay value -> jth model
end
step(Tsens)
grid, title('Closed-loop response for 5 delay values between 2.0 and 3.0')
This plot shows that uncertainty on the delay value has little effect on closed-loop
characteristics. Note that while you can change the values of internal delays, you cannot
change how many there are because this is part of the model structure. To eliminate some
internal delays, set their value to zero or use PADE with order zero:
8-42
Analyzing Control Systems with Delays
Tnotch0 = Tnotch;
Tnotch0.InternalDelay = 0;
bode(Tnotch,'b',Tnotch0,'r',{1e-2,3})
grid, legend('Delay = 2.6','No delay','Location','SouthWest')
Discretization
You can use C2D to discretize continuous-time delay systems. Available methods include
zero-order hold (ZOH), first-order hold (FOH), and Tustin. For models with internal
delays, the ZOH discretization is not always "exact," i.e., the continuous and discretized
step responses may not match:
8-43
8 Frequency Domain Analysis
Td = c2d(T,1);
step(T,'b',Td,'r')
grid, legend('Continuous','ZOH Discretization')
To correct such discretization gaps, reduce the sampling period until the continuous and
discrete responses match closely:
Td = c2d(T,0.05);
step(T,'b',Td,'r')
grid, legend('Continuous','ZOH Discretization')
8-44
Analyzing Control Systems with Delays
Note that internal delays remain internal in the discretized model and do not inflate the
model order:
order(Td)
Td.InternalDelay
ans =
8-45
8 Frequency Domain Analysis
ans =
52
The time and frequency responses of delay systems can look bizarre and suspicious to
those only familiar with delay-free LTI analysis. Time responses can behave chaotically,
Bode plots can exhibit gain oscillations, etc. These are not software quirks but real
features of such systems. Below are a few illustrations of these phenomena
Gain ripples:
G = exp(-5*s)/(s+1);
T = feedback(G,.5);
bodemag(T)
8-46
Analyzing Control Systems with Delays
Gain oscillations:
G = 1 + 0.5 * exp(-3*s);
bodemag(G)
8-47
8 Frequency Domain Analysis
G = exp(-s) * (0.8*s^2+s+2)/(s^2+s);
T = feedback(G,1);
step(T)
8-48
Analyzing Control Systems with Delays
Chaotic response:
G = 1/(s+1) + exp(-4*s);
T = feedback(1,G);
step(T)
8-49
8 Frequency Domain Analysis
See Also
margin | pade
Related Examples
• “Analyzing the Response of an RLC Circuit” on page 8-52
More About
• “Time Delays in Linear Systems” on page 2-40
8-50
See Also
8-51
8 Frequency Domain Analysis
The following figure shows the parallel form of a bandpass RLC circuit:
8-52
Analyzing the Response of an RLC Circuit
The product LC controls the bandpass frequency while RC controls how narrow the
passing band is. To build a bandpass filter tuned to the frequency 1 rad/s, set L=C=1 and
use R to tune the filter band.
The Bode plot is a convenient tool for investigating the bandpass characteristics of the
RLC network. Use tf to specify the circuit's transfer function for the values
%|R=L=C=1|:
R = 1; L = 1; C = 1;
G = tf([1/(R*C) 0],[1 1/(R*C) 1/(L*C)])
G =
s
-----------
s^2 + s + 1
bode(G), grid
8-53
8 Frequency Domain Analysis
As expected, the RLC filter has maximum gain at the frequency 1 rad/s. However, the
attenuation is only -10dB half a decade away from this frequency. To get a narrower
passing band, try increasing values of R as follows:
8-54
Analyzing the Response of an RLC Circuit
The resistor value R=20 gives a filter narrowly tuned around the target frequency of 1
rad/s.
We can confirm the attenuation properties of the circuit G2 (R=20) by simulating how this
filter transforms sine waves with frequency 0.9, 1, and 1.1 rad/s:
t = 0:0.05:250;
opt = timeoptions;
opt.Title.FontWeight = 'Bold';
subplot(311), lsim(G2,sin(t),t,opt), title('w = 1')
subplot(312), lsim(G2,sin(0.9*t),t,opt), title('w = 0.9')
subplot(313), lsim(G2,sin(1.1*t),t,opt), title('w = 1.1')
8-55
8 Frequency Domain Analysis
The waves at 0.9 and 1.1 rad/s are considerably attenuated. The wave at 1 rad/s comes
out unchanged once the transients have died off. The long transient results from the
poorly damped poles of the filters, which unfortunately are required for a narrow passing
band:
damp(pole(G2))
8-56
Analyzing the Response of an RLC Circuit
Interactive GUI
To analyze other standard circuit configurations such as low-pass and high-pass RLC
networks, click on the link below to launch an interactive GUI. In this GUI, you can
change the R,L,C parameters and see the effect on the time and frequency responses in
real time.
rlc_gui
8-57
8 Frequency Domain Analysis
See Also
bodeplot | lsim | stepplot
Related Examples
• “Joint Time-Domain and Frequency-Domain Analysis” on page 7-40
8-58
9
Sensitivity Analysis
Create an array of transfer functions representing the following low-pass filter at three
values of the roll-off frequency, a.
Create transfer function models representing the filter with roll-off frequency at a = 3, 5,
and 7.
F1 = tf(3,[1 3]);
F2 = tf(5,[1 5]);
F3 = tf(7,[1 7]);
Farray = stack(1,F1,F2,F3);
The first argument to stack specifies the array dimension along which stack builds an
array. The remaining arguments specify the models to arrange along that dimension.
Thus, Farray is a 3-by-1 array of transfer functions.
G = [F1;F2;F3];
When working with a model array that represents parameter variations, You can associate
the corresponding parameter value with each entry in the array. Set the SamplingGrid
property to a data structure that contains the name of the parameter and the sampled
parameter values corresponding with each model in the array. This assignment helps you
keep track of which model corresponds to which parameter value.
9-2
Model Array with Single Parameter Variation
Farray(:,:,1,1) [alpha=3] =
3
-----
s + 3
Farray(:,:,2,1) [alpha=5] =
5
-----
s + 5
Farray(:,:,3,1) [alpha=7] =
7
-----
s + 7
The parameter values in Farray.SamplingGrid are displayed along with the each transfer
function in the array.
Plot the frequency response of the array to examine the effect of parameter variation on
the filter behavior.
bodeplot(Farray)
9-3
9 Sensitivity Analysis
When you use analysis commands such as bodeplot on a model array, the resulting plot
shows the response of each model in the array. Therefore, you can see the range of
responses that results from the parameter variation.
See Also
stack
More About
• “Model Arrays” on page 2-100
9-4
See Also
9-5
9 Sensitivity Analysis
You can use the technique of this example to create higher-dimensional arrays with
variations of more parameters. Such arrays are useful for studying the effects of multiple-
parameter variations on system response.
depends on two parameters: the damping ratio, , and the natural frequency, . If both
and vary, you obtain multiple transfer functions of the form:
Preallocate memory for the model array. Preallocating memory is an optional step that
can enhance computation efficiency. To preallocate, create a model array of the required
size and initialize its entries to zero.
H = tf(zeros(1,1,3,3));
In this example, there are three values for each parameter in the transfer function H.
Therefore, this command creates a 3-by-3 array of single-input, single-output (SISO) zero
transfer functions.
9-6
Model Array with Variations in Two Parameters
H is a 3-by-3 array of transfer functions. varies as you move from model to model along
a single column of H. The parameter varies as you move along a single row.
Plot the step response of H to see how the parameter variation affects the step response.
stepplot(H)
9-7
9 Sensitivity Analysis
You can set the SamplingGrid property of the model array to help keep track of which
set of parameter values corresponds to which entry in the array. To do so, create a grid of
parameter values that matches the dimensions of the array. Then, assign these values to
H.SamplingGrid with the parameter names.
[zetagrid,wgrid] = ndgrid(zeta,w);
H.SamplingGrid = struct('zeta',zetagrid,'w',wgrid);
When you display H, the parameter values in H.SamplingGrid are displayed along with
the each transfer function in the array.
See Also
ndgrid
More About
• “Model Arrays” on page 2-100
• “Study Parameter Variation by Sampling Tunable Model” on page 9-9
9-8
Study Parameter Variation by Sampling Tunable Model
Sample this filter at varying values of the damping constant and the natural frequency
. Create a parametric model of the filter by using tunable elements for and .
wn = realp('wn',3);
zeta = realp('zeta',0.8);
F = tf(wn^2,[1 2*zeta*wn wn^2])
F =
Type "ss(F)" to see the current value, "get(F)" to see all properties, and "F.Blocks" t
F is a genss model with two tunable Control Design Blocks, the realp blocks wn and
zeta. The blocks wn and zeta have initial values of 3 and 0.8, respectively.
Here, sampleBlock samples the model independently over the two values and three
values. Thus, Fsample is a 2-by-3 array of state-space models. Each entry in the array
is a state-space model that represents F evaluated at the corresponding (wn, zeta) pair.
For example, Fsample(:,:,2,3) has wn = 5 and zeta = 1.0.
Set the SamplingGrid property of the model array to help keep track of which set of
parameter values corresponds to which entry in the array. To do so, create a grid of
9-9
9 Sensitivity Analysis
parameter values that matches the dimensions of the array. Then, assign these values to
Fsample.SamplingGrid in a structure with the parameter names.
[wngrid,zetagrid] = ndgrid(wnvals,zetavals);
Fsample.SamplingGrid = struct('wn',wngrid,'zeta',zetagrid);
The ndgrid command produces the full 2-by-3 grid of (wn, zeta) combinations. When
you display Fsample in the command window, the parameter values in
Fsample.SamplingGrid are displayed along with the each transfer function in the
array. The parameter information is also available in response plots. For instance,
examine the step response of Fsample.
stepplot(Fsample)
9-10
See Also
The step response plots show the variation in the natural frequency and damping
constant across the six models in the array. When you click on one of the responses in the
plot, the datatip includes the corresponding wn and zeta values as specified in
Fsample.SamplingGrid.
See Also
sampleBlock
More About
• “Models with Tunable Coefficients” on page 1-19
9-11
9 Sensitivity Analysis
Time delays are rarely known accurately, so it is often important to understand how
sensitive a control system is to the delay value. Such sensitivity analysis is easily
performed using LTI arrays and the InternalDelay property. For example, consider the
notched PI control system developed in "PI Control Loop with Dead Time" from the
example "Analyzing Control Systems with Delays." The following commands create an LTI
model of that closed-loop system, a third-order plant with an input delay, a PI controller
and a notch filter.
s = tf('s');
G = exp(-2.6*s)*(s+3)/(s^2+0.3*s+1);
C = 0.06 * (1 + 1/s);
T = feedback(ss(G*C),1);
notch = tf([1 0.2 1],[1 .8 1]);
C = 0.05 * (1 + 1/s);
Tnotch = feedback(ss(G*C*notch),1);
Tnotch.InternalDelay
ans = 2.6000
The 2.6-second input delay of the plant G becomes an internal delay of 2.6 s in the closed-
loop system. To examine the sensitivity of the responses of Tnotch to variations in this
delay, create an array of copies of Tnotch. Then, vary the internal delay across the array.
The array Tsens contains five models with internal delays that range from 2.0 to 3.0.
stepplot(Tsens)
9-12
See Also
The plot shows that uncertainty on the delay value has a small effect on closed-loop
characteristics.
See Also
More About
• “Time Delays in Linear Systems” on page 2-40
9-13
10
For example, a PID controller is passive because the control signal (the output) moves in
the same direction as the error signal (the input). But a PID controller with delay is not
passive, because the control signal can move in the opposite direction from the error, a
potential cause of instability.
Most physical systems are passive. The Passivity Theorem holds that the negative-
feedback interconnection of two strictly passive systems is passive and stable. As a result,
it can be desirable to enforce passivity of the controller for a passive system, or to
passivate the operator of a passive system, such as the driver of a car.
In practice, passivity can easily be destroyed by the phase lags introduced by sensors,
actuators, and communication delays. These problems have led to extension of the
Passivity Theorem that consider excesses or shortages of passivity, frequency-dependent
measures of passivity, and a mix of passivity and small-gain properties.
Passive Systems
where denotes the transpose of . For physical systems, the integral typically
represents the energy going into the system,. Thus passive systems are systems that only
consume or dissipate energy. As a result, passive systems are intrinsically stable.
10-2
About Passivity and Passivity Indices
For SISO systems, this is saying that at all frequencies, so the entire
Nyquist plot lies in the right-half plane.
Passive systems have the following important properties for control purposes:
10-3
10 Passivity and Conic Sectors
For stability, knowing whether a system is passive or not does not tell the full story. It is
often desirable to know by how much it is passive or fails to be passive. In addition, a
shortage of passivity in the plant can be compensated by an excess of passivity in the
controller, and vice versa. It is therefore important to measure the excess or shortage of
passivity, and this is where passivity indices come into play.
There are different types of indices with different applications. One class of indices
measure the excess or shortage of passivity in a particular direction of the input/output
space. For example, the input passivity index is defined as the largest such that:
for all trajectories and . The system G is input strictly passive (ISP)
when , and has a shortage of passivity when . The input passivity index is also
called the input feedforward passivity (IFP) index because it corresponds to the minimum
static feedforward action needed to make the system passive.
10-4
About Passivity and Passivity Indices
where denotes the smallest eigenvalue. In the SISO case, is the abscissa of the
leftmost point on the Nyquist curve.
Similarly, the output passivity index is defined as the largest such that:
for all trajectories and . The system G is output strictly passive (OSP)
when , and has a shortage of passivity when . The output passivity index is
also called the output feedback passivity (OFP) index because it corresponds to the
minimum static feedback action needed to make the system passive.
In the SISO case, is the abscissa of the leftmost point on the Nyquist curve of .
Combining these two notions leads to the I/O passivity index, which is the largest such
that:
10-5
10 Passivity and Conic Sectors
A system with is very strictly passive. More generally, we can define the index in
the direction as the largest such that:
The input, output, and I/O passivity indices all correspond to special choices of and
are collectively referred to as directional passivity indices. You can use
getPassiveIndex to compute any of these indices for linear systems in either
parametric or FRD form. You can also use passiveplot to plot the input, output, or I/O
passivity indices as a function of frequency. This plot provides insight into which
frequency bands have weaker or stronger passivity.
There are many results quantifying how the input and output passivity indices propagate
through parallel, series, or feedback interconnections. There are also results quantifying
the excess of input or output passivity needed to compensate a given shortage of passivity
in a feedback loop. For details, see:
10-6
About Passivity and Passivity Indices
for all trajectories and . When is minimum phase, you can use
passiveplot to plot the principal gains of . This plot is
entirely analogous to the singular value plot (see sigma), and shows how the degree of
passivity changes with frequency and direction.
The following result is analogous to the Small Gain Theorem for feedback loops. It gives a
simple condition on R-indices for compensating a shortage of passivity in one system by
an excess of passivity in the other.
Small-R Theorem: Let and be two linear systems with passivity R-indices
and , respectively. If , then the negative feedback interconnection of
and is stable.
10-7
10 Passivity and Conic Sectors
See Also
getPassiveIndex | isPassive | passiveplot
Related Examples
• “Passivity Indices” on page 10-19
• “Parallel Interconnection of Passive Systems” on page 10-24
• “Series Interconnection of Passive Systems” on page 10-27
• “Feedback Interconnection of Passive Systems” on page 10-31
• “About Sector Bounds and Sector Indices” on page 10-9
10-8
About Sector Bounds and Sector Indices
In its simplest form, a conic sector is the 2-D region delimited by two lines, and
.
where is a 2x2 symmetric indefinite matrix ( has one positive and one negative
eigenvalue). We call the sector matrix. This concept generalizes to higher dimensions.
In an N-dimensional space, a conic sector is a set:
10-9
10 Passivity and Conic Sectors
Sector Bounds
Sector bounds are constraints on the behavior of a system. Gain constraints and passivity
constraints are special cases of sector bounds. If for all nonzero input trajectories ,
the output trajectory of a linear system satisfies:
then the output trajectories of lie in the conic sector with matrix . Selecting different
matrices imposes different conditions on the system's response. For example, consider
trajectories and the following values:
In other words, passivity is a particular sector bound on the system defined by:
10-10
About Sector Bounds and Sector Indices
Frequency-Domain Condition
Because the time-domain condition must hold for all , deriving an equivalent
frequency-domain bound takes a little care and is not always possible. Let the following:
be (any) decomposition of the indefinite matrix into its positive and negative parts.
When is square and minimum phase (has no unstable zeros), the time-domain
condition:
It is therefore enough to check the sector inequality for real frequencies. Using the
decomposition of , this is also equivalent to:
Note that is square when has as many negative eigenvalues as input channels in
. If this condition is not met, it is no longer enough (in general) to just look at real
frequencies. Note also that if is square, then it must be minimum phase for the
sector bound to hold.
10-11
10 Passivity and Conic Sectors
For instance, examine the sector plot of a 2-output, 2-input system for a particular sector.
rng(4);
H = rss(3,4,2);
Q = [-5.12 2.16 -2.04 2.17
2.16 -1.22 -0.28 -1.11
-2.04 -0.28 -3.35 0.00
2.17 -1.11 0.00 0.18];
sectorplot(H,Q)
10-12
About Sector Bounds and Sector Indices
We can extend the notion of relative passivity index to arbitrary sectors. Let be an
LTI system, and let:
To understand the geometrical interpretation of the R-index, consider the family of cones
10-13
10 Passivity and Conic Sectors
In the diagram,
and
When is square and minimum phase, the R-index can also be characterized in
the frequency domain as the smallest such that:
10-14
About Sector Bounds and Sector Indices
In other words, the R-index is the peak gain of the (stable) transfer function
Similarly, we can extend the notion of directional passivity index to arbitrary sectors.
Given a conic sector with matrix , and a direction , the directional sector index is
the largest such that for all output trajectories :
The directional sector index measures by how much we need to deform the sector in the
direction to make it fit tightly around the output trajectories of . The sector bound
is satisfied if and only if the directional index is positive.
Common Sectors
There are many ways to specify sector bounds. Next we review commonly encountered
expressions and give the corresponding system and sector matrix for the standard
form used by getSectorIndex and sectorplot:
10-15
10 Passivity and Conic Sectors
Passivity
Gain constraint
Ratio of distances
The underlying conic sector is symmetric with respect to . Similarly, the "exterior"
constraint,
10-16
About Sector Bounds and Sector Indices
Double inequality
When dealing with static nonlinearities, it is common to consider conic sectors of the form
where is the nonlinearity output. While this relationship is not a sector bound
per se, it clearly implies:
along all I/O trajectories and for all . This condition in turn is equivalent to a sector
bound with:
Product form
correspond to:
10-17
10 Passivity and Conic Sectors
QSR dissipative
See Also
getSectorCrossover | getSectorIndex | sectorplot
Related Examples
• “About Passivity and Passivity Indices” on page 10-2
10-18
Passivity Indices
Passivity Indices
This example shows how to compute various measures of passivity for linear time-
invariant systems.
Passive Systems
•
The input passivity index is defined as the largest such that
The system G is "input strictly passive" (ISP) when . is also called the "input
feedforward passivity" (IFP) index and corresponds to the minimum feedforward action
needed to make the system passive.
•
The output passivity index is defined as the largest such that
10-19
10 Passivity and Conic Sectors
The system G is "output strictly passive" (OSP) when . is also called the "output
feedback passivity" (OFP) index and corresponds to the minimum feedback action needed
to make the system passive.
•
The I/O passivity index is defined as the largest such that
Circuit Example
Consider the following example. We take the current as the input and the voltage as
the output. Based on Kirchhoff's current and voltage law, we obtain the transfer function
for ,
Let , and .
R = 2; L = 1; C = 0.1;
s = tf('s');
G = (L*s+R)*(R*s+1/C)/(L*s^2 + 2*R*s+1/C);
10-20
Passivity Indices
PF = isPassive(G)
PF = logical
1
nu = 2
rho = 0.2857
tau = 0.2642
Frequency-Domain Characterization
The smallest eigenvalue of the left-hand-side is related to the input passivity index :
10-21
10 Passivity and Conic Sectors
Verify this for the circuit example. Plot the Nyquist plot of the circuit transfer function.
nyquist(G)
The entire Nyquist plot lies in the right-half plane so is positive real. The leftmost
point on the Nyquist curve is so the input passivity index is , the
same value we obtained earlier. Similarly, the leftmost point on the Nyquist curve for
gives the output passivity index value .
10-22
See Also
The relative passivity index (R-index) is the peak gain over frequency of
when is minimum phase, and otherwise:
The system is passive if and only if , and the smaller is, the more passive
the system is. Use getPassiveIndex to compute the R-index for the circuit example.
R = getPassiveIndex(G)
R = 0.5556
The resulting value indicates that the circuit is a very passive system.
See Also
getPassiveIndex | isPassive
Related Examples
• “About Passivity and Passivity Indices” on page 10-2
• “Parallel Interconnection of Passive Systems” on page 10-24
• “Series Interconnection of Passive Systems” on page 10-27
• “Feedback Interconnection of Passive Systems” on page 10-31
10-23
10 Passivity and Conic Sectors
If both systems and are passive, then the interconnected system is guaranteed
to be passive. Take for example
ans = logical
1
G2 = tf([1,2,1],[1,3,10]);
isPassive(G2)
ans = logical
1
10-24
Parallel Interconnection of Passive Systems
ans = logical
1
There is a relationship between the passivity indices of and and the passivity
indices of the interconnected system . Let and denote the input passivity indices
for and , and let and denote the output passivity indices. If all these indices
are nonnegative, then the input passivity index and the output passivity index for the
parallel interconnection satisfy
In other words, we can infer some minimum level of input and output passivity for the
parallel connection from the input and output passivity indices of and . For
details, see the paper by Yu, H., "Passivity and dissipativity as design and analysis tools
for networked control systems," Chapter 2, PhD Thesis, University of Notre Dame, 2012.
Verify the lower bound for the input passivity index .
% Input passivity index for G1
nu1 = getPassiveIndex(G1,'input');
% Input passivity index for G2
nu2 = getPassiveIndex(G2,'input');
% Input passivity index for H
nu = getPassiveIndex(H,'input')
nu = 0.3777
% Lower bound
nu1+nu2
ans = 0.1474
10-25
10 Passivity and Conic Sectors
Similarly, verify the lower bound for the output passivity index of .
rho = 0.6450
% Lower bound
rho1*rho2/(rho1+rho2)
ans = 0.2098
See Also
getPassiveIndex | isPassive
Related Examples
• “About Passivity and Passivity Indices” on page 10-2
• “Series Interconnection of Passive Systems” on page 10-27
• “Feedback Interconnection of Passive Systems” on page 10-31
10-26
Series Interconnection of Passive Systems
ans = logical
1
G2 = tf([1,1,5,.1],[1,2,3,4]);
isPassive(G2)
ans = logical
1
10-27
10 Passivity and Conic Sectors
H = G2*G1;
isPassive(H)
ans = logical
0
This is confirmed by verifying that the Nyquist plot of is not positive real.
nyquist(H)
10-28
Series Interconnection of Passive Systems
While the series interconnection of passive systems is not passive in general, there is a
relationship between the passivity indices of and and the passivity indices of
. Let and denote the input passivity indices for and , and let
and denote the output passivity indices. If all these indices are positive, then the input
passivity index and the output passivity index for the series interconnection
satisfy
In other words, the shortage of passivity at the inputs or outputs of is no worse than
the right-hand-side expressions. For details, see the paper by Arcak, M. and Sontag, E.D.,
"Diagonal stability of a class of cyclic systems and its connection with the secant
criterion," Automatica, Vol 42, No. 9, 2006, pp. 1531-1537. Verify these lower bounds for
the example above.
% Output passivity index for G1
rho1 = getPassiveIndex(G1,'output');
% Output passivity index for G2
rho2 = getPassiveIndex(G2,'output');
% Input passivity index for H=G2*G1
nu = getPassiveIndex(H,'input')
nu = -1.2875
% Lower bound
-0.125/(rho1*rho2)
ans = -2.4119
Similarly, verify the lower bound for the output passivity index of .
% Input passivity index for G1
nu1 = getPassiveIndex(G1,'input');
% Input passivity index for G2
nu2 = getPassiveIndex(G2,'input');
% Output passivity index for H=G2*G1
rho = getPassiveIndex(H,'output')
10-29
10 Passivity and Conic Sectors
rho = -0.6966
% Lower bound
-0.125/(nu1*nu2)
ans = -5.9420
See Also
getPassiveIndex | isPassive
Related Examples
• “About Passivity and Passivity Indices” on page 10-2
• “Parallel Interconnection of Passive Systems” on page 10-24
• “Feedback Interconnection of Passive Systems” on page 10-31
10-30
Feedback Interconnection of Passive Systems
If both systems and are passive, then the interconnected system is guaranteed
to be passive. Take for example
G1 = tf([1,1,1],[1,1,4]);
isPassive(G1)
ans = logical
1
G2 = tf([1,2],[1,5]);
isPassive(G2)
ans = logical
1
10-31
10 Passivity and Conic Sectors
H = feedback(G1,G2);
isPassive(H)
ans = logical
1
nyquist(H)
10-32
Feedback Interconnection of Passive Systems
There is a relationship between the passivity indices of and and the passivity
indices of the interconnected system . Let and denote the input passivity indices
for and , and let and denote the output passivity indices. If all these indices
are positive, then the input passivity index and the output passivity index for the
feedback interconnection satisfy
In other words, we can infer some minimum level of input and output passivity for the
closed-loop system from the input and output passivity indices of and . For
details, see the paper by Zhu, F. and Xia, M and Antsaklis, P.J., "Passivity analysis and
passivation of feedback systems using passivity indices," American Control Conference ,
2014, pp. 1833-1838. Verify the lower bound for the input passivity index .
% Input passivity index for G1
nu1 = getPassiveIndex(G1,'input');
% Output passivity index for G2
rho2 = getPassiveIndex(G2,'output');
% Input passivity index for H
nu = getPassiveIndex(H,'input')
nu = 0.1293
% Lower bound
nu1*rho2/(nu1+rho2)
ans = 2.4923e-06
Similarly, verify the lower bound for the output passivity index of .
% Output passivity index for G1
rho1 = getPassiveIndex(G1,'output');
% Input passivity index for G2
nu2 = getPassiveIndex(G2,'input');
% Output passivity index for H
rho = getPassiveIndex(H,'output')
10-33
10 Passivity and Conic Sectors
rho = 0.4485
% Lower bound
rho1+nu2
ans = 0.4000
See Also
getPassiveIndex | isPassive
Related Examples
• “About Passivity and Passivity Indices” on page 10-2
• “Parallel Interconnection of Passive Systems” on page 10-24
• “Series Interconnection of Passive Systems” on page 10-27
• “Passive Control with Communication Delays” on page 17-282
10-34
Control Design
35
11
As a first pass, create a model of the plant and design a simple PI controller for it.
C_pi =
1
Kp + Ki * ---
s
info =
Stable: 1
CrossoverFrequency: 0.5205
PhaseMargin: 60.0000
C_pi is a pid controller object that represents a PI controller. The fields of info show
that the tuning algorithm chooses an open-loop crossover frequency of about 0.52 rad/s.
Examine the closed-loop step response (reference tracking) of the controlled system.
11-2
PID Controller Design at the Command Line
To improve the response time, you can set a higher target crossover frequency than the
result that pidtune automatically selects, 0.52. Increase the crossover frequency to 1.0.
[C_pi_fast,info] = pidtune(sys,'PI',1.0)
C_pi_fast =
1
Kp + Ki * ---
s
11-3
11 PID Controller Design
info =
Stable: 1
CrossoverFrequency: 1
PhaseMargin: 43.9973
The new controller achieves the higher crossover frequency, but at the cost of a reduced
phase margin.
T_pi_fast = feedback(C_pi_fast*sys,1);
step(T_pi,T_pi_fast)
axis([0 30 0 1.4])
legend('PI','PI,fast')
11-4
PID Controller Design at the Command Line
This reduction in performance results because the PI controller does not have enough
degrees of freedom to achieve a good phase margin at a crossover frequency of 1.0 rad/s.
Adding a derivative action improves the response.
Design a PIDF controller for Gc with the target crossover frequency of 1.0 rad/s.
[C_pidf_fast,info] = pidtune(sys,'PIDF',1.0)
C_pidf_fast =
1 s
Kp + Ki * --- + Kd * --------
s Tf*s+1
11-5
11 PID Controller Design
info =
Stable: 1
CrossoverFrequency: 1
PhaseMargin: 60.0000
The fields of info show that the derivative action in the controller allows the tuning
algorithm to design a more aggressive controller that achieves the target crossover
frequency with a good phase margin.
Compare the closed-loop step response and disturbance rejection for the fast PI and PIDF
controllers.
T_pidf_fast = feedback(C_pidf_fast*sys,1);
step(T_pi_fast, T_pidf_fast);
axis([0 30 0 1.4]);
legend('PI,fast','PIDF,fast');
11-6
PID Controller Design at the Command Line
You can compare the input (load) disturbance rejection of the controlled system with the
fast PI and PIDF controllers. To do so, plot the response of the closed-loop transfer
function from the plant input to the plant output.
S_pi_fast = feedback(sys,C_pi_fast);
S_pidf_fast = feedback(sys,C_pidf_fast);
step(S_pi_fast,S_pidf_fast);
axis([0 50 0 0.4]);
legend('PI,fast','PIDF,fast');
11-7
11 PID Controller Design
This plot shows that the PIDF controller also provides faster disturbance rejection.
See Also
pid | pidtune
More About
• “Choosing a PID Controller Design Tool”
• “Designing Cascade Control System with PI Controllers” on page 11-9
• “PID Controller Design for Fast Reference Tracking”
11-8
Designing Cascade Control System with PI Controllers
Controller C1 in the outer loop is the primary controller that regulates the primary
controlled variable y1 by setting the set-point of the inner loop. Controller C2 in the inner
loop is the secondary controller that rejects disturbance d2 locally before it propagates to
P1. For a cascade control system to function properly, the inner loop must respond much
faster than the outer loop.
In this example, you will design a single loop control system with a PI controller and a
cascade control system with two PI controllers. The responses of the two control systems
are compared for both reference tracking and disturbance rejection.
Plant
11-9
11 PID Controller Design
P2 = zpk([],-2,3);
P1 = zpk([],[-1 -1 -1],10);
Use pidtune command to design a PI controller in standard form for the whole plant
model P = P1 * P2.
The desired open loop bandwidth is 0.2 rad/s, which roughly corresponds to the response
time of 10 seconds.
C =
1 1
Kp * (1 + ---- * ---)
Ti s
11-10
Designing Cascade Control System with PI Controllers
The best practice is to design the inner loop controller C2 first and then design the outer
loop controller C1 with the inner loop closed. In this example, the inner loop bandwidth is
selected as 2 rad/s, which is ten times higher than the desired outer loop bandwidth. In
order to have an effective cascade control system, it is essential that the inner loop
responds much faster than the outer loop.
C2 = pidtune(P2,pidstd(1,1),2);
C2
C2 =
1 1
Kp * (1 + ---- * ---)
Ti s
Tune outer-loop controller C1 with the same bandwidth as the single loop system.
C1 =
1 1
Kp * (1 + ---- * ---)
Ti s
Performance Comparison
First, plot the step reference tracking responses for both control systems.
11-11
11 PID Controller Design
Secondly, plot the step disturbance rejection responses of d2 for both control systems.
11-12
Designing Cascade Control System with PI Controllers
11-13
11 PID Controller Design
From the two response plots you can conclude that the cascade control system performs
much better in rejecting disturbance d2 while the set-point tracking performances are
almost identical.
See Also
pidstd | pidtune
More About
• “Choosing a PID Controller Design Tool”
• “PID Controller Design at the Command Line” on page 11-2
• “Tune PID Controller to Favor Reference Tracking or Disturbance Rejection
(Command Line)”
11-14
Tune 2-DOF PID Controller (Command Line)
2-DOF PID controllers include setpoint weighting on the proportional and derivative
terms. Compared to a 1-DOF PID controller, a 2-DOF PID controller can achieve better
disturbance rejection without significant increase of overshoot in setpoint tracking. A
typical control architecture using a 2-DOF PID controller is shown in the following
diagram.
For this example, design a 2-DOF controller for the plant given by:
Suppose that your target bandwidth for the system is 1.5 rad/s.
wc = 1.5;
G = tf(1,[1 0.5 0.1]);
C2 = pidtune(G,'PID2',wc)
C2 =
1
u = Kp (b*r-y) + Ki --- (r-y) + Kd*s (c*r-y)
s
11-15
11 PID Controller Design
Using the type 'PID2' causes pidtune to generate a 2-DOF controller, represented as a
pid2 object. The display confirms this result. pidtune tunes all controller coefficients,
including the setpoint weights b and c, to balance performance and robustness.
To compute the closed-loop response, note that a 2-DOF PID controller is a 2-input, 1-
output dynamic system. You can resolve the controller into two channels, one for the
reference signal and one for the feedback signal, as shown in the diagram. (See
“Continuous-Time 2-DOF PID Controller Representations” on page 2-16 for more
information.)
Decompose the controller into the components Cr and Cy, and use them to compute the
closed-loop response from r to y.
C2tf = tf(C2);
Cr = C2tf(1);
Cy = C2tf(2);
T2 = Cr*feedback(G,Cy,+1);
S2 = feedback(G,Cy,+1);
11-16
Tune 2-DOF PID Controller (Command Line)
For comparison, design a 1-DOF PID controller with the same bandwidth and compute the
corresponding transfer functions. Then compare the step responses.
C1 = pidtune(G,'PID',wc);
T1 = feedback(G*C1,1);
S1 = feedback(G,C1);
subplot(2,1,1)
stepplot(T1,T2)
title('Reference Tracking')
subplot(2,1,2)
stepplot(S1,S2)
title('Disturbance Rejection')
legend('1-DOF','2-DOF')
11-17
11 PID Controller Design
The plots show that adding the second degree of freedom eliminates the overshoot in the
reference-tracking response without any cost to disturbance rejection. You can improve
disturbance rejection too using the DesignFocus option. This option causes pidtune to
favor disturbance rejection over setpoint tracking.
opt = pidtuneOptions('DesignFocus','disturbance-rejection');
C2dr = pidtune(G,'PID2',wc,opt)
C2dr =
1
u = Kp (b*r-y) + Ki --- (r-y) + Kd*s (c*r-y)
s
With the default balanced design focus, pidtune selects a b value between 0 and 1. For
this plant, when you change design focus to favor disturbance rejection, pidtune sets b
= 0 and c = 0. Thus, pidtune automatically generates an I-PD controller to optimize for
disturbance rejection. (Explicitly specifying an I-PD controller without setting the design
focus yields a similar controller.)
C2dr_tf = tf(C2dr);
Cdr_r = C2dr_tf(1);
Cdr_y = C2dr_tf(2);
T2dr = Cdr_r*feedback(G,Cdr_y,+1);
S2dr = feedback(G,Cdr_y,+1);
subplot(2,1,1)
stepplot(T1,T2,T2dr)
title('Reference Tracking')
subplot(2,1,2)
stepplot(S1,S2,S2dr);
title('Disturbance Rejection')
legend('1-DOF','2-DOF','2-DOF rejection focus')
11-18
Tune 2-DOF PID Controller (Command Line)
The plots show that the disturbance rejection is further improved compared to the
balanced 2-DOF controller. This improvement comes with some sacrifice of reference-
tracking performance, which is slightly slower. However, the reference-tracking response
still has no overshoot.
Thus, using 2-DOF control can improve disturbance rejection without sacrificing as much
reference tracking performance as 1-DOF control. These effects on system performance
depend strongly on the properties of your plant. For some plants and some control
11-19
11 PID Controller Design
bandwidths, using 2-DOF control or changing the design focus has less or no impact on
the tuned result.
See Also
pid2 | pidtune
More About
• “Designing PID Controllers with PID Tuner”
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
• “Tune 2-DOF PID Controller (PID Tuner)” on page 11-21
• “Analyze Design in PID Tuner”
11-20
Tune 2-DOF PID Controller (PID Tuner)
In this example, you represent the plant as an LTI model on page 1-13. For information
about using PID Tuner to tune a PID Controller (2DOF) block in a Simulink model, see
“Design Two-Degree-of-Freedom PID Controllers” (Simulink Control Design).
2-DOF PID controllers include setpoint weighting on the proportional and derivative
terms. Compared to a 1-DOF PID controller, a 2-DOF PID controller can achieve better
disturbance rejection without significant increase of overshoot in setpoint tracking. A
typical control architecture using a 2-DOF PID controller is shown in the following
diagram.
For this example, first design a 1-DOF controller for the plant given by:
1
G ( s) = .
2
s + 0 .5 s + 0 .1
11-21
11 PID Controller Design
Suppose for this example that your application requires a faster response than the PID
Tuner initial design. In the text box next to the Response Time slider, enter 2.
11-22
Tune 2-DOF PID Controller (PID Tuner)
The resulting response is fast, but has a considerable amount of overshoot. Design a 2-
DOF controller to improve the overshoot. First, set the 1-DOF controller as the baseline
controller for comparison. Click the Export arrow and select Save as Baseline.
11-23
11 PID Controller Design
11-24
Tune 2-DOF PID Controller (PID Tuner)
PID Tuner generates a 2-DOF controller with the same target response time. The
controller parameters displayed at the bottom right show that PID Tuner tunes all
11-25
11 PID Controller Design
Adding the second degree of freedom eliminates the overshoot in the reference tracking
response. Next, add a step response plot to compare the disturbance rejection
performance of the two controllers. Select Add Plot > Input Disturbance Rejection.
11-26
Tune 2-DOF PID Controller (PID Tuner)
PID Tuner tiles the disturbance-rejection plot side by side with the reference-tracking
plot.
11-27
11 PID Controller Design
You can improve disturbance rejection too by changing the PID Tuner design focus. First,
click the Export arrow and select Save as Baseline again to set the 2-DOF
controller as the baseline for comparison.
Change the PID Tuner design focus to favor reference tracking without changing the
11-28
Tune 2-DOF PID Controller (PID Tuner)
PID Tuner automatically retunes the controller coefficients with a focus on disturbance-
rejection performance.
11-29
11 PID Controller Design
With the default balanced design focus, PID Tuner selects a b value between 0 and 1. For
this plant, when you change design focus to favor disturbance rejection, PID Tuner sets
b = 0 and c = 0. Thus, PID Tuner automatically generates an I-PD controller to optimize
for disturbance rejection. (Explicitly specifying an I-PD controller without setting the
design focus yields a similar controller.)
The response plots show that with the change in design focus, the disturbance rejection is
further improved compared to the balanced 2-DOF controller. This improvement comes
with some sacrifice of reference-tracking performance, which is slightly slower. However,
the reference-tracking response still has no overshoot.
Thus, using 2-DOF control can improve disturbance rejection without sacrificing as much
reference tracking performance as 1-DOF control. These effects on system performance
depend strongly on the properties of your plant and the speed of your controller. For some
plants and some control bandwidths, using 2-DOF control or changing the design focus
has less or no impact on the tuned result.
See Also
pidTuner
More About
• “Designing PID Controllers with PID Tuner”
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
• “Tune 2-DOF PID Controller (Command Line)” on page 11-15
• “Analyze Design in PID Tuner”
11-30
PID Controller Types for Tuning
• For command-line tuning, provide the type argument to the pidtune command. For
example, C = pidtune(G,'PI') tunes a PI controller for plant G.
• For tuning in PID Tuner:
• Provide the type argument to the pidTuner command when you open PID Tuner.
For example, pidTuner(G,'PIDF2') opens PID Tuner with an initial design that
is a 2-DOF PID controller with a filter on the derivative term.
• Provide the baseline-controller Cbase argument to the pidTuner command when
you open PID Tuner. PID Tuner designs a controller of the same type as Cbase.
For example, suppose C0 is a pid controller object that has proportional and
derivative action only (PD controller). Then, pidTuner(G,C0) opens PID Tuner
with an initial design that is a PD controller.
• In PID Tuner, use the Type menu to change controller types.
11-31
11 PID Controller Design
11-32
PID Controller Types for Tuning
1-DOF Controllers
The following table summarizes the available 1-DOF PID controller types and provides
representative controller formulas for parallel form. The standard-form and discrete-time
formulas are analogous.
2-DOF Controllers
PID Tuner can automatically design 2-DOF PID controller types with free setpoint
weights. The following table summarizes the 2-DOF controller types in PID Tuner. The
standard-form and discrete-time formulas are analogous. For more information about 2-
11-33
11 PID Controller Design
Ki
u = K p ( br - y ) + ( r - y) + K d s ( cr - y)
s
11-34
See Also
If you set b = 0 and c = 0, then changes in the setpoint r do not feed through directly to
either the proportional or the derivative terms in u. The b = 0, c = 0 controller is called
an I-PD type controller. I-PD controllers are also useful for improving disturbance
rejection.
Use PID Tuner to design the fixed-setpoint-weight controller types summarized in the
following table. The standard-form and discrete-time formulas are analogous.
See Also
pidTuner | pidtune
11-35
11 PID Controller Design
More About
• “Designing PID Controllers with PID Tuner”
• “Proportional-Integral-Derivative (PID) Controllers” on page 2-13
• “Two-Degree-of-Freedom PID Controllers” on page 2-16
• “PID Controller Design at the Command Line” on page 11-2
• “PID Controller Design for Fast Reference Tracking”
• “Tune 2-DOF PID Controller (Command Line)” on page 11-15
• “Tune 2-DOF PID Controller (PID Tuner)” on page 11-21
11-36
12
12-2
See Also
See Also
More About
• “PID Controller Tuning”
• “Classical Control Design”
• “Tuning with Control System Tuner”
• “Programmatic Tuning”
12-3
12 Classical Control Design
When using graphical tuning, you can modify the compensator either directly from the
editor plots or using the compensator editor. A common design approach is to roughly
tune your compensator using the editor plots, and then use the compensator editor to
fine-tune the compensator parameters. For more information, see “Edit Compensator
Dynamics” on page 12-93
The graphical tuning methods are not mutually exclusive. For example, you can tune your
compensator using both the Bode editor and root locus editor simultaneously. This option
is useful when designing to both time-domain and frequency-domain specifications.
12-4
Control System Designer Tuning Methods
12-5
12 Classical Control Design
A common design approach is to generate an initial compensator using PID tuning, LQG
synthesis, loop shaping, or IMC tuning. You can then improve the compensator
performance using either optimization-based tuning or graphical tuning.
For more information on automated tuning methods, see “Design Compensator Using
Automated Tuning Methods” on page 12-100.
Knowing the properties of the effective plant seen by your compensator can help you
understand which tuning methods work for your system. For example, some automated
Ÿ
tuning methods apply only to compensators whose open loops ( L = C P ) have stable
12-6
Control System Designer Tuning Methods
Ÿ
effective plants ( P ). Also, for tuning methods such as IMC and loop shaping, the
maximum controller order depends on the dynamics of the effective plant.
12-7
12 Classical Control Design
See Also
Control System Designer
Related Examples
• “Bode Diagram Design” on page 12-49
• “Root Locus Design” on page 12-64
• “Nichols Plot Design” on page 12-80
• “Design Compensator Using Automated Tuning Methods” on page 12-100
12-8
Design Requirements
Design Requirements
This topic describes time-domain and frequency-domain design requirements available in
Control System Designer. Each requirement defines an exclusion region, indicated by a
yellow shaded area. To satisfy a requirement, a response plot must remain outside of the
associated exclusion region.
12-9
12 Classical Control Design
If you have Simulink Design Optimization software installed, you can use response
optimization techniques to find a compensator that meets your specified design
requirements. For examples of optimization-based control design using design
requirements, see “Optimize LTI System to Meet Frequency-Domain Requirements”
(Simulink Design Optimization) and “Design Optimization-Based PID Controller for
Linearized Simulink Model (GUI)” (Simulink Design Optimization).
For other Control System Designer tuning methods, you can use the specified design
requirements as visual guidelines during the tuning process.
To add a design requirement to a plot, in Control System Designer, right-click the plot,
and select Design Requirements > New.
12-10
Design Requirements
In the New Design Requirement dialog box, in the Design requirement type drop-down
list, select the type of requirement to add. You can select any valid requirement for the
associated plot type.
To create the specified requirement and add it to the plot, click OK.
When using optimization-based tuning, you can add design requirements from the
Response Optimization dialog box.
To do so, on the Design Requirements tab, click Add new design requirement.
12-11
12 Classical Control Design
In the New Design Requirement dialog box, select a Design requirement type from the
drop-down list.
In the Requirement for response drop-down list, specify the response to which to apply
the design requirement. You can select any response in Data Browser.
To create the specified design requirement, click OK. In the Response Optimization dialog
box, on the Design Requirements tab, the new requirement is added to the table.
The app also adds the design requirement to a corresponding editor or analysis plot. The
plot type used depends on the selected design requirement type.
Otherwise, if the requirement is for a different plot type, the requirement is added to an
appropriate analysis plot. For example, a Step requirement bound is added to a new
step analysis plot.
12-12
Design Requirements
In the Edit Design Requirement dialog box, in the Design requirement drop-down list,
select a design requirement to edit. You can select any existing design requirement from
the current plot.
You can also interactively adjust design requirements by dragging the edges or vertices of
the shaded exclusion region in the associated plot.
12-13
12 Classical Control Design
Specifying a settling time for a continuous-time system adds a vertical boundary line to
the root locus or pole-zero plot. This line represents pole locations associated with the
specified settling time. This boundary is exact for a second-order system with no zeros.
For higher order systems, the boundary is an approximation based on second-order
dominant systems.
To satisfy this requirement, your system poles must be to the left of the boundary line.
For a discrete-time system, the design requirement boundary is a curved line centered on
the origin. In this case, your system poles must be within the boundary line to satisfy the
requirement.
Percent Overshoot
Specifying percent overshoot for a continuous-time system adds two rays to the plot that
start at the origin. These rays are the locus of poles associated with the specified
overshoot value. In the discrete-time case, the design requirement adds two curves
originating at (1,0) and meeting on the real axis in the left-hand plane.
12-14
Design Requirements
Note The percent overshoot (p.o.) design requirement can be expressed in terms of the
damping ratio, ζ:
Ê pz ˆ
p.o. = 100 exp Á - ˜
Á 1- z 2 ˜
Ë ¯
Damping Ratio
Specifying a damping ratio for a continuous-time system adds two rays to the plot that
start at the origin. These rays are the locus of poles associated with the specified
overshoot value. This boundary is exact for a second-order system and, for higher order
systems, is an approximation based on second-order dominant systems.
To meet this requirement, your system poles must be to the left of the boundary lines.
For discrete-time systems, the design requirement adds two curves originating at (1,0)
and meeting on the real axis in the left-hand plane. In this case, your system poles must
be within the boundary curves.
Natural Frequency
Specifying a natural frequency bound adds a semicircle to the plot that is centered
around the origin. The radius of the semicircle equals the natural frequency.
If you specify a natural frequency lower bound, the system poles must remain outside this
semicircle. If you specify a natural frequency upper bound, the system poles must remain
within this semicircle.
Region Constraint
To specify a region constraint, define two or more vertices of a piece-wise linear boundary
line. For each vertex, specify Real and Imaginary components. This requirement adds a
shaded exclusion region on one side of the boundary line. To switch the exclusion region
to the opposite side of the boundary, in the response plot, right-click the requirement, and
select Flip.
To satisfy this requirement, your system poles must be outside of the exclusion region.
12-15
12 Classical Control Design
You can specify upper gain limits for both open-loop and closed-loop Bode responses.
A gain limit consists of one or more line segments. For the start and end points of each
segment, specify a frequency, Freq, and magnitude, Mag. You can also specify the slope
of the line segment in dB/decade. When you change the slope, the magnitude for the end
point updates.
If you are using optimization-based tuning, you can assign a tuning Weight to each
segment to indicate their relative importance.
In the Type drop-down list you can select whether to constrain the magnitude to be above
or below the specified boundary.
You can specify lower gain limits in the same way as upper gain limits.
You can specify a lower bound for the gain margin, the phase margin, or both. The
specified bounds appear in text on the Bode magnitude plot.
Note Gain and phase margin requirements are only applicable to open-loop Bode
diagrams.
12-16
Design Requirements
Gain Margin
Specify a minimum gain margin value. Graphically, Control System Designer displays
this requirement as a region of exclusion along the -180 degree open-loop phase axis.
Specify a minimum closed-loop peak gain value. The specified dB value can be positive or
negative. The design requirement follows the curves of the Nichols plot grid. As a best
practice, have the grid on when using a closed-loop peak gain requirement.
Display Location
When editing a phase margin, gain margin, or closed-loop peak gain requirement, you can
specify the display location as -180 ± k360 degrees, where k is an integer value.
12-17
12 Classical Control Design
If you enter an invalid location, the closest valid location is selected. While displayed
graphically at only one location, these requirements apply regardless of actual phase; that
is, they are applied for all values of k.
You can specify upper time response bounds for both step and impulse responses.
A time-response bound consists of one or more line segments. For the start and end
points of each segment, specify a Time and Amplitude value. You can also specify the
slope of the line segment. When you change the slope, the amplitude for the end point
updates.
12-18
Design Requirements
If you are using optimization-based tuning, you can assign a tuning Weight to each
segment to indicate its relative importance.
In the Type drop-down list, you can select whether to constrain the response to be above
or below the specified boundary.
You can specify lower time response bounds for both step and impulse responses in the
same way as upper gain limits.
For a step response plot, you can also specify a step response bound design requirement.
To define a step response bound requirement, specify the following step response
parameters:
In Control System Designer, step response plots always use an Initial value and a
Step time of 0
12-19
12 Classical Control Design
See Also
More About
• “Optimize LTI System to Meet Frequency-Domain Requirements” (Simulink Design
Optimization)
• “Design Optimization-Based PID Controller for Linearized Simulink Model (GUI)”
(Simulink Design Optimization)
12-20
Feedback Control Architectures
12-21
12 Classical Control Design
If your control application does not match one of the supported control architectures, you
can use block diagram algebra to convert your system to match an architecture. For an
example of such an application, see “Design Multiloop Control System” on page 12-24.
Note If you are unable to match your application to one of the supported control
architectures, consider using the Control System Tuner app to design your control
system.
12-22
See Also
See Also
Control System Designer | sisoinit
More About
• “Design Multiloop Control System” on page 12-24
12-23
12 Classical Control Design
The typical workflow is to tune the compensator for the inner loop first, by isolating the
inner loop from the rest of the control system. Once the inner loop is satisfactorily tuned,
tune the outer loop to achieve your desired closed-loop response.
System Model
For this example develop a position control system for a DC motor. A single-loop angular
velocity controller is designed in “Bode Diagram Design” on page 12-49. To design an
angular position controller, add an outer loop that contains an integrator.
12-24
Design Multiloop Control System
B = [1/L; 0];
C = [0 1];
D = [0];
sys_dc = ss(A,B,C,D);
Design Objectives
The design objective is to minimize the closed-loop step response settling time, while
maintaining an inner-loop phase margin of at least 65 degrees with maximum bandwidth:
Control System Designer has six possible control architectures from which you can
choose. For more information on these architectures, see “Feedback Control
Architectures” on page 12-21.
For this example use Configuration 4, which has an inner and outer control loop.
Currently, the control system structure does not match Configuration 4. However, using
block diagram algebra, you can modify the system model by adding:
12-25
12 Classical Control Design
At the MATLAB command line, add the integrator to the motor plant model.
plant = sys_dc*tf(1,[1,0]);
Create an initial model of the inner-loop compensator that contains the feedback
differentiator.
Cdiff = tf('s');
controlSystemDesigner
In Control System Designer, on the Control System tab, click Edit Architecture.
In the Edit Architecture dialog box, under Select Control Architecture, click the fourth
architecture.
12-26
Design Multiloop Control System
Import the plant and controller models from the MATLAB workspace.
Click OK.
The app updates the control architecture and imports the specified models for the motor
plant and the inner-loop controller.
12-27
12 Classical Control Design
• Bode Editor for LoopTransfer_C1 — Open-loop Bode Editor for the outer loop
• Root Locus Editor for LoopTransfer_C1 — Open-loop Root Locus Editor for the
outer loop
• Bode Editor for LoopTransfer_C2 — Open-loop Bode Editor for the inner loop
• Root Locus Editor for LoopTransfer_C2 — Open-loop root Locus Editor for the
inner loop
• IOTransfer_r2y: step — Overall closed-loop step response from input r to output y
For this example, close the Bode Editor for LoopTransfer_C1 and Root Locus Editor
for LoopTransfer_C2 plots.
Since the inner loop is tuned first, configure the plots to view just the inner-loop Bode
editor plot. On the View tab, click Single, and click Bode Editor for LoopTransfer_C2.
To isolate the inner loop from the rest of the control system architecture, add a loop
opening to the open-loop response of the inner loop. In the Data Browser, right-click
LoopTransfer_C2, and select Open Selection.
To add a loop opening at the output of outer-loop compensator, C1, in the Open-Loop
Transfer Function dialog box, click Add loop opening location to list. Then, select
uC1.
12-28
Design Multiloop Control System
Click OK.
The app adds a loop opening at the selected location. This opening removes the effect of
the outer control loop on the open-loop transfer function of the inner loop.
The Bode Editor response plot updates to reflect the new open-loop transfer function.
12-29
12 Classical Control Design
To increase the bandwidth of the inner loop, increase the gain of compensator C2.
In the Bode Editor plot, drag the magnitude response upward until the phase margin is
65 degrees. This corresponds to a compensator gain of 107. Increasing the gain further
reduces the phase margin below 65 degrees.
12-30
Design Multiloop Control System
Alternatively, you can adjust the gain value using the compensator editor. For more
information, see “Edit Compensator Dynamics” on page 12-93.
With the inner loop tuned, you can now tune the outer loop to reduce the closed-loop
settling time.
12-31
12 Classical Control Design
In Control System Designer, on the View tab, select Left/Right. Arrange the plots to
display the Root Locus for LoopTransfer_C1 and IOTransfer_r2y_step plots
simultaneously.
To view the current settling time, right-click in the step response plot and select
Characteristics > Settling Time.
In the Root Locus Editor, increase the gain of compensator C1. As the gain increases,
the complex pole pair moves toward a slower time constant and the real pole moves
12-32
See Also
toward a faster time constant. A gain of 600 produces a good compromise between rise
time and settling time.
With a closed-loop settling time below 0.8 seconds and an inner-loop phase margin of 65
degrees, the design satisfies the design requirements.
See Also
Control System Designer
12-33
12 Classical Control Design
More About
• “Feedback Control Architectures” on page 12-21
12-34
Multimodel Control Design
Any controller you design for such a system must satisfy the design requirements for all
potential system dynamics.
Model Arrays
In Control System Designer, you can specify multiple models for any plant or sensor in
the current control architecture using an array of LTI models (see “Model Arrays” on
page 2-100). If you specify model arrays for more than one plant or sensor, the lengths of
the arrays must match.
• Create multiple LTI models using the tf, ss, zpk, or frd commands.
12-35
12 Classical Control Design
To import models as arrays, you can pass them as input arguments when opening Control
System Designer from the MATLAB command line. For more information, see Control
System Designer.
You can also import model arrays into Control System Designer when configuring the
control architecture. In the Edit Architecture dialog box:
• In the Value text box, specify the name of an LTI model from the MATLAB workspace.
• To import block data from the MATLAB workspace or from a MAT-file in your current
12-36
Multimodel Control Design
Nominal Model
What Is a Nominal Model?
The nominal model is a representative model in the array of LTI models that you use to
design the controller in Control System Designer. Use the editor and analysis plots to
visualize and analyze the effect of the controller on the remaining plants in the array.
You can select any model in the array as your nominal model. For example, you can
choose a model that:
12-37
12 Classical Control Design
Tip You can plot and analyze the open-loop dynamics of the system on a Bode plot to
determine which model to choose as nominal.
To select a nominal model from the array of LTI models, in Control System Designer,
click Multimodel Configuration. Then, in the Multimodel Configuration dialog box,
select a Nominal model index. The default index is 1.
For each plant or sensor that is defined as a model array, the app selects the model at the
specified index as the nominal model. Otherwise, the app uses scalar expansion to apply
the single LTI model for all model indices.
if G and H are both three-element arrays and the nominal model index is 2, the software
uses the second element in both the arrays to compute the nominal model:
Nominal Model
r y
2
12-38
Multimodel Control Design
CG2
T=
1 + CG2 H2
The app also computes and plots the responses showing the effect of C on the remaining
pairs of plant and sensor models — G1H1 and G3H3.
If only G is an array of LTI models, and the specified nominal model is 2, then the control
architecture for nominal response is:
Nominal Model
r y
2
CG2
T=
1 + CG2 H
The app also computes and plots the responses showing the effect of C on the remaining
pairs of plant and sensor model — G1H and G3H.
Frequency Grid
The frequency response of a system is computed at a series of frequency values, called a
frequency grid. By default, Control System Designer computes a logarithmically equally
spaced grid based on the dynamic range of each model in the array.
• The automatic grid has more points than you require. To improve computational
efficiency, specify a less dense grid spacing.
• The automatic grid is not sufficiently dense within a particular frequency range. For
example, if the response does not capture the resonant peak dynamics of an
underdamped system, specify a more dense grid around the corner frequency.
12-39
12 Classical Control Design
• You are only interested in the response within specific frequency ranges. To improve
computational efficiency, specify a grid that covers only the frequency ranges of
interest.
Note Modifying the frequency grid does not affect the frequency response computation
for the nominal model. The app always uses the Auto select option to compute the
nominal model frequency response.
H1 = tf(1,[1/0.1,1]);
H2 = tf(1,[1/0.15,1]);
H3 = tf(1,[1/0.2,1]);
H = stack(1,H1,H2,H3);
3 Open Control System Designer
Open Control System Designer, and import the plant and sensor model arrays.
controlSystemDesigner(G,1,H)
12-40
Multimodel Control Design
The app opens and imports the plant and sensor model arrays.
4 Configure Analysis Plot
To view the closed-loop step response in a larger plot, in Control System Designer,
on the View tab, click Single. Then, click IOTransfer_r2y: step.
12-41
12 Classical Control Design
By default the step response shows only the nominal response. To display the
individual responses for the other model indices, right-click the plot area, and select
Multimodel Configuration > Individual Responses.
12-42
Multimodel Control Design
Note To view an envelope of all model responses, right-click the plot area, and select
Multimodel Configuration > Bounds
The plot updates to display the responses for the other models.
12-43
12 Classical Control Design
Click Close.
12-44
Multimodel Control Design
To design a compensator using the nominal model, you can use any of the supported
“Control System Designer Tuning Methods” on page 12-4.
12-45
12 Classical Control Design
For this example, use the Compensator Editor to manually specify the compensator
dynamics. Add an integrator to the compensator and set the compensator gain to
0.4. For more information, see “Edit Compensator Dynamics” on page 12-93.
7 Analyze Results
The tuned controller produces a step response with minimal overshoot for the
nominal models and a worst-case overshoot less than 10%.
12-46
See Also
See Also
Control System Designer
Related Examples
• “Model Arrays” on page 2-100
12-47
12 Classical Control Design
12-48
Bode Diagram Design
To adjust the loop shape, you can add poles and zeros to your compensator and adjust
their values directly in the Bode Editor, or you can use the Compensator Editor. For
more information, see “Edit Compensator Dynamics” on page 12-93.
For information on all of the tuning methods available in Control System Designer, see
“Control System Designer Tuning Methods” on page 12-4.
The transfer function of the DC motor plant, as described in “SISO Example: The DC
Motor”, is:
1.5
G=
2
s + 14 s + 40 .02
12-49
12 Classical Control Design
At the MATLAB command line, create a transfer function model of the plant, and open
Control System Designer in the Bode Editor configuration.
G = tf(1.5,[1 14 40.02]);
controlSystemDesigner('bode',G);
The app opens and imports G as the plant model for the default control architecture,
Configuration 1.
• Open-loop Bode Editor for the LoopTransfer_C response. This response is the
open-loop transfer function GC, where C is the compensator and G is the plant.
• Step Response for the IOTransfer_r2y response. This response is the input-output
transfer function for the overall closed-loop system.
Tip To open the open-loop Bode Editor when Control System Designer is already
open, on the Control System tab, in the Tuning Methods drop-down list, select Bode
Editor. In the Select Response to Edit dialog box, select an existing response to plot, or
create a New Open-Loop Response.
To view the open-loop frequency response and closed-loop step response simultaneously,
on the Views tab, click Left/Right.
The app displays the Bode Editor and Step Response plots side-by-side.
12-50
Bode Diagram Design
Adjust Bandwidth
Since the design requires a rise time less than 0.5 seconds, set the open-loop DC
crossover frequency to about 3 rad/s. To a first-order approximation, this crossover
frequency corresponds to a time constant of 0.33 seconds.
To make the crossover easier to see, turn on the plot grid. Right-click the Bode Editor
plot area, and select Grid. The app adds a grid to the Bode response plots.
To adjust the crossover frequency increase the compensator gain. In the Bode Editor
plot, in the Magnitude response plot, drag the response upward. Doing so increases the
gain of the compensator.
12-51
12 Classical Control Design
As you drag the magnitude plot, the app computes the compensator gain and updates the
response plots.
Drag the magnitude response upward until the crossover frequency is about 3 rad/s.
12-52
Bode Diagram Design
To add the rise time to the Step Response plot, right-click the plot area, and select
Characteristics > Rise Time.
To view the rise time, move the cursor over the rise time indicator.
12-53
12 Classical Control Design
The rise time is around 0.23 seconds, which satisfies the design requirements.
Similarly, to add the peak response to the Step Response plot, right-click the plot area,
and select Characteristics > Peak Response.
12-54
Bode Diagram Design
To meet the 5% steady-state error requirement, eliminate steady-state error from the
closed-loop step response by adding an integrator to your compensator. In the Bode
Editor right-click in the plot area, and select Add Pole/Zero > Integrator.
12-55
12 Classical Control Design
To return the crossover frequency to around 3 rad/s, increase the compensator gain
further. Right-click the Bode Editor plot area, and select Edit Compensator.
In the Compensator Editor dialog box, in the Compensator section, specify a gain of 99,
and press Enter.
12-56
Bode Diagram Design
The rise time is around 0.4 seconds, which satisfies the design requirements. However,
the peak overshoot is around 32%. A compensator consisting of a gain and an integrator
is not sufficient to meet the design requirements. Therefore, the compensator requires
additional dynamics.
In the Bode Editor, review the gain margin and phase margin for the current
compensator design. The design requires a gain margin greater than 20 dB and phase
12-57
12 Classical Control Design
margin greater than 40 degrees. The current design does not meet either of these
requirements.
In the Bode Editor, right-click and select Add Pole/Zero > Lead.
To specify the location of the lead network pole, click on the magnitude response. The app
adds a real pole (red X) and real zero (red O) to the compensator and to the Bode Editor
plot.
In the Bode Editor, drag the pole and zero to change their locations. As you drag them,
the app updates the pole/zero values and updates the response plots.
To decrease the magnitude of a pole or zero, drag it towards the left. Since the pole and
zero are on the negative real axis, dragging them to the left moves them closer to the
origin in the complex plane.
12-58
Bode Diagram Design
Tip As you drag a pole or zero, the app displays the new value in the status bar, on the
right side.
As an initial estimate, drag the zero to a location around -7 and the pole to a location
around -11.
12-59
12 Classical Control Design
The phase margin meets the design requirements; however, the gain margin is still too
low.
In the Compensator Editor dialog box, in the Dynamics section, click the Lead row.
12-60
Bode Diagram Design
In the Edit Selected Dynamics section, in the Real Zero text box, specify a location of
-4.3, and press Enter. This value is near the slowest (left-most) pole of the DC motor
plant.
In the Real Pole text box, specify a value of -28, and press Enter.
When you modify a lead network parameters, the Compensator and response plots
update automatically.
In the app, in the Bode Editor, the gain margin of 20.5 just meets the design
requirement.
12-61
12 Classical Control Design
To add robustness to the system, in the Compensator Editor dialog box, decrease the
compensator gain to 84.5, and press Enter. The gain margin increases to 21.8, and the
response plots update.
In Control System Designer, in the response plots, compare the system performance to
the design requirements. The system performance characteristics are:
12-62
See Also
• Overshoot is 3.39%.
• Gain margin is 21.8 dB.
• Phase margin is 65.6 degrees.
See Also
Control System Designer | bodeplot
More About
• “Edit Compensator Dynamics” on page 12-93
• “Control System Designer Tuning Methods” on page 12-4
• “Root Locus Design” on page 12-64
• “Nichols Plot Design” on page 12-80
12-63
12 Classical Control Design
As the open-loop gain, k, of a control system varies over a continuous range of values, the
root locus diagram shows the trajectories of the closed-loop poles of the feedback system.
For example, in the following tracking system:
P(s) is the plant, H(s) is the sensor dynamics, and k is an adjustable scalar gain The
closed-loop poles are the roots of
q ( s) = 1 + kP ( s) H ( s )
The root locus technique consists of plotting the closed-loop pole trajectories in the
complex plane as k varies. You can use this plot to identify the gain value associated with
a desired set of closed-loop poles.
Plant Model
12-64
Root Locus Design
The force on the spool is proportional to the current in the electromagnet coil. As the
spool moves, the valve opens, allowing the high-pressure hydraulic fluid to flow through
the chamber. The moving fluid forces the piston to move in the opposite direction of the
spool. For more information on this model, including the derivation of a linearized model,
see [1].
You can use the input voltage to the electromagnet to control the ram position. When
measurements of the ram position are available, you can use feedback for the ram
position control, as shown in the following, where Gservo represents the
servomechanism:
12-65
12 Classical Control Design
Design Requirements
For this example, tune the compensator, C(s) to meet the following closed-loop step
response requirements:
At the MATLAB command line, load a linearized model of the servomechanism, and open
Control System Designer in the root locus editor configuration.
The app opens and imports Gservo as the plant model for the default control
architecture, Configuration 1.
In Control System Designer, a Root Locus Editor plot and input-output Step
Response open.
To view the open-loop frequency response and closed-loop step response simultaneously,
on the Views tab, click Left/Right.
12-66
Root Locus Design
The app displays Bode Editor and Step Response plots side-by-side.
12-67
12 Classical Control Design
In the closed-loop step response plot, the rise time is around two seconds, which does not
satisfy the design requirements.
To make the root locus diagram easier to read, zoom in. In the Root Locus Editor, right-
click the plot area and select Properties.
In the Property Editor dialog box, on the Limits tab, specify Real Axis and Imaginary
Axis limits from -500 to 500.
12-68
Root Locus Design
Click Close.
To create a faster response, increase the compensator gain. In the Root Locus Editor,
right-click the plot area and select Edit Compensator.
12-69
12 Classical Control Design
In the Root Locus Editor plot, the closed-loop pole locations move to reflect the new
gain value. Also, the Step Response plot updates.
12-70
Root Locus Design
The closed-loop response does not satisfy the settling time requirement and exhibits
unwanted ringing.
Increasing the gain makes the system underdamped and further increases lead to
instability. Therefore, to meet the design requirements, you must specify additional
compensator dynamics. For more information on adding and editing compensator
dynamics, see “Edit Compensator Dynamics” on page 12-93.
12-71
12 Classical Control Design
To add a complex pole pair to the compensator, in the Root Locus Editor, right-click the
plot area and select Add Pole/Zero > Complex Pole. Click the plot area where you want
to add one of the complex poles.
The app adds the complex pole pair to the root locus plot as red X’s, and updates the step
response plot.
In the Root Locus Editor, drag the new poles to locations near –140 ± 260i. As you drag
one pole, the other pole updates automatically.
12-72
Root Locus Design
Tip As you drag a pole or zero, the app displays the new value in the status bar, on the
right side.
12-73
12 Classical Control Design
To add a complex zero pair to your compensator, in the Compensator Editor dialog box,
right-click the Dynamics table, and select Add Pole/Zero > Complex Zero
12-74
Root Locus Design
In the Dynamics table, click the Complex Zero row. Then in the Edit Selected
Dynamics section, specify a Real Part of -170 and an Imaginary Part of 430.
12-75
12 Classical Control Design
The compensator and response plots automatically update to reflect the new zero
locations.
12-76
Root Locus Design
In the Step Response plot, the settling time is around 0.1 seconds, which does not
satisfy the design requirements.
The compensator design process can involve some trial and error. Adjust the compensator
gain, pole locations, and zero locations until you meet the design criteria.
One possible compensator design that satisfies the design requirements is:
12-77
12 Classical Control Design
• Compensator gain of 10
• Complex poles at –110 ± 140i
• Complex zeros at –70 ± 270i
In the Compensator Editor dialog box, configure your compensator using these values. In
the Step Response plot, the settling time is around 0.05 seconds.
To verify the exact settling time, right-click the Step Response plot area and select
Characteristics > Settling Time. A settling time indicator appears on the response plot.
12-78
See Also
To view the settling time, move the cursor over the settling time indicator.
The settling time is about 0.043 seconds, which satisfies the design requirements.
References
[1] Clark, R. N. Control System Dynamics, Cambridge University Press, 1996.
See Also
Control System Designer | rlocusplot
More About
• “Edit Compensator Dynamics” on page 12-93
• “Control System Designer Tuning Methods” on page 12-4
• “Bode Diagram Design” on page 12-49
• “Nichols Plot Design” on page 12-80
12-79
12 Classical Control Design
The transfer function of the DC motor plant, as described in “SISO Example: The DC
Motor”, is:
1.5
G=
2
s + 14 s + 40 .02
At the MATLAB command line, create a transfer function model of the plant, and open
Control System Designer in the Nichols Editor configuration.
G = tf(1.5,[1 14 40.02]);
controlSystemDesigner('nichols',G);
12-80
Nichols Plot Design
The app opens and imports G as the plant model for the default control architecture,
Configuration 1.
• Open-loop Nichols Editor for the LoopTransfer_C response. This response is the
open-loop transfer function GC, where C is the compensator and G is the plant.
• Step Response for the IOTransfer_r2y response. This response is the input-output
transfer function for the overall closed-loop system.
Tip To open the open-loop Nichols Editor when Control System Designer is already
open, on the Control System tab, in the Tuning Methods drop-down list, select Nichols
Editor. In the Select Response to Edit dialog box, select an existing response to plot, or
create a New Open-Loop Response.
To view the open-loop frequency response and closed-loop step response simultaneously,
on the Views tab, click Left/Right.
The app displays the Nichols Editor and Step Response plots side-by-side.
Adjust Bandwidth
Since the design requires a rise time less than 0.5 seconds, set the open-loop DC
crossover frequency to about 3 rad/s. To a first-order approximation, this crossover
frequency corresponds to a time constant of 0.33 seconds.
12-81
12 Classical Control Design
To adjust the crossover frequency increase the compensator gain. In the Nichols Editor,
drag the response upward. Doing so increases the gain of the compensator.
As you drag the Nichols plot, the app computes the compensator gain and updates the
response plots.
Drag the Nichols plot upward until the crossover frequency is about 3 rad/s.
12-82
Nichols Plot Design
To add the rise time to the Step Response plot, right-click the plot area, and select
Characteristics > Rise Time.
To view the rise time, move the cursor over the rise time indicator.
12-83
12 Classical Control Design
The rise time is around 0.23 seconds, which satisfies the design requirements.
Similarly, to add the peak response to the Step Response plot, right-click the plot area,
and select Characteristics > Peak Response.
12-84
Nichols Plot Design
To meet the 5% steady-state error requirement, eliminate steady-state error from the
closed-loop step response by adding an integrator to your compensator. In the Nichols
Editor right-click in the plot area, and select Add Pole/Zero > Integrator.
12-85
12 Classical Control Design
To return the crossover frequency to around 3 rad/s, increase the compensator gain
further. Right-click the Nichols Editor plot area, and select Edit Compensator.
In the Compensator Editor dialog box, in the Compensator section, specify a gain of 99,
and press Enter.
12-86
Nichols Plot Design
The rise time is around 0.4 seconds, which satisfies the design requirements. However,
the peak overshoot is around 32%. A compensator consisting of a gain and an integrator
is not sufficient to meet the design requirements. Therefore, the compensator requires
additional dynamics.
In the Nichols Editor, review the gain margin and phase margin for the current
compensator design. The design requires a gain margin greater than 20 dB and phase
12-87
12 Classical Control Design
margin greater than 40 degrees. The current design does not meet either of these
requirements.
In the Nichols Editor, right-click and select Add Pole/Zero > Lead.
To specify the location of the lead network pole, click on the magnitude response. The app
adds a real pole (red X) and real zero (red O) to the compensator and to the Nichols
Editor plot.
In the Nichols Editor, drag the pole and zero to change their locations. As you drag
them, the app updates the pole/zero values and updates the response plots.
To decrease the magnitude of a pole or zero, drag it towards the left. Since the pole and
zero are on the negative real axis, dragging them to the left moves them closer to the
origin in the complex plane.
Tip As you drag a pole or zero, the app displays the new value in the status bar, on the
right side.
12-88
Nichols Plot Design
As an initial estimate, drag the zero to a location around -7 and the pole to a location
around -11.
The phase margin meets the design requirements; however, the gain margin is still too
low.
In the Compensator Editor dialog box, in the Dynamics section, click the Lead row.
12-89
12 Classical Control Design
In the Edit Selected Dynamics section, in the Real Zero text box, specify a location of
-4.3, and press Enter. This value is near the slowest (left-most) pole of the DC motor
plant.
In the Real Pole text box, specify a value of -28, and press Enter.
When you modify a lead network parameters, the Compensator and response plots
update automatically.
In the app, in the Nichols Editor, the gain margin of 20.5 just meets the design
requirement.
12-90
Nichols Plot Design
To add robustness to the system, in the Compensator Editor dialog box, decrease the
compensator gain to 84.5, and press Enter. The gain margin increases to 21.8, and the
response plots update.
In Control System Designer, in the response plots, compare the system performance to
the design requirements. The system performance characteristics are:
12-91
12 Classical Control Design
• Overshoot is 3.39%.
• Gain margin is 21.8 dB.
• Phase margin is 65.6 degrees.
See Also
Control System Designer | nicholsplot
More About
• “Edit Compensator Dynamics” on page 12-93
• “Control System Designer Tuning Methods” on page 12-4
• “Bode Diagram Design” on page 12-49
• “Root Locus Design” on page 12-64
12-92
Edit Compensator Dynamics
You can add dynamics and modify compensator parameters using the Compensator Editor
or using the graphical Bode Editor, Root Locus Editor, or Nichols Editor plots.
Compensator Editor
To open the Compensator Editor dialog box, in Control System Designer, in an editor
plot area, right-click and select Edit Compensator. Alternatively, in the Data Browser,
in the Controllers section, right-click the compensator you want to edit and click Open
Selection.
12-93
12 Classical Control Design
The Compensator Editor displays the transfer function for the currently selected
compensator. You can select a different compensator to edit using the drop-down list. By
default, the compensator transfer function displays in the time constant format. You can
select a different format by changing the corresponding Control System Designer
preference.
In Control System Designer, on the Control System tab, click Preferences. In the
Control System Designer Preferences dialog box, on the Options tab, select a
Compensator Format.
To add poles and zeros to your compensator, in the Compensator Editor, right-click in the
Dynamics table and, under Add Pole/Zero, select the type of pole/zero you want to add.
12-94
Edit Compensator Dynamics
The app adds a pole or zero of the selected type with default parameters.
To edit a pole or zero, in the Dynamics table, click on the pole/zero type you want to edit.
Then, in the Edit Selected Dynamics section, in the text boxes, specify the pole and zero
locations.
To delete poles and zeros, in the Dynamics table, click on the pole/zero type you want to
delete. Then, right-click and select Delete Pole/Zero.
12-95
12 Classical Control Design
To add poles and zeros directly from an editor plot, right-click the plot area and, under
Add Pole/Zero, select the type of pole/zero you want to add. In the editor plot, the app
displays the editable compensator poles and zeros as red X’s and O’s respectively.
In the editor plots, you can drag poles and zeros to adjust their locations. As you drag a
pole or zero, the app displays the new value in the status bar, on the right side.
To delete a pole or zero, right-click the plot area and select Delete Pole/Zero. Then, in
the editor plot, click the pole or zero you want to delete.
12-96
Edit Compensator Dynamics
To configure a lead or lag network for your compensator, use one of the following options:
• Specify the pole and zero locations. Placing the pole and zero further apart increases
the amount of phase angle change.
12-97
12 Classical Control Design
• Specify the maximum amount of phase angle change and the frequency at which this
change occurs. The app automatically computes the pole and zero locations.
When graphically changing pole and zero locations for a lead or lag compensator, in the
editor plot, you can drag the pole and zeros independently.
Notch Filters
If you know that your system has disturbances at a particular frequency, you can add a
notch filter to attenuate the gain of the system at that frequency. The notch filter transfer
function is:
s2 + 2x1wn s + wn2
s2 + 2x2wn s + wn2
where
To configure a notch filter for your compensator, in the Compensator Editor dialog box,
you can specify the: