Computation
Visualization
Programming
For Use with MATLAB
?
Model Predictive Control
Toolbox
User’s Guide
Version 1
Manfred Morari
N. Lawrence Ricker
How to Contact The MathWorks:
508-647-7000 Phone
508-647-7001 Fax
The MathWorks, Inc. Mail
3AppleHillDrive
Natick, MA 01760-2098
http://www.mathworks.com Web
ftp.mathworks.com Anonymous FTP server
comp.soft-sys.matlab Newsgroup
support@mathworks.com Technical support
suggest@mathworks.com Product enhancement suggestions
bugs@mathworks.com Bug reports
doc@mathworks.com Documentation error reports
subscribe@mathworks.com Subscribing user registration
service@mathworks.com Order status, license renewals, passcodes
info@mathworks.com Sales, pricing, and general information
Model Predictive Control Toolbox User’s Guide
COPYRIGHT 1995 - 1998 by The MathWorks, Inc.
The software described in this document is furnished under a license agreement. The software may be used
or copied only under the terms of the license agreement. No part of this manual may be photocopied or repro-
duced in any form without prior written consent from The MathWorks, Inc.
FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation by
or for the federal government of the United States. By accepting delivery of the Program, the government
hereby agrees that this software qualifies as "commercial" computer software within the meaning of FAR
Part 12.212, DFARS Part 227.7202-1, DFARS Part 227.7202-3, DFARS Part 252.227-7013, and DFARS Part
252.227-7014. The terms and conditions of The MathWorks, Inc. Software License Agreement shall pertain
to the government’s use and disclosure of the Program and Documentation, and shall supersede any
conflicting contractual terms or conditions. If this license fails to meet the government’s minimum needs or
is inconsistent in any respect with federal procurement law, the government agrees to return the Program
and Documentation, unused, to MathWorks.
MATLAB, Simulink, Stateflow, Handle Graphics, and Real-Time Workshop are registered trademarks, and
Target Language Compiler is a trademark of The MathWorks, Inc.
Other product or brand names are trademarks or registered trademarks of their respective holders.
Printing History: January 1995 First printing
October 1998 Revised (Online only)
i
Contents
Preface
1
Tutorial
Introduction ......................................... 1-2
TargetAudiencefortheMPCToolbox .................... 1-3
SystemRequirements ................................. 1-3
2
MPC Based on Step Response Models
Step Response Models ................................. 2-2
Model Identification .................................. 2-6
Unconstrained Model Predictive Control .............. 2-11
Closed-Loop Analysis ................................ 2-18
Constrained Model Predictive Control ................. 2-20
Application: Idle Speed Control ....................... 2-22
ProcessDescription .................................. 2-22
ControlProblemFormulation.......................... 2-22
Simulations ........................................ 2-24
Application: Control of a Fluid Catalytic Cracking Unit . 2-31
ProcessDescription .................................. 2-31
ControlProblemFormulation.......................... 2-33
ii Contents
Simulations ........................................ 2-34
StepResponseModel............................... 2-34
AssociatedVariables ............................... 2-36
UnconstrainedControlLaw ......................... 2-36
ConstrainedControlLaw ........................... 2-36
3
MPC Based on State-Space Models
State-Space Models .................................... 3-2
ModFormat ......................................... 3-3
SISO Continuous-Time Transfer Function to Mod Format . . . . 3-3
SISO Discrete-Time Transfer Function to Mod Format . . . . . . 3-6
MIMOTransferFunctionDescriptiontoModFormat ....... 3-7
ContinuousorDiscreteState-SpacetoModFormat ......... 3-9
Identification Toolbox (“Theta”) Format to Mod Format . . . . . . 3-9
CombinationofModelsinModFormat .................. 3-10
ConvertingModFormattoOtherModelFormats .......... 3-10
Unconstrained MPC Using State-Space Models ......... 3-12
State-Space MPC with Constraints .................... 3-20
Application: Paper Machine Headbox Control .......... 3-26
MPCDesignBasedonNominalLinearModel............. 3-27
MPCofNonlinearPlant .............................. 3-38
4
Command Reference
Commands Grouped by Function ....................... 4-2
Index
Preface
Preface
iv
Acknowledgments
The toolbox was developed in cooperation with: Douglas B. Raven and Alex
Zheng
The contributions of the following people are acknowledged: Yaman Arkun,
Nikolaos Bekiaris, Richard D. Braatz, Marc S. Gelormino, Evelio Hernandez,
TylerR.Holcomb,IftikharHuq,SameerM.Jalnapurkar,JayH.Lee,Yusha
Liu, Simone L. Oliveira, Argimiro R. Secchi, and Shwu-Yien Yang
We would like to thank Liz Callanan, Jim Tung and Wes Wang from the
MathWorks for assisting us with the project, and Patricia New who did such
an excellent job putting the manuscript into LATEX.
About the Authors
v
About the Authors
Manfred Morari
Manfred Morari received his diploma from ETH Zurich in 1974 and his Ph.D.
from the University of Minnesota in 1977, both in chemical engineering.
Currently he is the McCollum-Corcoran Professor and Executive Officer for
Control and Dynamical Systems at the California Institute of Technology.
Morari’s research interests are in the areas of process control and design. In
recognition of his numerous contributions, he has received the Donald P.
Eckman Award of the Automatic Control Council, the Allan P. Colburn Award
of the AIChE, the Curtis W. McGraw Research Award of the ASEE, was a Case
Visiting Scholar, the Gulf Visiting Professor at Carnegie Mellon University
and was recently elected to the National Academy of Engineering. Dr. Morari
hasheldappointmentswithExxonR&EandICIandhasconsulted
internationally for a number of major corporations. He has coauthored one
book on Robust Process Control with another on Model Predictive Control in
preparation.
N. Lawrence Ricker
Larry Ricker received his B.S. degree from the University of Michigan in 1970,
and his M.S. and Ph.D. degrees from the University of California, Berkeley, in
1972/78. All are in Chemical Engineering. He is currently Professor of
Chemical Engineering at the University of Washington, Seattle. Dr. Ricker has
over 80 publications in the general area of chemical plant design and operation.
He has been active in Model Predictive Control research and teaching for more
than a decade. For example, he published one of the first nonproprietary
studies of the application of MPC to an industrial process, and is currently
involved in a large-scale MPC application involving more than 40 decision
variables.
Preface
vi
1
Tutorial
1 Tutorial
1-2
Introduction
The Model Predictive Control (MPC) Toolbox is a collection of functions
(commands) developed for the analysis and design of model predictive control
(MPC) systems. Model predictive control was conceived in the 1970s primarily
by industry. Its popularity steadily increased throughout the 1980s. At
present, there is little doubt that it is the most widely used multivariable
control algorithm in the chemical process industries and in other areas. While
MPC is suitable for almost any kind of problem, it displays its main strength
when applied to problems with:
? A large number of manipulated and controlled variables
? Constraints imposed on both the manipulated and controlled variables
? Changing control objectives and/or equipment (sensor/actuator) failure
? Time delays
Some of the popular names associated with model predictive control are
Dynamic Matrix Control (DMC), IDCOM, model algorithmic control, etc. While
these algorithms differ in certain details, the main ideas behind them are very
similar. Indeed, in its basic unconstrained form MPC is closely related to linear
quadratic optimal control. In the constrained case, however, MPC leads to an
optimization problem which is solved on-line in real time at each sampling
interval. MPC takes full advantage of the power available in today’s control
computer hardware.
This software and the accompanying manual are not intended to teach the user
the basic ideas behind MPC. Background material is available in standard
textbooks like those authored by Seborg, Edgar and Mellichamp (1989)
1
,
Deshpande and Ash (1988)
2
and the monograph devoted solely to this topic
authored by Morari and coworkers (Morari et al., 1994)
3
. This section provides
a basic introduction to the main ideas behind MPC and the specific form of
implementation chosen for this toolbox. The algorithms used here are
consistent with those described in the monograph by Morari et al. Indeed, the
software is meant to accompany the monograph and vice versa. The routines
included in the MPC Toolbox fall into two basic categories: routines which use
1. D.E. Seborg, T.F. Edgar, D.A. Mellichamp; Process Dynamics and Control; JohnWiley & Sons,
1989
2. P.B. Deshpande, R.H. Ash; Computer Process Control with Advanced Control Applications,
2nd ed., ISA, 1988
3. M. Morari, C.E. Garcia, J.H. Lee, D.M. Prett; Model Predictive Control; Prentice Hall, 1994
(In the process of being written.)
Introduction
1-3
a step response model description and routines which use a state-space model
description. In addition, simple identification tools are provided for identifying
step response models from plant data. Finally, there are also various
conversion routines which convert between different model formats and
analysis routines which can aid in determining the stability of the
unconstrained system, etc. All MPC Toolbox commands have a built-in usage
display. Any command called with no input arguments results in a brief
description of the command line. For example, typing mpccon at the command
line gives the following:
usage: Kmpc = mpccon(model,ywt,uwt,M,P)
The following sections include several examples. They are available in the
tutorial programs mpctut.m, mpctutid.m, mpctutst.m,andmpctutss.m. You
can copy these demo files from the mpctools/mpcdemos source into a local
directory and examine the effects of modifying some of the commands.
Target Audience for the MPC Toolbox
The package is intended for the classroom and for the practicing engineer. It
can assist in communicating the concepts of MPC to a student in an
introductory control course. At the same time it is sophisticated enough to
allow an engineer in industry to evaluate alternate control strategies in
simulation studies.
System Requirements
The MPC Toolbox assumes the following operating system requirements:
? MATLAB
?
is running on your system.
? If nonlinear systems are to be simulated, Simulink
?
is required for the
functions nlcmpc and nlmpcsim.
? If the theta format from the System Identification Toolbox is to be used to
create models in the MPC mod format (using the MPC Toolbox function,
th2mod), then the System Identification Toolbox function polyform and the
Control System Toolbox function append are required.
1 Tutorial
1-4
The MPC Toolbox analysis and simulation algorithms are numerically
intensive and require approximately 1MB of memory, depending on the
number of inputs and outputs. The available memory on your computer may
limit the size of the systems handled by the MPC Toolbox.
Note: there is a pack command in MATLAB that can help free memory space
by compacting fragmented memory locations. For reasonable response times,
a computer with power equivalent to an 80386 machine is recommended
unless only simple tutorial example problems are of interest.
2
MPC Based on Step
Response Models
2 MPC Based on Step Response Models
2-2
Step Response Models
Step response models are based on the following idea. Assume that the system
is at rest. For a linear time-invariant single-input single-output (SISO) system
let the output change for a unit input change D v be given by
{0, s
1
, s
2
,...,s
n
, s
n
,...}
Here we assume that the system settles exactly after n steps. The step response
{s
1
, s
2
,...,s
n
} constitutes a complete model of the system, which allows us to
compute the system output for any input sequence:
Step response models can be used for both stable and integrating processes. For
an integrating process it is assumed that the slope of the response remains
constant after n steps, i.e.,
s
n
– s
n –1
= s
n +1
– s
n
= s
n +2
– s
n +1
=...
For a multi-input multi-output (MIMO) process with n
v
inputs and n
y
outputs,
one obtains a series of step response coefficient matrices
where s
l,m,i
is the i
th
step response coefficient relating the m
th
input to the l
th
output.
yk() s
i
vD ki–()
i 1=
n
Ge5
s
n
vk n– 1–()+=
S
i
s
11i,,
s
12i,,
… s
1 n
v
i,,
s
21i,,
s
n
y
1 i,,
s
n
y
2 i,,
… s
n
y
n
v
i,,
=
…
Step Response Models
2-3
The MPC Toolbox stores step response models in the following format:
where delt2 isthesamplingtimeandthevectornout indicates if a particular
output is integrating or not:
nout(i) = 0 if output i is integrating.
nout(i) = 1 if output i is stable.
The step response can be obtained directly from identification experiments, or
generated from a continuous or discrete transfer function or state-space model.
For example, if the discrete system description (sampling time T =0.1)is
y(k)=–0.5y(k –1)+v(k –3)
then the transfer function is
plant
S
1
S
2
S
n
nout(1) 0 … 0
nout(2) 0 … 0
nout n
y
() 0 … 0
n
y
0 … 0
delt2 0 … 0
nn
y
n
y
2++()n
v
·
=
…
…… …
gz()
z
3–
10.5z
1–
+
--------------------------=
2 MPC Based on Step Response Models
2-4
The following commands (see mpctut.m) generate the step response model for
this system and plot it:
num = 1;
den = [1 0.5];
delt1 = 0.1;
delay = 2;
g = poly2tfd(num,den,delt1,delay);
% Set up the model in tf format
tfinal = 1.6;
delt2 = delt1;
nout = 1;
plant = tfd2step(tfinal,delt2,nout,g);
% Calculate the step response
plotstep(plant) % Plot the step response
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
u1 step response : y1
TIME
Step Response Models
2-5
Alternatively, we could first generate a state-space description applying the
command tf2ss and then generate the step response with ss2step.Inthis
case, we need to pad the numerator and denominator polynomials to account
for the time delay.
num = [0 0 0 num];
den = [den 0 0];
[phi,gam,c,d] = tf2ss(num,den); % Convert to state-space
plant = ss2step(phi,gam,c,d,tfinal,delt1,delt2,nout);
% Calculate step response
We can get some information on the contents of a matrix in the MPC Toolbox
via the command mpcinfo.Forourexample,mpcinfo(plant) returns:
This is a matrix in MPC Step format.
sampling time = 0.1
number of inputs = 1
number of outputs = 1
number of step response coefficients = 16
All outputs are stable.
2 MPC Based on Step Response Models
2-6
Model Identification
The identification routines available in the MPC Toolbox are designed for
multi-input single-output (MISO) systems. Based on a historical record of the
output y
l
(k) and the inputs v
1
(k); v
2
(k),..., (k),
thestepresponsecoefficients
are estimated. For the estimation of the step response coefficients we write the
SISO model in the form
where
D y(k)=y(k)–y(k –1)
and
h
i
= s
i
–s
i –1
v
n
v
yy
l
y
l
1()
y
l
2()
y
l
3()
=
…
v
v
1
1()v
2
1() … v
n
v
1()
v
1
2()v
2
2() … v
n
v
2()
v
1
3()v
2
3() … v
n
v
3()
=
… …
s
l 11,,
s
l 21,,
… s
l n
v
1,,
s
l 12,,
s
l 22,,
… s
l n
v
2,,
s
l 1 i,,
s
l 2 i,,
… s
l n
v
i,,
… … …
…
D yk() h
i
vD ki–()
i 1=
n
Ge5
=
Model Identification
2-7
h
i
are the impulse response coefficients. This model allows the designer to
present all the input (v) and output(y
l
) information in deviation form, which is
often desirable. If the particular output is integrating, then the model
where
D (D y(k)) = D y(k)–D y(k –1)
D h
i
= h
i
– h
i –1
should be used to estimate D h
i
, and thus h
i
and s
i
are given by
For parameter estimation it is usually recommended to scale all the variables
such that they are the same order of magnitude. This may be done via the MPC
Toolbox functions autosc or scal. Then the data has to be arranged into the
form
Y = XQ
where Y contains all the output information (D y(k) for stable and D (D y(k)) for
integrating outputs) and X all the input information (D v(k)) appropriately
arranged. Q is a vector including all the parameters to be estimated (h
i
for
stable and D h
i
for integrating outputs). This rearrangement of the input and
output information is handled by wrtreg. The parameters Q can be estimated
via multivariable least squares regression (mlr) or partial least squares (plsr).
Finally, the step response model is generated from the impulse response
coefficients via imp2step. The following example (seempctutid) illustrates this
procedure.
DD yk()() Dh
i
vD ki–()
i 1=
n
Ge5
=
h
i
D h
k
k 1=
i
Ge5
=
s
i
h
j
j 1=
i
Ge5
D h
k
k 1=
j
Ge5
j 1=
i
Ge5
==
2 MPC Based on Step Response Models
2-8
Example:
% Load the input and output data. The input and output
% data are generated from the following transfer
% functions and random zero-mean noises.
% TF from input 1 to output 1: g11=5.72exp(-14s)/
% (60s+1)
% TF from input 2 to output 1: g12=1.52exp(-15s)/
% (25s+1)
% Sampling time of 7 minutes was used.
%
% load mlrdat
%
% Determine the standard deviations for input data using
% autosc.
[ax,mx,stdx]=autosc(x);
%
% Scale the input data by their standard deviations only.
mx=0*mx;
sx=scal(x,mx,stdx);
%
% Put the input and output data in a form such that they
% can be used to determine the impulse response
% coefficients. 35 coefficients (n) are used.
n=35;
[xreg,yreg]=wrtreg(sx,y,n);
%
% Determine the impulse response coefficients via mlr.
% By specifying plotopt=2, two plots—plot of predicted
% output and actual output, and plot of the output
% residual (or prediction error)—are produced.
ninput=2; plotopt=2;
[theta,yres]=mlr(xreg,yreg,ninput,plotopt);
Model Identification
2-9
% Scale theta based on the standard deviations used in
% scaling the input.
theta=scal(theta,mx,stdx);
%
% Convert the impulse model to a step model to be used
% in MPC design.
% Sampling time of 7 minutes was used in determining
% the impulse model.
% Number of outputs (1 in this case) must be specified.
nout=1;
delt=7;
model=imp2step(delt,nout,theta);
%
% Plot the step response.
plotstep(model)
0 20 40 60 80 100 120 140 160 180
?6
?4
?2
0
2
4
Actual value (o) versus Predicted Value (+)
Sample Number
Output
0 20 40 60 80 100 120 140 160 180
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
Output Residual or Prediction Error
Sample Number
Residual
2 MPC Based on Step Response Models
2-10
0 50 100 150 200 250
0
1
2
3
4
5
6
u1 step response : y1
TIME
0 50 100 150 200 250
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
u2 step response : y1
TIME
Unconstrained Model Predictive Control
2-11
Unconstrained Model Predictive Control
The MPC control law can be most easily derived by referring to the following
figure.
For any assumed set of present and future control moves D u(k), D u(k +1),
...,D u(k + m – 1) the future behavior of the process outputs y(k +1|k),
y(k +2|k),...,y(k + p|k) can be predicted over a horizon p.Them present and
future control moves (m ￡ p) are computed to minimize a quadratic objective of
the form
Here and are weighting matrices to penalize particular components of
y or u at certain future time intervals. r(k + l) is the (possibly time-varying)
vector of future reference values (setpoints). Though m control moves are
calculated, only the first one (D u(k)) is implemented. At the next sampling
interval, new values of the measured output are obtained, the control horizon
is shifted forward by one step, and the same computations are repeated. The
resulting control law is referred to as “moving horizon” or “receding horizon.”
futurepast
target
u(k)
k+pk
k+1
k+2
horizon
y(k+1|k)
^
y
Projected
Manipulated
Variables
Outputs
min
D uk()…Duk m 1–+()
G
l
y
yk l|k+()rk l+()–[]
2
l=1
p
Ge5
+
G
l
u
D uk l+ 1–()[]
2
l=1
m
Ge5
G
l
y
G
l
u
2 MPC Based on Step Response Models
2-12
The predicted process outputs y(k +1|k),...,y(k + p|k) depend on the current
measurement and assumptions we make about the unmeasured
disturbances and measurement noise affecting the outputs. The MPC Toolbox
assumes that the unmeasured disturbances for each output are steps passing
through a first order lag with time constant tfilter(2,:).
1
For rejecting
measurement noise, the time constant of an exponential filter tfilter(1,:)
can be specified by the user. (It can be shown that this procedure is optimal for
white noise disturbances passed through an integrator and a first order lag,
and white measurement noise). For conventional Dynamic Matrix Control
(DMC) the disturbance time constant is assumed to be zero
(tfilter(2,:) = zeros(1,ny)), i.e., the unmeasured disturbances have the
form of steps, and the noise filter time constant is also set to zero
(tfilter(1,:) = zeros(1,ny)),i.e.,thereisnomeasurementnoisefiltering
for doing the prediction.
Under the stated assumptions, it can be shown that a linear time-invariant
feedback control law results
D u(k)=K
MP C
E
p
(k +1|k)
where E
p
(k +1|k) is the vector of predicted future errors over the horizon p
which would result if all present and future manipulated variable moves were
equal to zero D u(k)=D u(k +1)=...=0.
For open-loop stable plants, nominal stability of the closed-loop system
depends only on K
MP C
whichinturnisaffectedbythehorizonp, the number
of moves m and the weighting matrices and . No precise conditions on m,
p, and exist which guarantee closed-loop stability. In general,
decreasing m relative to p makes the control action less aggressive and tends
to stabilize a system. For p = 1, nominal stability of the closed-loop system is
guaranteed for any finite m, and time-invariant input and output weights.
More commonly, is used as a tuning parameter. Increasing always has
the effect of making the control action less aggressive.
The noise filter time constanttfilter(1,:) and the disturbance time constant
tfilter(2,:) do not affect closed-loop stability or the response of the system
to setpoint changes or measured disturbances. They do, however, affect the
robustness and the response to unmeasured disturbances.
1. See cmpc in the “Reference” section for details on how to specify tfilter.
y? k()
G
l
y
G
l
u
G
l
y
G
l
u
G
l
u
G
l
u
Unconstrained Model Predictive Control
2-13
Increasing the noise filter time constant makes the system more robust and the
unmeasured disturbance response more sluggish. Increasing the disturbance
time constant increases the lead in the loop, making the control action
somewhat more aggressive, and is recommended for disturbances which have
a slow effect on the output.
All controllers designed with the MPC Toolbox track steps asymptotically
error-free (Type 1). If the unmeasured disturbance model or the system itself
is integrating, ramps are also tracked error-free (Type 2).
Example: (see mpctutst.m)
% Plant transfer function: g=5.72exp(-14s)/(60s+1)
% Disturbance transfer function: gd=1.52exp(-15s)/
% (25s+1)
%
% Build the step response models for a sampling period
% of 7.
delt1=0;
delay1=14;
num1=5.72;
den1=[60 1];
g=poly2tfd(num1,den1,delt1,delay1);
tfinal=245;
delt2=7;
nout1=1;
plant=tfd2step(tfinal,delt2,nout1,g);
delay2=15;
num2=1.52;
den2=[25 1];
gd=poly2tfd(num2,den2,delt1,delay2);
delt2=7;
nout2=1;
dplant=tfd2step(tfinal,delt2,nout2,gd);
%
% Calculate the MPC controller gain matrix for
% No plant/model mismatch,
% Output Weight=1, Input Weight=0
% Input Horizon=5, Output Horizon=20
model=plant;
ywt=1; uwt=0;
2 MPC Based on Step Response Models
2-14
M=5; P=20;
Kmpc1=mpccon(model,ywt,uwt,M,P);
%
% Simulate and plot response for unmeasured and measured
% step disturbance through dplant.
tend=245;
r=[ ]; usat=[ ]; tfilter=[ ];
dmodel=[ ];
dstep=1;
[y1,u1]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
dmodel=dplant; % measured disturbance
[y2,u2]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
plotall([y1,y2],[u1,u2],delt2);
pause; % Perfect rejection for measured disturbance case.
%
% Calculate a new MPC controller gain matrix for
% No plant/model mismatch,
% Output Weight=1, Input Weight=10
0 50 100 150 200 250
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 50 100 150 200 250
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
Manipulated Variables
Time
Unconstrained Model Predictive Control
2-15
% Input Horizon=5, Output Horizon=20
model=plant;
ywt=1; uwt=10;
M=5; P=20;
mpc2=mpccon(model,ywt,uwt,M,P);
%
% Simulate and plot response for unmeasured and measured
% step disturbance through dplant.
tend=245;
r=[ ]; usat=[ ]; tfilter=[ ];
dmodel=[ ];
dstep=1;
[y3,u3]=mpcsim(plant,model,Kmpc2,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
dmodel=dplant; % measured disturbance
[y4,u4]=mpcsim(plant,model,Kmpc2,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
plotall([y3,y4],[u3,u4],delt2);
pause;
0 50 100 150 200 250
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 50 100 150 200 250
?0.35
?0.3
?0.25
?0.2
?0.15
?0.1
?0.05
0
Manipulated Variables
Time
2 MPC Based on Step Response Models
2-16
%
% Simulate and plot response for unmeasured
% step disturbance through dplant with uwt=0,
% with and without noise filtering.
tend=245;
r=[ ]; usat=[ ]; dmodel=[ ];
tfilter=[ ];
dstep=1;
[y5,u5]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
tfilter=20; % noise filtering time constant=20
[y6,u6]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
plotall([y5,y6],[u5,u6],delt2);
pause;
0 50 100 150 200 250
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 50 100 150 200 250
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
Manipulated Variables
Time
Unconstrained Model Predictive Control
2-17
%
% Simulate and plot response for unmeasured
% step disturbance through dplant with uwt=0,
% with and without unmeasured disturbance time constant
% being specified.
tend=245;
r=[ ]; usat=[ ]; dmodel=[ ];
tfilter=[ ];
dstep=1;
[y7,u7]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
tfilter=[0;25]; % unmeasured disturbance time constant=25
[y8,u8]=mpcsim(plant,model,Kmpc1,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
plotall([y7,y8],[u7,u8],delt2);
pause;
0 50 100 150 200 250
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 50 100 150 200 250
?1.5
?1
?0.5
0
Manipulated Variables
Time
2 MPC Based on Step Response Models
2-18
Closed-Loop Analysis
Apart from simulation, other tools are available in the MPC Toolbox to analyze
the stability and performance of a closed-loop system. We can obtain the
state-space description of the closed-loop system with the command mpccl and
then determine the pole locations with smpcpole.
Example: (mpctutst.m)
% Construct a closed-loop system for no disturbances
% and uwt = 0. Determine the poles of the system.
clmod = mpccl(plant,model,Kmpc1);
poles = smpcpole(clmod);
maxpole = max(poles)
Result is: maxpole = 1.0
The closed-loop system is stable if all the poles are inside or on the unit-circle.
Furthermore we can examine the frequency response of the closed-loop system.
For multivariable systems, singular values as a function of frequency can be
obtained using svdfrsp.
Example: (mpctutst.m)
% Calculate and plot the frequency response of the
% sensitivity and complementary sensitivity functions.
freq = [ -3,0,30];
ny = 1;
in = [1:ny]; % input is r for comp. sensitivity
out = [1:ny]; % output is yp for comp. sensitivity
[frsp,eyefrsp] = mod2frsp(clmod,freq,out,in);
plotfrsp(eyefrsp); % sensitivity
pause;
plotfrsp(frsp); % complementary sensitivity
pause; % Magnitude = 1 for all frequencies chosen.
Closed-Loop Analysis
2-19
10
?3
10
?2
10
?1
10
0
10
?2
10
?1
10
0
10
1
Frequency (radians/time)
Log Magnitude
BODE PLOTS
10
?3
10
?2
10
?1
10
0
?100
?50
0
50
100
Frequency (radians/time)
Phase (degrees)
10
?3
10
?2
10
?1
10
0
10
0
10
1
Frequency (radians/time)
Log Magnitude
BODE PLOTS
10
?3
10
?2
10
?1
10
0
?800
?600
?400
?200
0
Frequency (radians/time)
Phase (degrees)
2 MPC Based on Step Response Models
2-20
Constrained Model Predictive Control
The control action can also be computed subject to hard constraints on the
manipulated variables and the outputs.
Manipulated variable constraints:
u
min
(l) ￡ u(k +l) ￡ u
max
(l)
Manipulated variable rate constraints:
|D u(k + l)| ￡D u
max
(l)
Output variable constraints:
y
min
(l) ￡ y(k +l|k) ￡ y
max
(l)
When hard constraints of this form are imposed, a quadratic program has to be
solved at each time step to determine the control action and the resulting
control law is generally nonlinear. The performance of such a control system
has to be evaluated via simulation.
Example: (mpctutst.m)
% Simulate and plot response for unmeasured step
% disturbance through dplant with and without
% input constraints.
% No plant/model mismatch,
% Output Weight = 1, Input Weight = 0
% Input Horizon = 5, Output Horizon = 20
% Minimum Constraint on Input = -0.4
% Maximum Constraint on Input = inf
% Delta Constraint on Input = 0.1
model = plant;
ywt = 1; uwt = 0;
M = 5; P = 20;
tend = 245;
r = 0;
ulim = [ ];
ylim = [ ]; tfilter = [ ]; dmodel = [ ];
dstep = 1;
[y9,u9] = cmpc(plant,model,ywt,uwt,M,P,tend,r,...
ulim,ylim,tfilter,dplant,dmodel,dstep);
ulim = [ -0.4, inf, 0.1]; % impose constraints
Constrained Model Predictive Control
2-21
[y10,u10] = cmpc(plant,model,ywt,uwt,M,P,tend,r,...
ulim,ylim,tfilter,dplant,dmodel,dstep);
plotall([y9,y10],[u9,u10],delt2);
0 50 100 150 200 250
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 50 100 150 200 250
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
0.1
Manipulated Variables
Time
2 MPC Based on Step Response Models
2-22
Application: Idle Speed Control
Process Description
An idle speed control
2
system should maintain the desired engine rpm with
automatic transmission in neutral or drive. Despite sudden load changes due
to the actions of air conditioning, power steering, etc., the control system
should maintain smooth stable operation. Because of varying operating
conditions and engine-to-engine variability inevitable in mass production, the
system dynamics may change. The controller must be designed to be robust
with respect to these changes. Two control inputs, bypass valve (or throttle)
opening and spark advance, are used to maintain the engine rpm at a desired
idle speed level. For safe operation, spark advance should not change by more
than 20 degrees. Also, spark advance should be at the design point at
steady-state for fuel economy. Thus, spark advance is viewed both as a
manipulated input and a controlled output.
Control Problem Formulation
Here we consider two different operating conditions (transmission in neutral
and drive positions) and the models for the two plants are taken from the paper
by Hrovat and Bodenheimer.
3
Thegoalistodesignamodelpredictive
controller such that the closed loop performance at both operating conditions is
good in the presence of the input constraint specified above. There is no
synthesis method available which systematically generates a controller design
which guarantees robust performance (or even just robust stability) in the
presence of constraints. Thus, we must rely on both design and simulation tools
to determine achievable performance objectives when there are both
constraints and robustness requirements. The toolbox helps us toward this
objective.
Consider the following system:
2. More detail on the problem formulation can be found in the paper by Williams et al., “Idle
Speed Control Design Using an H
￥
Approach,” Proceedings of American Control Conference,
1989, pages 1950–1956.
3. D. Hrovat and B. Bodenheimer, “Robust Automotive Idle Speed Control Design Based on
μ-Synthesis,” Proceedings of American Control Conference, 1993, pages 1778–1783.
y
1
y
2
G
11
G
21
01
u
1
u
2
G
d
0
w+=
Application: Idle Speed Control
2-23
where y
1
is engine rpm, y
2
and u
2
are spark advance, u
1
is bypass valve, w is
torque load (unmeasured disturbance), and G
11
, G
21
and Gd are the
corresponding transfer functions. After scaling, the constraints on spark
advance become ±0.7, i.e., |u
2
| ￡ 0.7.
Plant #1 corresponds to operation in drive at 800 rpm and a load of 30 Nm and
the transfer functions are given by
Plant #2 corresponds to operation at 800 rpm in neutral with zero load and the
transfer functions are given by
The goal is to design a model predictive controller such that the closed-loop
performance is good for plants #1 and #2 when subjected to an unmeasured
torque load disturbance.
G
11
9.62e
0.16s–
s
2
2.4s 5.05++
----------------------------------------=
G
21
15.9 s( 3 )e+
0.04s–
s
2
2.4s 5.05++
-----------------------------------------------=
G
d
19.1– s 3+()
s
2
2.4s 5.05++
-----------------------------------------=
G
11
20.5e
0.16s–
s
2
2.2s 12.8++
-----------------------------------------=
G
21
47.6 s( 3.5 )e+
0.04s–
s
2
2.2s 12.8++
----------------------------------------------------=
G
d
19.1– s 3.5+()
s
2
2.2s 12.8++
----------------------------------------=
2 MPC Based on Step Response Models
2-24
Simulations
Since the toolbox handles only discrete-time systems, the models are
discretizedusingasamplingtimeof0.1.Weapproximateeachofthediscrete
transfer functions with 40 step response coefficients. The function cmpc is used
to generate the controller and to simulate the closed-loop response because it
determines optimal changes of the manipulated variables subject to
constraints. For comparison (Simulation # 4), we also use the functions mpccon
for controller design and mpcsim for simulating closed-loop responses. On-line
computations are simpler, but the resulting controller is linear and the
constraints are not handled in an optimal fashion. The following additional
functions from the toolbox are also used: tfd2step and plotall.TheMATLAB
code for the following simulations can be found in the file idlectr.m in the
directory mpcdemos.
Figure 2-1 Responses to a Unit Torque Disturbance for Plant #1 (no model/
plant mismatch)
Simulation #1. No model/plant mismatch. The following parameters are used:
M = 10, P = inf, ywt = [5 1], uwt = [0.5 0.5],
tfilter = [ ]
The larger weight on the first output (engine rpm) is to emphasize that
controlling engine rpm is more important than controlling spark advance.
Figure 2-1 and Figure 2-2 show the closed-loop response for a unit step torque
123456789100
5
0
-5
123456789100
5
0
-5
10
Outputs
Time
Manipulated Variables
Time
y2
y1
u1
u2
Application: Idle Speed Control
2-25
load change. No model/plant mismatch is introduced, i.e., we use Plant #1 and
Plant #2 as the nominal model for simulating the closed loop response for Plant
#1 and Plant #2, respectively.
As we can see, both controllers perform well for their respective plants.
Because of the infinite output horizon, i.e., P = inf, nominal stability is
guaranteed for both systems. In some sense, the performance observed in
Figure 2-1 and Figure 2-2 is the best which can be expected, when the spark
advance constraint is invoked and there is no model/plant mismatch.
Obviously, if we want to control Plant #1 and Plant #2 with the same controller
the nominal performance for each plant will deteriorate.
Figure 2-2 Responses to a Unit Torque Disturbance for Plant #2
(no model/plant mismatch)
123456789100
123456789100
0
-2
-4
2
Outputs
Time
Manipulated Variables
Time
5
0
-5
y2
y1
u1
u2
2 MPC Based on Step Response Models
2-26
Figure 2-3 Responses to a Unit Torque Disturbance for Plant #1
(nominal model = Plant #2)
Simulation #2. Model/plant mismatch. All parameters are kept the same as in
Simulation #1. Shown in Figure 2-3 is the response to a unit torque disturbance
for Plant #1 using Plant #2 as the nominal model. Figure 2-4 depicts the
response to a unit torque disturbance for Plant #2 using Plant #1 as the
nominal model. As one can see, both systems are unstable. Therefore, the
controllers must be detuned to improve robustness if one wants to control both
plants with the same controller.
123456789100
123456789100
10
0
-10
20
Outputs
Time
Manipulated Variables
Time
20
0
-20
y2
y1
u1
u2
Application: Idle Speed Control
2-27
Figure 2-4 Responses to a Unit Torque Disturbance for Plant #2 (nominal
model = Plant #1)
Simulation #3. The input weight is increased to [10 20] to improve robustness.
All other parameters are kept the same as in Simulation #1. Plant #1 is used
as the nominal model. The simulation results depicted in Figure 2-5 and Figure
2-6 seem to indicate that with an input weight of [10 20] the controller
stabilizes both plants. However, we must point out that the design procedure
guarantees global asymptotic stability only for the nominal system, i.e., Plant
# 1. Because of the input constraints, the system is nonlinear. The observed
stability for Plant # 2 in Figure 2-6 should not be mistaken as an indication of
global asymptotic stability.
123456789100
123456789100
2
0
-2
4
Outputs
Time
Manipulated Variables
Time
5
0
-5
y2
y1
u1
u2
2 MPC Based on Step Response Models
2-28
Figure 2-5 Responses to a Unit Torque Disturbance for Plant #1 (nominal
model = Plant #1)
Figure 2-6 Responses to a Unit Torque Disturbance for Plant #2 (nominal
model = Plant #1)
123456789100
123456789100
0
10
Outputs
Time
Manipulated Variables
Time
5
-5
-10
y2
y1
u1
u2
0
123456789100
123456789100
0
-2
-4
2
Outputs
Time
Manipulated Variables
Time
4
2
0
y2
y1
u1
u2
Application: Idle Speed Control
2-29
As expected, the nominal performance for both Plant #1 and Plant #2 has
deteriorated when compared to the simulations shown in Figure 2-1 and Figure
2-2. A similar effect would be observed if we had detuned the controller which
uses Plant #2 as the nominal model.
Simulation #4. The parameter values are the same as in Simulation #3. Instead
of using cmpc,weusempccon and mpcsim for simulating the closed loop
responses. Figure 2-2 compares the responses for Plant #1 using mpccon and
mpcsim,andcmpc. As we can see, for this example and these tuning
parameters, the improvement obtained through the on-line optimization in
cmpc is small. However, the difference could be large, especially for
ill-conditioned systems and other tuning parameters. For example, by reducing
the output horizon to P = 80 while keeping the other parameters the same, the
responses for Plant # 1 found with mpccon and mpcsim are significantly slower
than those obtained with cmpc (Figure 2-8).
Figure 2-7 Comparison of Responses From cmpc, and mpccon and mpcsim for
Plant #1 P = inf
2468100
-4
-5
-6
-3
Solid: cmpc
Time
0
-1
-2
1
Dashed: mpccon and mpcsim
Ou
t
p
u
t
y1
y2
2 MPC Based on Step Response Models
2-30
Figure 2-8 Comparison of Responses From cmpc, and mpccon and mpcsim for
Plant #1 (P = 80)
2468100
-4
-5
-6
-3
Solid: cmpc
Time
0
-1
-2
1
Dashed: mpccon and mpcsim
Ou
tp
u
t
y1
y2
-7
Application: Control of a Fluid Catalytic Cracking Unit
2-31
Application: Control of a Fluid Catalytic Cracking Unit
Process Description
Fluid Catalytic Cracking Units (FCCUs) are widely used in the petroleum
refining industry to convert high boiling oil cuts (of low economic value) to
lighter more valuable hydrocarbons including gasoline. Cracking refers to the
catalyst enhanced thermal breakdown of high molecular weight hydrocarbons
into lower molecular weight materials. A schematic of the FCCU studied
4
is
given in Figure 2-9. Fresh feed is contacted with hot catalyst at the base of the
riser and travels rapidly up the riser where the cracking reactions occur. The
desirable products of reaction are gaseous (lighter) hydrocarbons which are
passed to a fractionator and subsequently to separation units for recovery and
purification. The undesirable byproduct of cracking is coke which is deposited
on the catalyst particles, reducing their activity. Catalyst coated with coke is
transported to the regenerator section where the coke is burned off thereby
restoring catalytic activity and raising catalyst temperature. The regenerated
catalyst is then transported to the riser base where it is contacted with more
fresh feed. Regenerated catalyst at the elevated temperature provides the heat
required to vaporize the fresh feed as well as the energy required for the
endothermic cracking reaction.
4. A detailed problem description and the model used for this study can be found in the paper by
McFarlane et al., “Dynamic Simulator for a Model IV Fluid Catalytic Cracking Unit,” Comp. &
Chem. Eng., 17(3), 1993, pages 275–300
2 MPC Based on Step Response Models
2-32
Figure 2-9 Fluid Catalytic Cracking Unit Schematic
Product composition, and therefore the economic viability of the process, is
determined by the cracking temperature. The bulk of the combustion air in the
regenerator section is provided by the main air compressor which is operated
at full capacity. Additional combustion air is provided by the lift air
compressor, the throughput of which is adjustable by altering compressor
speed. By maintaining excess flue gas oxygen concentration, it is possible to
ensure essentially complete coke removal from the catalyst.
Application: Control of a Fluid Catalytic Cracking Unit
2-33
Control Problem Formulation
The open loop system is modeled as follows:
where
G is the plant model and G
d
is the disturbance model. The variables are:
? Controlled variables
- Cracking temperature — (T
r
)
- Flue gas oxygen concentration — ( )
? Associated variable
- LiftAirCompressorSurgeIndicator—(D F
la
)
? Manipulated variables
- Liftaircompressorspeed—(V
lift
)
- Flue gas valve opening — (V
fg
)
? Modeled disturbances
- Variations in ambient temperature affect compressor throughput — (d
1
)
- Fluctuations in heavy oil feed composition to the FCCU — (d
2
)
- Pressure upset in down stream units propagating back to the FCCU — (d
3
)
In addition to the controlled variables there are many process variables that
need not be maintained at specific setpoints but which need to be within
bounds for safety or economic reasons. These variables are called associated
variables. For example, compressors must not surge during process upsets i.e.,
the suction flow rate must be greater than some minimum flow rate (surge flow
rate).
The control objective is to maintain the controlled variables (cracking
temperature and flue gas oxygen concentration) at pre-determined setpoints in
the presence of typical process disturbances while maintaining safe plant
operation.
yGuG
d
d+=
uV
fg
V
lift
[]
T
= yC
O
2 sg,
T
r
D F
la
[]
T
= dd
1
d
2
d
3
[]
T
=
C
O
2 sg,
2 MPC Based on Step Response Models
2-34
Simulations
State-space realizations of the plant and disturbance models are available in
the MATLAB file fcc_dat.mat in the directory mpcdemos. A MATLAB script
detailing the simulations is also included (fcc_demo.m). The following table
gives the parameters used for controller design and examination of the closed
loop response:
Table 2-1 FCCU Simulation Parameters
Step Response Model
Figure 2-10A shows the plant open loop step response to a unit step in V
fg
.
Although the plant is stable the settling time is large (1 day). The time scale of
interest for control purposes is on the order of one hour — which corresponds
to the initial plant response, Figure 2-10B. For time scales of one hour, the
process can be approximated by an integrating system. In deriving the step
response model, the plant is therefore assumed to be an integrating process.
Simulation Time(s) tend = 2500
# Step Response Coefficients 60
Process Sampling Time delt2 = 100
Output Weights ywt = [3 3 0]
Input Weights uwt = [0 2]
Prediction Horizon (steps) P = 12
# manipulated variable moves M = 3
input constraints u
i
? [-1,1],i =1,2
output constraints y
i
? [ -1,1] , i =1,2
y
3
?-1(hard constraint)
Application: Control of a Fluid Catalytic Cracking Unit
2-35
Figure 2-10 Open Loop Step Response to u = [1 0]
Figure 2-11 Unconstrained Closed Loop Response to d = [0 -0.5 -0.75]
0 2 4 6
x 10
4
?600
?500
?400
?300
?200
?100
0
time (s)
cO2sg
A
0 2000 4000 6000 8000 10000
?300
?250
?200
?150
?100
?50
0
time (s)
cO2sg
B
0 500 1000 1500 2000 2500
?2
?1
0
1
2
time (s)
output variables
cO2sg
Tr
surge indicator
Disturbance Rejection, UNCONSTRAINED DESIGN
0 500 1000 1500 2000 2500
?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
time (s)
manipulated variables
V
f
g
V
l
ift
2 MPC Based on Step Response Models
2-36
Associated Variables
As mentioned previously, associated variables need not be at any setpoint as
long as they are within acceptable bounds. Deviations of associated variables
from the nominal value do not appear in overall objective function to be
minimized and the output weight corresponding to the associated variable is
set to zero in Table 2-1.
Unconstrained Control Law
Figure 2-11 shows the closed loop response to a disturbance
d = [0 -0.5 -0.75] at t = 0. The controller gain matrix is derived using
mpccon and the closed loop response is examined using mpcsim.Notethe
following:
? At the first time step (t = 100s) the controlled variables are outside their
allowed limits. The onset of the disturbance at t = 0 is unknown to the
controller at t = 0 since there is no disturbance feedforward loop. Thus from
t = 0 to t = 100s there is no control action and the process response is the open
loop response with no control action. Only after t = 100s is corrective action
implemented.
? At t = 200s (2
nd
time step) riser temperature is outside the allowed limits.
? The lift air compressor surges during the interval 200s ￡ t ￡ 800s which is
unacceptable. Compressor surging will result in undesirable vibrations in
the compressor leading to rapid wear and tear.
Constrained Control Law
It is clear that the unconstrained control law generated using mpcsim is
physically unacceptable since hard output constraints are violated. Figure 2-12
shows the closed loop response of the nominal plant to the same disturbance
taking process constraints explicitly into account. The closed loop response is
determined using the command cmpc.
Application: Control of a Fluid Catalytic Cracking Unit
2-37
Figure 2-12 Constrained Closed Loop Response to d = [0 -0.5 -0.75]
The output limit vector is:
0 500 1000 1500 2000 2500
?2
?1
0
1
2
time (s)
output variables
Disturbance Rejection, CONSTRAINED DESIGN
cO2sg
Tr
surge indicator
0 500 1000 1500 2000 2500
?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
time (s)
manipulated variables
V
f
g
V
l
ift
ylim
1–1–1–11Inf
1–1–1–11Inf
Inf–Inf–Inf– Inf Inf Inf
=
2 MPC Based on Step Response Models
2-38
Note the following:
? At the first time step (t = 100s) the controlled variables are outside their
allowed limits. In fact the outputs are identical to the outputs for the
unconstrained case at t = 100s. This should be expected as there is no control
action from t = 0 to t = 100s for both constrained and unconstrained designs.
? At t = 200s (2
nd
time step) riser temperature (y
2
) is still outside the allowed
limits. This is because the constrained QP solved at t = 100s assumes that
disturbances are constant for t > 100s which is not the case for this process.
Thus while the manipulated variable move made at t = 100s ensures that the
predicted y
2
= 1 at t = 200s, the actual output at t = 200s exceeds one.
? The lift air compressor does not surge during the disturbance transient,
Figure 2-13.
The constrained control law therefore ensures safe operation while rejecting
process disturbances. If no constraints are violated, mpccon and mpcsim,and
cmpc will give identical closed loop responses. Note that the disturbance
d = [0 -0.5 -0.75] was specifically chosen to illustrate the possibility of
constraint violations during disturbance transient periods.
Figure 2-13 Comparison of Constrained and Unconstrained Response of D F
la
to d = [0 -0.5 -0.75]
0 500 1000 1500 2000 2500
?1.2
?1
?0.8
?0.6
?0.4
?0.2
0
0.2
Unconstrained Response
Constrained Response
time (s)
associated variable
Associated Variable Response to Constrained & Unconstrained Control Laws
3
MPC Based on
State-Space Models
3 MPC Based on State-Space Models
3-2
State-Space Models
Consider the process shown in the above block diagram. The general
discrete-time linear time invariant (LTI) state-space representation used in
the MPC Toolbox is:
where x is a vector of n state variables, u represents the n
u
manipulated
variables, d represents n
d
measured but freely-varying inputs (i.e., measured
disturbances), w represents n
w
unmeasured disturbances, y is a vector of n
y
plant outputs, z is measurement noise, and F, G
u
, etc., are constant matrices of
appropriate size. The variable (k) represents the plant output before the
addition of measurement noise. Define:
S
y y
Unmeasured
Disturbances
Manipulated
Variables
Measured
Disturbances
Measurement
Noise
w
u
d
++
z
Plant
Measured
Outputs
xk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() zk()+=
Cx k() D
u
uk() D
d
dk() D
w
wk() zk()+++ +=
y
G
G
u
G
d
G
w
=
D
D
u
D
d
D
w
=
State-Space Models
3-3
Inmanyapplications,alloutputsaremeasured.Insomecases,however,one
has n
ym
measured and n
yu
unmeasured outputs in y,wheren
ym
+ n
yu
= n
y
.If
so, the MPC Toolbox assumes that the y vector and the C and D matrices are
arranged such that the measured outputs come first, followed by the
unmeasured outputs.
Mod Format
The MPC Toolbox works with state-space models in a special format, called the
mod format. The mod format is a single matrix that contains the state-space
F, G , C,andD matrices plus some additional information (see mod format in
the "Command Reference" chapter for details). The MPC Toolbox includes a
number of commands that make it easy to generate models in the mod format.
The following sections illustrate the use of these commands.
SISO Continuous-Time Transfer Function to Mod
Format
The MPC Toolbox uses a format called the tf format. Let the continuous-time
transfer function be
where T
d
is the time delay. The tf format is a matrix consisting of three rows:
The tf matrix will always have at least two columns, since that is the
minimum width of the third row.
row 1: The n coefficients of the numerator polynomial, b
0
to b
n
.
row 2: The n coefficients of the denominator polynomial, a
0
to a
n
.
row 3: column 1: The sampling period. This must be zero for a
continuous system. (It must be positive for discrete transfer
functions — see next section).
column 2: The time delay in time units. It must satisfy Td ? 0.
Gs()
b
0
s
n
b
1
s
n 1–
… b
n
+++
a
0
s
n
a
1
s
n 1–
… a
n
+++
----------------------------------------------------------------e
T–
d
s
=
3 MPC Based on State-Space Models
3-4
You can either define a model in the tf format directly or use the command
poly2tfd. The general form of this command is
g = poly2tfd(num,den,delt,delay)
For example, consider a SISO system modeled by the transfer function
To create the tf format directly you could use the command
G = [0 -13.6 1; 54.3 113.5 1; 0 5.3 0];
which defines a matrix consisting of three rows and three columns. Note that
all rows must have the same number of columns so you must be careful to
insert zeros where appropriate. The poly2tfd command is more convenient
since it does that for you automatically:
G = poly2tfd([-13.6 1],[54.3 113.5 1],0,5.3);
Either command would define a variable G in your workspace, containing the
matrix
To convert this to the mod format, use the command tfd2mod, which has the
form
model = tfd2mod(delt,ny,g1,g2,g3,...,gN)
G=
0 -13.6000 1.0000
54.3000 113.5000 1.0000
0 5.3000 0
Gs()
13.6s– 1+
54.3s
2
113.5s 1++
--------------------------------------------------- e
5.3s–
=
State-Space Models
3-5
where:
Suppose you want to convert the SISO model defined above to the mod format
with a sampling period of 2.1 time units. The appropriate command would be
mod = tfd2mod(2.1,1,G);
This would define a variable called mod in your workspace that would contain
the discrete-time state-space description of your system.
delt The sampling period. tfd2mod will convert your
continuous time transfer function(s) g1, ..., gN to
discrete-time using this sampling period.
ny is the number of output variables in the plant you are
modeling.
g1,g2,...,gN A sequence of N transfer functions in the tf format, where
N ? 1. tfd2mod assumes that these are the individual
elements of a transfer-function matrix:
Thus N must be an integer multiple (n
u
) of the number of
outputs, n
y
. You supply the transfer functions in
column-wise order. In other words, you first give the n
y
transfer functions for input 1 ( ), then the n
y
transfer functions for input 2 ( ), etc.
g
11,
g
12,
… g
1 n
u
,
g
21,
g
22,
… g
2 n
u
,
g
n
y
1,
g
n
y
2,
… g
n
y
n
u
,
… … …
.
.
.
g
11,
to g
n
y
1,
g
12,
to g
n
y
2,
3 MPC Based on State-Space Models
3-6
SISO Discrete-Time Transfer Function to Mod Format
Suppose you have a transfer function in discrete-time format (in terms of the
forward-shift operator, z):
where d is an integer (? 0) and represents the sampling periods of pure delay.
The corresponding tf formatisthesameasforthecontinuous-timecaseexcept
for the definition of the third row:
As in the previous section, you can use poly2tfd followed by tfd2mod to get
such a transfer function in mod format. For example, the discrete-time
representation of the SISO system considered in the previous section is
If you had this to begin with, you could convert it to the mod format as follows:
G = poly2tfd([-0.1048 0.1215 0.0033],[1 -0.9882 0.0082], 2.1,3);
mod = tfd2mod(2.1,1,G);
Note that both the poly2tfd and tfd2mod commands specify the same
sampling period (delt=2.1). This would be the usual case, but you have the
option of converting a discrete-time model in the tf format to a different
sampling period in the mod format.
column 1 is the sampling period for which G(z) was created. It must
be positive (in contrast to the continuous-time case
described above).
column 2 is the periods of pure delay, d, which must be an integer ? 0.
Contrast this to the continuous case, where the delay is
given in time units.
Gq()
b
0
b
1
z
1–
… b
n
z
n–
+++
a
0
a
1
q
1–
… a
n
z
n–
+++
-------------------------------------------------------------- z
d–
=
Gz()
0.1048– 0.1215z
1–
0.0033z
2–
++
1 0.9882z
1–
0.0082z
2–
+–
----------------------------------------------------------------------------------------- z
3–
=
State-Space Models
3-7
MIMO Transfer Function Description to Mod Format
Suppose you have a transfer-function matrix description of your system in the
form
where g
i,j
is the transfer function of the i
th
output with respect to the j
th
input.
If all n
y
outputs are measured and all n
u
inputs are manipulated variables, the
default mode of tfd2mod will give you the correct mod format. For example,
consider the 2-output, 3-input system:
The following sequence of commands would convert this to the equivalent mod
format with a sampling period of T =4:
g11 = poly2tfd(12.8,[16.7 1],0,1);
g21 = poly2tfd(6.6,[10.9 1],0,7);
g12 = poly2tfd(-18.9,[21.0 1],0,3);
g22 = poly2tfd(-19.4,[14.4 1],0,3);
g13 = poly2tfd(3.8,[14.9 1],0,8);
g23 = poly2tfd(4.9,[13.2 1],0,3);
pmod = tfd2mod(4,2,g11,g21,g12,g22,g13,g23);
g
11,
g
12,
… g
1 n
u
,
g
21,
g
22,
… g
2 n
u
,
g
n
y
1,
g
n
y
2,
… g
n
y
n
u
,
… … …
.
.
.
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
3.8e
8s–
14.9s 1+
------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
4.9e
3s–
13.2s 1+
------------------------
u
1
s()
u
2
s()
u
3
s()
=
3 MPC Based on State-Space Models
3-8
Suppose, however, that the third input were actually an unmeasured
disturbance, i.e., the system were
In this case you would need to override the default mode of tfd2mod by
specifying the number of inputs in each of the three categories described at the
beginning of this section, i.e., manipulated variables, measured disturbances,
and unmeasured disturbances. This and other information about the system is
contained in the first 7 columns of row 1 of the mod format, as follows:
Forexample,ifyouhaddefinedpmod using the default mode of tfd2mod as
shown above, the contents of row 1, columns 1 to 7 of pmod would be:
4 13 3 0 0 2 0
You could override this to set n
u
=2andn
w
=1asfollows:
pmod(1,3)=2;
pmod(1,5)=1;
column 1 T, the sampling period.
2 n, the number of states.
3 n
u
, the number of manipulated variable inputs.
4 n
d
, the number of measured disturbances.
5 n
w
, the number of unmeasured disturbances.
6 n
ym
, the number of measured outputs.
7 n
yu
, the number of unmeasured outputs.
y
1
s
y
2
s
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s
u
2
s
=
3.8e
8s–
14.9s 1+
------------------------
4.9e
3s–
13.2s 1+
------------------------
ws()+
State-Space Models
3-9
Note that in the original transfer function matrix description, the first n
u
columns must be for the manipulated variables, the next n
d
for the measured
disturbances (if any), and the last n
w
for the unmeasured disturbances (if any).
Similarly, the first n
ym
outputs must be measured and the last n
yu
(? 0)
unmeasured.
Continuous or Discrete State-Space to Mod Format
If you have a continuous-time state-space model, you may convert it to mod
format by first using the function c2dmp (continuous to discrete-time
state-space), followed by ss2mod(discrete-time state-space to mod format). Of
course, if you are starting with a discrete-time state-space model you can skip
the first step.
For example, suppose a, b, c,andd are matrices describing a continuous-time
system. To convert to the mod format using a sampling period of T =1.5,you
could use the following commands:
[phi,gam] = c2dmp(a,b,1.5);
mod = ss2mod(phi,gam,c,d,1.5);
If your system is complicated, i.e., it contains disturbance inputs and/or
unmeasured outputs, you will need to override the default mode of ss2mod.See
the “Command Reference” section for more details.
Identification Toolbox (“Theta”) Format to Mod
Format
The System Identification Toolbox identifies discrete-time transfer-function
models from input/output data. The result is a model in a special form called
the theta format (see the System Identification Toolbox User’s Guide for
details).Ingeneral,eachtheta model describes the response of a single output
tooneormoreinputs(MISOmodel).
The MPC Toolbox function, th2mod, converts one or more such models to the
mod format. Suppose, for example, that
th1 is the theta model describing the response of output y
1
to inputs u
1
and u
2
.
th2 is the theta model describing the response of output y
2
to inputs u
1
and u
2
.
3 MPC Based on State-Space Models
3-10
Then the following command would provide the equivalent mod format with
n
y
=2andn
u
=2:
mod = th2mod(th1,th2);
Combination of Models in Mod Format
The functions addmod, addmd, addumd, appmod, paramod,andsermod allow you
to combine simple models in the mod format to generate more complex plant
structures. For example, addmd includes the effect of one or more measured
disturbances in an existing model, as shown in the following schematic:
pmod gives the effect of one or more manipulated variables, u(k), and optional
unmeasured disturbance(s), w(k), on the output(s), y(k). dmod gives the effect of
the measured disturbance(s), d(k), on the same outputs. Once you have defined
pmod and dmod (e.g., starting from transfer functions as illustrated above), you
can use the command addmd to generate the composite, model:
model = addmd(pmod,dmod);
Please see Chapter 4, “Command Reference” for more details on the various
model-building commands.
Converting Mod Format to Other Model Formats
The function mod2ss converts a model in the mod format to the standard
discrete-time state-space format:
[phi,gam,c,d,minfo] = mod2ss(mod);
Here, phi, gam, c,andd are the coefficient matrices of
dmod
w
d
Model
S
y
+
+
u
y
d
pmod
w
u
xk 1+()Fxk() Guk()+=
State-Space Models
3-11
The vectorminfocontains the first 7 columns of the first row inmod.Thesection
“MIMO Transfer Function Description to Mod Format” gives the significance
of this information.
Once you have phi, gam, c,andd,youcanused2cmp, ss2tf2, and other
functions to convert from discrete state-space to other model forms.
The function mod2step uses a model in the mod format to generate a step-
response model in the step format as required by the functions mpccon,mpcsim,
etc., discussed in Chapter 2, “MPC Based on Step Response Models”. See the
Chapter 4, “Command Reference” for details on the use of mod2step.
yk() Cx k() Du k()+=
3 MPC Based on State-Space Models
3-12
Unconstrained MPC Using State-Space Models
Once you have described your system by generating state-space models in the
mod format you can use the commands:
In addition, you can analyze certain properties of the closed-loop system using
the commands:
Note: smpcgain, smpcpole and mod2frsp also work with open-loop models in
the mod format.
smpccon to calculate the unconstrained controller gain matrix.
smpcest to design a state estimator (optional).
smpcsim tosimulatetheresponseoftheclosed-loopsystemtooneor
more specified inputs.
plotall
(or ploteach) to plot the response(s).
smpccl to generate a model of the closed-loop system (plant plus
controller).
smpcgain to calculate the closed-loop gain matrix.
smpcpole to calculate the closed-loop poles.
mod2frsp (and plotfrsp)tocalculateandplottheclosed-loop
frequency response.
svdfrsp to calculate the singular values of the frequency response.
Unconstrained MPC Using State-Space Models
3-13
Example: (see mpctutss.m)
The following example (mpctutss.m) illustrates the basic procedures. The
example process has 2 measured outputs, 2 manipulated variables, and an
unmeasured disturbance:
We first define the model in the mod format. The following commands use a
sampling period of T = 2 time units (chosen arbitrarily):
delt = 2;
ny = 2;
g11 = poly2tfd(12.8,[16.7 1],0,1);
g21 = poly2tfd(6.6,[10.9 1],0,7);
g12 = poly2tfd(-18.9,[21.0 1],0,3);
g22 = poly2tfd(-19.4,[14.4 1],0,3);
umod = tfd2mod(delt,ny,g11,g21,g12,g22);
% Defines the effect of u inputs
g13 = poly2tfd(3.8,[14.9 1],0,8);
g23 = poly2tfd(4.9,[13.2 1],0,3);
dmod = tfd2mod(delt,ny,g13,g23);
% Defines the effect of w input
pmod = addumd(umod,dmod); % Combines the two models.
We now design an unconstrained MPC controller. The design parameters are
essentially the same as for the functions based on step response models (see
Chapter 2). In this case, start by choosing design parameters such that we get
the perfect controller:
imod = pmod;% assume perfect modeling
ywt = [ ]; % default (unity) weights on both outputs
uwt = [ ]; % default (zero) weights on both inputs
P = 5; % prediction horizon
M = P; % control horizon
Ks = smpccon(imod,ywt,uwt,M,P);
y
1
s
y
2
s
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s
u
2
s
=
3.8e
8s–
14.9s 1+
------------------------
4.9e
3s–
13.2s 1+
------------------------
ws()+
3 MPC Based on State-Space Models
3-14
We check the design by running a simulation for a step increase in the setpoint
of output y
1
:
tend=30; % time period for simulation.
r = [1 0]; % setpoints for the two outputs.
[y,u] = smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,delt)
Notethatthereisnomodelerrorsinceweusedthesamemodeltorepresent
the plant (pmod) as that used to design the controller (imod). The results are:
Note that we get perfect tracking of the specified setpoint change (y
1
=1,
y
2
= 0), but the manipulated variables are ringing. You could have anticipated
this by calculating the poles of the controller:
[clmod,cmod] = smpccl(pmod,imod,Ks);
smpcpole(cmod)
The result shows that one of the poles is at -0.9487 and another is at -0.9223.
In general, such negative-real poles cause ringing.
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?1.5
?1
?0.5
0
0.5
1
1.5
Manipulated Variables
Time
Unconstrained MPC Using State-Space Models
3-15
One way to minimize ringing is to make the prediction horizon significantly
larger than the control horizon:
P = 10;
M = 3;
Ks = smpccon(imod,ywt,uwt,M,P);
[y,u] = smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,delt)
This results in the following improved responses:
Another (often more effective) way is to use blocking.Inthecaseofblocking,
each element of the vector M indicates the number of steps over which D u =0
during the optimization. For example, M = [2 3] indicates that u(k +1)=u(k)
or D u(k +1)=0andu(k +4)=u(k +3)=u(k +2)(orD u(k +3)=D u(k +4)=0):
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.4
?0.2
0
0.2
0.4
0.6
0.8
1
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-16
M = [2 3 4]; % Defines 3 blocks of control moves
Ks = smpccon(imod,ywt,uwt,M,P);
[y,u] = smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,delt)
pause
This completely eliminates the ringing, as shown in the following responses, at
the expense of a more sluggish servo response and a larger disturbance in y
2
.
A third approach is to increase the weights on the manipulated variables:
uwt = [1 1]; % increase input weighting
P = 5; % original prediction horizon
M = P; % original control horizon
Ks = smpccon(imod,ywt,uwt,M,P);
[y,u] = smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,delt)
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
Manipulated Variables
Time
Unconstrained MPC Using State-Space Models
3-17
for which the response is:
In general, you must choose the horizons and weights by trial-and-error, using
simulations to judge their effectiveness.
The servo-response of the last controller looks good. Let’s see how it responds
to a unit-step in the unmeasured disturbance, w(k):
ulim = [ ]; % default (no) constraints on u variables.
Kest = [ ]; % default (DMC) state estimator.
r = [0 0]; % Both output setpoints at zero.
z = [ ]; % default (zero) measurement noise.
v = [ ]; % default (zero) measured disturbances.
w = [1]; % unit-step in unmeasured disturbance.
[y,u] = smpcsim(pmod,imod,Ks,tend,r,ulim,Kest,z,v,w);
plotall(y,u,delt)
pause
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-18
Note that both setpoints are zero. The resulting response is:
The regulatory response has a rather long transient. Let’s see if we can improve
it by using a state estimator other than the default (DMC) estimator:
[Kest,newmod] = smpcest(imod,[15 15],[3 3]);
Ks = smpccon(newmod,ywt,uwt,M,P);
[y,u] = smpcsim(pmod,newmod,Ks,tend,r,ulim,Kest,z,v,w);
plotall(y,u,delt)
0 5 10 15 20 25 30
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
Manipulated Variables
Time
Unconstrained MPC Using State-Space Models
3-19
See the detailed description of the smpcest function for a discussion of the
estimator design parameters. The results show that the controller now
compensates for the disturbance much more rapidly:
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-20
State-Space MPC with Constraints
The function scmpc handles problems with inequality constraints on the
manipulated variables and/or outputs. The recommended procedure is to first
use the tools described in the section section Unconstrained MPC Using
State-Space Models to find values of the prediction horizon, P, control horizon,
M, input and output weights, and a state-estimation strategy that work well for
the unconstrained version of your problem. Then define the constraints and
solve the problem using scmpc. The following example illustrates the use of
scmpc.
Example: (see mpctutss.m)
We use the same example process as in the previous section, but use a sampling
period of 1 and omit the unmeasured disturbance input:
T = 1;
g11 = poly2tfd(12.8,[16.7 1],0,1);
g21 = poly2tfd(6.6,[10.9 1],0,7);
g12 = poly2tfd(-18.9,[21.0 1],0,3);
g22 = poly2tfd(-19.4,[14.4 1],0,3);
imod = tfd2mod(2,T,g11,g21,g12,g22);
The following statements specify parameters required in both the constrained
and unconstrained cases, and calculate the gain for the unconstrained
controller.
nhor = 10; % Prediction horizon.
ywt = [ ]; % Unity weighting on output tracking errors
% (default).
uwt = [ ]; % Zero weighting on man. variable moves
% (default).
blks = [2 3 5]; % Allows 3 moves of manipulated variables.
K = [ ]; % DMC-type state estimation (default).
Ks = smpccon(imod,ywt,uwt,blks,nhor);
State-Space MPC with Constraints
3-21
Let’s first verify that the constrained and unconstrained solutions will be the
same when the constraints are loose i.e. inactive. The following statements
define upper and lower bounds on u(k)at-￥ and ￥ , respectively, and bounds on
D u(k)at10(bothu
1
and u
2
).
1
Also, bounds on y(k) are set at the default values
of ±￥ .
ulim = [-inf -inf inf inf 10 10];
ylim = [ ]; % Default -- no limits on outputs.
For the simulation we will make a step change of 0.8 in the setpoint for y
1
.We
will also assume a perfect model, i.e., use the same model for the plant as was
used to design the controller.
setpts = [0.8 0]; % Define the step in the setpoint.
plant = imod;
tend = 20; % Duration of the simulation
[y1,u1] = smpcsim(plant,imod,Ks,tend,setpts,ulim,K);
[y,u] = scmpc(plant,imod,ywt,uwt,blks,nhor,tend,...
setpts,ulim,ylim,K);
plotall([y y1],[u u1],T)
The above plotall command plots the results from smpcsim and scmpc on the
same graph. Since the constraints were loose, there should be no difference. In
the following plots, you can only distinguish two curves, i.e., the two
simulations give the same values of y and u,asexpected.
1. Finite bounds on D u are required by scmpc. Here they are chosen large enough so that they
have no effect.
3 MPC Based on State-Space Models
3-22
Now let’s add some constraints to the problem. Suppose we want the maximum
value of y
2
to be -0.1. In the previous case it goes slightly above zero (dashed
line in the Outputs plot). The following statements define a hard upper limit of
y
2
= -0.1, starting at the 4th step in the prediction horizon. This accounts for the
minimum delay of 3 sampling periods before y
2
can be affected by either u
1
or
u
2
, i.e., it is important to leave y
2
unconstrained for the first 3 steps in the
prediction horizon. In this case, since the initial condition is y
2
=0,itis
impossible to make y
2
￡ -0.1 prior to t = 4. If you were to attempt to do so, you
would get an error message stating that the problem is infeasible. Note also
that the upper bound on y
2
supersedes the setpoint, which is still specified as
zero. The controller thus maximizes the value of y
2
at steady state.
ylim = [-inf -inf inf inf
-inf -inf inf inf
-inf -inf inf inf
-inf -inf inf -0.1];
[y,u] = scmpc(plant,imod,ywt,uwt,blks,nhor,tend,...
setpts,ulim,ylim,K);
0 2 4 6 8 10 12 14 16 18 20
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 2 4 6 8 10 12 14 16 18 20
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Manipulated Variables
Time
State-Space MPC with Constraints
3-23
The following plot shows that the constraints are satisfied for t ? 4, as expected.
If you use a long prediction horizon with constraints, the calculations can be
time-consuming. You can minimize this by turning off constraints that are far
out in the prediction horizon. The following example defines bounds on y
1
and
y
2
, then turns them off beyond the 4th point in the prediction horizon. The
calculations are much faster than would be the case if only the first 4 rows of
ylim had been used (try it). Also, since there is neither model error nor
unmeasured disturbances, the solution satisfies all constraints for t 4 in any
case. In general, output constraints must be chosen carefully to avoid
infeasibilities and maximize the speed of the calculations.
ylim = [-inf -inf inf inf
-inf -inf 0.8 inf
-inf -inf inf inf
-inf 0.10 inf inf]
-inf -inf inf inf
% Turns off remaining bounds.
0 2 4 6 8 10 12 14 16 18 20
?0.4
?0.2
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 2 4 6 8 10 12 14 16 18 20
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-24
[y,u] = scmpc(plant,imod,ywt,uwt,blks,nhor,tend,...
setpts,ulim,ylim,K);
plotall(y,u,T)
As a final example we impose bounds on the manipulated variables:
ulim = [-0.5 -0.5 0.5 0 0.3 0.3
-inf -inf inf inf 0.3 0.3];
[y,u] = scmpc(plant,imod,ywt,uwt,blks,nhor,tend,...
setpts,ulim,ylim,K);
plotall(y,u,T)
Again, to save computer time the constraints apply only for the first block in
the prediction horizon, i.e., the constraints are turned off for periods 2 through
P = 10. The following plot shows that the upper bound of u
2
￡ 0andD u
1
￡ 0.3
are the most restrictive. The former prevents y
2
from coming back to the
minimum allowed value of 0.1.
0 2 4 6 8 10 12 14 16 18 20
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 2 4 6 8 10 12 14 16 18 20
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Manipulated Variables
Time
State-Space MPC with Constraints
3-25
0 2 4 6 8 10 12 14 16 18 20
0
0.2
0.4
0.6
0.8
1
Outputs
Time
0 2 4 6 8 10 12 14 16 18 20
?0.1
0
0.1
0.2
0.3
0.4
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-26
Application: Paper Machine Headbox Control
Ying et al. (1992)
2
studied the control of composition and liquid level in a paper
machine headbox, a schematic of which is shown in Figure 3-1. The process
model is given by a set of ordinary differential equations (ODEs) in bilinear
form. Using their nomenclature, the states are x
T
=[H
1
H
2
N
1
N
2
], where H
1
is the liquid level in the feed tank, H
2
is that in the headbox, N
1
is the
consistency (percentage of pulp fibers in suspension) in the feed tank, and N
2
is that in the headbox. All states except H
1
are measured, i.e., the measured
outputs are y
T
=[H
2
N
1
N
2
]. The primary control objective is to hold H
2
and N
2
(y
1
and y
3
)atspecifiedsetpoints.
Therearetwomanipulatedvariables:u
T
=[G
p
G
w
], where G
p
is the flowrate of
stock entering the feed tank, and G
w
is that of the recycled white water.There
is a single measured disturbance: v =[N
p
], the consistency of the stock entering
the feed tank, and a single unmeasured disturbance: d =[N
w
], the consistency
of the white water. All variables are normalized such that they are zero at the
nominal steady state, and variations about the steady-state are of the same
order of magnitude. The process is open-loop stable.
2. Ying, Y., M. Rao, and Y. Sun, “Bilinear Control Strategy for Paper-Making Process,” Chem. Eng.
Comm. 1992, 111, 13–28.
Application: Paper Machine Headbox Control
3-27
Figure 3-1 Schematic of Paper Machine Headbox Control Problem
MPC Design Based on Nominal Linear Model
The standard MPC design methods require a linear model of the plant. We
therefore linearize the bilinear model at the nominal steady-state condition
(x =0;u =0;v=0;d = 0). Since the model is simple, one can linearize it
analytically to obtain:
where x
m
,y
m
,andd
m
are the model states, outputs, and disturbances,
respectively. The desired closed-loop response time is of the order of 10
minutes, so we choose a sampling period of T
s
= 2 minutes. The file pm_lin.m
in the directory mpcdemos contains the code for all the computations in this
section of the manual. The following commands define the linear model and
plot the response of the outputs to a unit step in each manipulated variable:
Feed
Tank
Headbox
Wire
White Water
Stock
Wet Paper
H
1
N
1
H
2
N
2
N
w
N
p
G
p
G
w
x
·
m
Ax
m
B
0
uB
v
vB
d
d
m
+++=
3 MPC Based on State-Space Models
3-28
Figure 3-2 Responses of Paper Machine Outputs to Unit Step in u
1
% Matrices= of the linearized paper machine model
A = [-1.93 0 0 0; .394 -.426 0 0; 0 0 .63 0; .82 -.784
.413 -.426];
B = [1.274 1.274 0 0;0 0 0 0;1.34 -.65 .203 .406;0 0 0 0];
C = [0 1 0 0; 0 0 1 0; 0 0 0 1];
D = zeros(3,4);
% Discretize the linear model and save in mod form.
dt = 2;
0 10 20 30
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
u1 step response : y1
TIME
0 10 20 30
0
0.5
1
1.5
2
2.5
u1 step response : y2
TIME
0 10 20 30
0
0.5
1
1.5
2
2.5
u1 step response : y3
TIME
Application: Paper Machine Headbox Control
3-29
[PHI,GAM] = c2dmp(A,B,dt);
minfo = [dt,4,2,1,1,3,0];
imod = ss2mod(PHI,GAM,C,D,minfo);
plotstep(mod2step(imod,30))
Figure 3-3 Responses of Paper Machine Outputs to Unit Step in u
2
0 10 20 30
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
u2 step response : y1
TIME
0 10 20 30
?1.4
?1.2
?1
?0.8
?0.6
?0.4
?0.2
0
u2 step response : y2
TIME
0 10 20 30
?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
u2 step response : y3
TIME
3 MPC Based on State-Space Models
3-30
The step responses (Figure 3-2 and Figure 3-3) show that there are large
interactions in the open-loop system. Adjustments in u
1
and u
2
have strong
effects on both y
1
and y
3
.Also,theu
2
fi y
3
step exhibits an inverse response.
We begin by attempting a controller design for good response of both y
1
and y
3
:
% Define controller parameters
P = 10; % Prediction horizon
M = 3; % Control horizon
ywt = [1,0,1]; % Equal weighting of y(1) and y(3),
% no control of y(2)
uwt = 0.6*[1 1]; % Equal weighting of u(1) and u(2).
ulim = [-10*[1 1] 10*[1 1] 2*[1 1]];% Constraints on u
ylim = [ ]; % No constraints on y
Kest = [ ]; % Use default estimator
% Simulation using scmpc -- no model error
pmod=imod; % plant and internal model are identical
setpts = [1 0 0];
% servo response to step in y(1) setpoint
tend = 30; % duration of simulation
[y,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest);
plotall(y,u,dt)
The prediction horizon of 10 sampling periods (20 minutes) extends well past
the desired closed-loop response time. Preliminary trials suggested that longer
horizons increased the computational load but provided no advantage in
setpoint tracking. The use of M < P is not required in this case, but helps to
reduce the computational load and inhibit ringing of the manipulated
variables. Note the equal penalties on setpoint tracking errors for y
1
and y
3
(ywt variable), reflecting our desire to track both setpoints accurately. There is
no penalty on y
2
, since it does not have a setpoint. The listeduwtpenalties were
determined by running several trials. Figure 3-4 shows smooth control of y
1
with the desired 10-minute response time, but there is a noticeable disturbance
in y
3
.
Application: Paper Machine Headbox Control
3-31
Figure 3-4 Response of closed-loop system to unit step in y
1
setpoint for equal
output weighting. Output y
2
is uncontrolled, and the y
3
setpoint is zero
One could repeat the simulation for a step change in the y
3
setpoint as follows:
setpts = [0 0 1]; % servo response to step in y(3) setpoint
[y,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest);
plotall(y,u,dt)
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
1.5
Outputs
Time
0 5 10 15 20 25 30
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-32
Normally, control of consistency (y
3
) is more important than control of liquid
level, y
1
. We can achieve better control of y
3
if we allow larger tracking errors
in y
1
. For example, an alternative controller design uses unequal weights on
the controlled outputs:
ywt = [0.2,0,1]; % Unequal weighting of y(1) and y(3),
% no control of y(2)
setpts = [1 0 0];
% servo response to step in y(1) setpoint
[y,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest);
plotall(y,u,dt)
AsshowninFigure3-5,ay
1
setpoint change causes a much smaller
disturbance in y
3
than before (compare with Figure 3-4). The disadvantage is
thattheresponsetimeofy
1
has increased from about 8 to 25 minutes.
Similarly, a step change in the y
3
setpoint would cause a larger disturbance in
y
1
than in the original design. Overall, however, the controller with unequal
weighting gives better nominal performance and will be used as the basis for
subsequent designs.
Application: Paper Machine Headbox Control
3-33
Figure 3-5 Response of closed-loop system to Unit step in y
1
setpoint for
unequal output weighting. Output y
2
is uncontrolled, and the y
3
setpoint is
zero.
We now evaluate the response of the above controller to a unit step in the
measured disturbance, v (i.e., feedforward control). The commands required for
this are:
setpts = [0 0 0]; % output
setpoints z = [ ]; % measurement noise
v = 1; % measured disturbance
d = 0; % unmeasured disturbance [y,u,ym] =
scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest,z,v,d);
plotall(y,u,dt)
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
Outputs
Time
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-34
As shown in Figure 3-6, both controlled outputs are held near their setpoints,
with larger deviations in y
1
,asexpected.
Figure 3-6 Response of closed-loop system to unit step in measured
disturbance, v. unequal output weighting with y
1
and y
3
setpoints at zero.
Finally, we check the response to a step in the unmeasured disturbance. The
required commands are:
setpts = [0 0 0]; % output setpoints
v = 0; % measured disturbance
d = 1; % unmeasured disturbance
[y,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest,z,v,d);
plotall(y,u,dt)
0 5 10 15 20 25 30
?0.04
?0.02
0
0.02
0.04
0.06
0.08
0.1
Outputs
Time
0 5 10 15 20 25 30
?0.2
?0.15
?0.1
?0.05
0
0.05
0.1
0.15
Manipulated Variables
Time
Application: Paper Machine Headbox Control
3-35
As shown in Figure 3-7, the unmeasured disturbance causes significant
deviations in both controlled outputs. In fact, the higher-priority output, y
3
,
exhibits the larger tracking error of the two controlled variables.
Figure 3-7 Response of closed-loop system (with default state estimator) to
unit step in unmeasured disturbance, d. unequal output weighting with y
1
and y
3
setpoints at zero.
The closed-loop response to unmeasured disturbances can often be improved by
a change in the state estimator. In the previous trials, we were using the
default estimator, which assumes that disturbances are independent, random
steps at each output. In fact, the known unmeasured disturbance, d,hasno
effect on y
1
, and its effects on y
2
and y
3
are approximately first order with time
constants of 3 and 5 minutes, respectively. One way to exploit this knowledge
is to specify an expected covariance for d and a measurement noise covariance
for y,thenusetheKalmangainforthemodeleddisturbancecharacteristics.
0 5 10 15 20 25 30
?0.2
?0.1
0
0.1
0.2
0.3
0.4
0.5
Outputs
Time
0 5 10 15 20 25 30
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
Manipulated Variables
Time
3 MPC Based on State-Space Models
3-36
Consider the following sequence of commands, leading to the responses shown
in Figure 3-8:
Figure 3-8 Response of closed-loop system to unit step in unmeasured
disturbance, d. with Kalman Estimator, unequal output weighting with y
1
and
y
3
setpoints at zero.
% Estimator design
Q = 30;
R = 1*eye(3);
Kest = smpcest(imod,Q,R);
% Simulation using scmpc -- no model error
setpts = [0 0 0]; % servo response to step in y1 setpoint
d = 1; % unmeasured disturbance
[y,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest,z,v,d);
plotall(y,u,dt)
0 5 10 15 20 25 30
?0.2
?0.1
0
0.1
0.2
0.3
0.4
0.5
Outputs
Time
0 5 10 15 20 25 30
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
Manipulated Variables
Time
Application: Paper Machine Headbox Control
3-37
The resulting Kest matrix is:
Kest =
0 0 0
0.0000 -0.0000 0.0000
0.0000 0.7446 0.0732
-0.0000 0.2985 0.0851
0.0000 -0.0000 0.0000
-0.0000 0.7777 0.2934
0.0000 0.2934 0.1977
We have specified equal measurement noise for each output (R is diagonal with
a rank equal to the number of outputs, and equal elements on the diagonal).
This makes the Kalman estimator give equal weight to each output
measurement.
3
The dimension of Q must equal the number of elements in d
(unity in this case). A relatively large value of Q(i,i) signifies an important
disturbance mode. In practice, the elements of Q and R are tuning parameters,
and one adjusts the relative magnitudes to achieve the desired balance of fast
disturbance rejection (usually promoted by making Q relatively large) and
robustness.
For the chosen Q and R, and the disturbance model in imod, the elements of
column 1 of Kest (shown above) are essentially zero. Thus, the measurement of
y
1
provides no information regarding the effect of d on the process states.
Output y
2
, on the other hand, provides large corrections to the state estimates.
If it were not available, rejection of d would degrade.
4
Figure 3-8 shows that although the revised estimator reduced the disturbance
in y
3
, it is still significant (compare to Figure 3-7). A key limiting factor is the
use of a 2-minute sampling period. As shown in Figure 3-8, the controller does
not respond to the disturbance until it is first detected at t = 2 minutes. You can
verify that reducing the sampling period to 0.25 minutes (holding all other
parameters constant) greatly reduces the disturbance in y
3
.Suchachange
would also speed up the setpoint tracking in the nominal case. It may cause
robustness problems, however, so we defer further consideration of the
samplingperiodtotestswiththenonlinearplant(seenextsection).
This application is unusual in that the characteristics of the unmeasured
disturbance are known. When this is not the case, the output disturbance form
3. If a measurement were known to be inaccurate, its R(i,i) value should be relatively large.
4. You can see how serious the degradation would be by setting R(2,2) to a large value,
e.g.,10000.
3 MPC Based on State-Space Models
3-38
of the estimator simplifies the design procedure. It requires only a rough idea
of the characteristic times for the disturbances, and the signal-to-noise ratio for
each output. For example, you can verify that the following design rejects the
d disturbance almost as well as the optimal Kalman design:
% Alternative estimator design -- output disturbances
taus = [5 5 5];
signoise = [10 10 10];
[Kest, newmod] = smpcest(imod,taus,signoise);
% Simulation using scmpc -- no model error
[y,u,ym] = scmpc(pmod,newmod,ywt,uwt,M,P,tend,...
setpts,ulim,ylim,Kest,z,v,d);
plotall(y,u,dt)
MPC of Nonlinear Plant
We are now ready to test the controller design on the real (nonlinear) plant. A
special version of the scmpc function (called scmpcnl) is available for this
purpose. It uses a nonlinear plant model in the S-function format required by
Simulink. (See the Simulink documentation for more information on how to
write such models.) The model of the paper machine is in the file pap_mach.m.
Simulations with Simulink involving nonlinear models usually take much
longer (by an order of magnitude) than linear simulations of a plant of
comparable complexity. This is especially likely if the plant model is in the form
of an .M file, as is the case here. If such models are to be used extensively, it
may be worthwhile to code them as a .mex file (see MATLAB documentation).
ToseehowwelltheMPCdesignrejectsthed disturbance of Figure 3-8, we
could use the commands found in the file pm_nonl.m in the directory
mpcdemos. The only differences between these commands and those for the
original linear simulation are:
? We have defined the initial values of the plant state and manipulated
variables (x0 and u0, respectively).
? A step size for numerical integration has been specified. The value of 0.05
minutes provides reasonable accuracy in this application. In general, one
mustchoosethestepsizetofittheproblem(oruseavariablestep-size
integration method, as provided by Simulink).
You can verify that the results are nearly identical to those shown in Figure
3-8. In other words, the nonlinearities in the plant have caused negligible
Application: Paper Machine Headbox Control
3-39
performance degradation. Very similar results are also obtained for the
setpoint change of Figure 3-4.
As the magnitude of the disturbance (or setpoint change) increases, nonlinear
effects become significant. For example, Figure 3-9 is for a step in d of 7 units.
If the plant were linear, the curves in Figure 3-9 would be the same shape as
those in Figure 3-8, but scaled by a factor of 7. Although this is approximately
true, there are some qualitative differences. For example, at t = 8 minutes in
Figure 3-9, y
2
has gone below y
1
, whereas in Figure 3-8, y
2
> y
1
at all times.
Figure 3-9 As for Figure 3-8, but With Nonlinear Plant, and Step in d of 7
Units
5101520250
-1
-2
-3
0
3
2
1
4
-4
30
y2
y3
y1
Outputs
Time
Manipulated Variables
-3
-4
-2
2
0
-1
1
51015202503
u2
u1
Time
3 MPC Based on State-Space Models
3-40
If d is increased to 8, control quality degrades dramatically and the maximum
tracking error in y
3
goes to about -10
5
(not shown). This is caused by changes
in the plant characteristics as it moves away from the nominal state (i.e.,
causing errors in the MPC’s linear model).
Sensitivity to modeling error can often be reduced by de-tuning the controller.
A common approach is to increase the magnitudes of theuwtparameters. When
nonlinear effects are severe, however, it may be impossible for any
time-invariant, linear controller to provide stable, offset-free performance. In
that case, if the nonlinear effects are predictable, one might try MPC based on
a nonlinear model (e.g., Gattu and Zafiriou, 1992).
5
Scripts for this purpose can
be developed using the functions in this toolbox.
As a final test, let’s repeat the simulation of Figure 3-8 with a controller
sampling period of 0.25 minutes (recall that the original sampling period was
2 minutes). Results appear in Figure 3-10. Compared to Figure 3-8, which had
no model error (i.e., linear plant), we reduced the disturbance in y
3
by a factor
of 3. Thus, a reduction in sampling period may not lead to robustness problems,
and should be tested more thoroughly. You can verify that it works well for
other combinations of small disturbances and setpoint changes.
5. Gattu, G. and E. Zafiriou, “Nonlinear Quadratic Dynamic Matrix Control with State
Estimation,” Ind. Eng. Chem. Research, 1992, 31, 1096–1104.
Application: Paper Machine Headbox Control
3-41
Figure 3-10 As for Figure 3-8, (d =1) but With Nonlinear Plant, Sampling
Period of 0.25 Minutes
0
-0.05
-0.1
0.05
0.2
0.15
0.1
0.25
-0.15
51015202503
y2
y3
y1
Outputs
Time
Manipulated Variables
Time
51015202503
-0.1
-0.2
0.2
0.1
0
0.3
-0.3
u2
u1
3 MPC Based on State-Space Models
3-42
4
Command Reference
4 Command Reference
4-2
Commands Grouped by Function
Identification
autosc Automatically scales a matrix by its means and standard
deviations.
imp2step Combines MISO impulse response models to form MIMO
models in MPC step format.
mlr Calculates MISO impulse response model via multi-
variable linear regression.
plsr Calculates MISO impulse response model via partial least
squares regression.
rescal Converts scaled data back to its original form.
scal Scales a matrix by specified means and standard
deviations.
validmod Validates a MISO impulse response model using new data.
wrtreg Writes data matrices used for regression.
Plotting and Matrix Information
mpcinfo Outputs matrix type and attributes of system
representation.
plotall Plots outputs and inputs from a simulation run on one
graph.
plotfrsp Plots the frequency response of a system as a Bode plot.
ploteach Makes separate plots of outputs and/or inputs from a
simulation run.
plotstep Plots the coefficients of a model in MPC step form.
Commands Grouped by Function
4-3
Model Conversions
c2dmp Converts state-space model from continuous time to
discrete-time. (Equivalent to c2d in Control System
Toolbox)
cp2dp Converts from a continuous to a discrete transfer function
in poly format.
d2cmp Converts state-space model from discrete-time to
continuous time. (Equivalent to d2c in Control System
Toolbox)
mod2mod Changes sampling period of a model in MPC mod format.
mod2ss Converts a model in MPC mod format to a state-space
model.
mod2step Converts a model in MPC mod format to MPC step
format.
poly2tfd Converts a transfer function in poly format to MPC tf
format.
ss2mod Converts a state-space model to MPC mod format.
ss2step Converts a state-space model to MPC step format.
ss2tf2 Converts state-space model to transfer function.
(Equivalent to ss2tf in Control System Toolbox)
tf2ssm Converts transfer function to state-space model.
(Equivalent to tf2ss in Control System Toolbox)
tfd2mod Converts a model in MPC tf format to MPC mod format.
tfd2step Converts a model in MPC tf format to MPC step format.
th2mod Converts a model in theta format (System Identification
Toolbox) into MPC mod format.
4 Command Reference
4-4
Model Building — MPC mod format
addmd Adds one or more measured disturbances to a plant model.
addmod Combines two models such that the output of one adds to the
input of the other.
addumd Adds one or more unmeasured disturbances to a plant model.
appmod Appends two models in an unconnected, parallel structure.
paramod Puts two models in parallel such that they share a common
output.
sermod Puts two models in series.
Controller Design and Simulation — MPC step format
cmpc Solves the quadratic programming problem to simulate
performance of a closed-loop system with input and output
constraints.
mpccl Creates a model in MPC mod format of a closed-loop
system with an unconstrained MPC controller.
mpccon Calculates the unconstrained controller gain matrix for
MPC.
mpcsim Simulates a closed-loop system with optional saturation
constraints on the manipulated variables.
nlcmpc Simulink S-function block for MPC controller with input
and output constraints (solves quadratic program).
nlmpcsim Simulink S-function block for MPC controller with optional
saturation constraints.
Commands Grouped by Function
4-5
Controller Design and Simulation — MPC mod format
scmpc Solves the quadratic programming problem to simulate
performance of a closed-loop system with input and output
constraints.
smpccl Creates a model in MPC mod format of a closed-loop system
with an unconstrained MPC controller.
smpccon Calculates the unconstrained controller gain matrix for MPC.
smpcest Designs a state estimator for use in MPC.
smpcsim Simulates a closed-loop system with optional saturation
constraints on the manipulated variables.
Analysis
mod2frsp Calculates frequency response for a system in MPC mod
format.
smpcgain Calculates steady-state gain matrix of a system in MPC mod
format.
smpcpole Calculates poles of a system in MPC mod format.
svdfrsp Calculates singular values of a frequency response.
4 Command Reference
4-6
Utility Functions
abcdchkm Checks dimensional consistency of (A,B,C,D) set.
(Equivalent to abcdchk in Control System Toolbox)
dantzgmp Solves quadratic programs.
dareiter Solves discrete Riccati equation by an iterative method.
dimpulsm Generates impulse response of discrete-time system.
(Equivalent to dimpulse in Control System Toolbox)
dlqe2 Calculates state-estimator gain matrix for discrete systems.
dlsimm Simulates discrete-time systems. (Equivalent to dlsim in
Control System Toolbox)
mpcaugss Augments a state-space model with its outputs.
mpcparal Puts two state-space models in parallel.
nargchkm Checks number of M-file arguments. (Equivalent to nargchk
in Control System Toolbox)
mpcstair Creates the stairstep format used to plot manipulated
variables.
vec2mat Converts a vector to a matrix.
addmd
4-7
1addmd
Purpose Adds one or more measured disturbances to a plant model in the MPC mod
format. Used to allow for feedforward compensation in MPC.
Syntax model = addmd(pmod,dmod)
Description The disturbance model contained in dmod adds to the plant model contained in
pmod to form a composite, model, with the structure given in the following block
diagram:
pmod, dmod and model are in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally
create pmod and dmod using either the tfd2mod, ss2mod or th2mod functions.
addmd is a specialized version of paramod. Its main advantage over paramod is
that it assumes all the inputs to dmod aretobemeasureddisturbances.This
saves you the trouble of designating the input types in a separate step.
Example See ss2mod for an example of the use of this function.
Algorithm addmd converts pmod and dmod into their state-space form, then uses the
mpcparal function to build the composite model.
Restrictions ?pmod and dmod must have been created with equal sampling periods and
number of output variables.
?pmod must not include measured disturbances, i.e., its mod format must
specify n
d
=0.
? All inputs to dmod must be classified as manipulated variables. (They will be
reclassified automatically as measured disturbances in model.) So the mod
format of dmod must specify n
d
= n
w
= 0 (which is the default for all model
creation functions).
See Also addmod, addumd, appmod, paramod, sermod
dmod
w
d
model
S
y
+
+
u
y
d
pmod
w
u
addmod
1-8
1addmod
Purpose Combines two models in the MPC mod format such that the output of one
combines with the manipulated inputs of the other. This function is specialized
and rarely needed. Its main purpose is to build up a model of a complex
structure that includes the situation shown in the diagram below.
Syntax pmod = addmod(mod1,mod2)
Description The output(s) of mod2 addtothemanipulatedvariable(s)ofmod1 to form a
composite system, pmod, with the structure given in the following block
diagram:
pmod, mod1 and mod2 are in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally
create mod1 and mod2 using either the tfd2mod, ss2mod or th2mod functions.
The different input types associated withmod1andmod2will be retained in pmod
andwillbeorderedasshowninthediagram.
Example See mod2ss for an example of the use of this function.
Restrictions ?mod1 and mod2 must have been created with equal sampling periods.
? The number of manipulated variables in mod1 must equal the number of
output variables in mod2.
See Also addmd, addumd, appmod, paramod, sermod
mod2
pmod
S
y
+
+
u
1
y
mod1
u
2
d
1
d
2
w
1
w
2
u
1
u
2
d
1
d
2
w
1
w
2
addumd
4-9
1addumd
Purpose Adds one or more unmeasured disturbances to a plant model in MPC mod
format. Used for simulation of disturbances and for design of state estimators
in MPC.
Syntax model = addumd(pmod,dmod)
Description The disturbance model contained in dmod adds to the plant model contained in
pmod to form a composite,model, with the structure given in the following block
diagram:
pmod, dmod and model are in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally
create pmod and dmod using either the tfd2mod, ss2mod or th2mod functions.
addumd is a specialized version of paramod. Its main advantage over paramod is
that it assumes all the inputs to dmod are to be unmeasured disturbances. This
saves you the trouble of designating the input types in a separate step.
Example See ss2mod for an example of the use of this function.
Algorithm addumd converts pmod and dmod into their state-space form, then uses the
mpcparal function to build the composite model.
Restrictions ?pmod and dmod must have been created with equal sampling periods and
number of output variables.
?pmod must not include unmeasured disturbances, i.e., its mod format must
specify n
w
=0.
? All inputs to dmod must be classified as manipulated variables. (They will be
reclassified automatically as unmeasured disturbances in model.) So the
mod format of dmod must specify n
d
= n
w
= 0 (which is the default for all
model creation functions).
See Also addmod, addmd, appmod, paramod, sermod, smpcest
dmod
w
d model
S
y
+
+
u
y
d
pmod
d
u
appmod
1-10
1appmod
Purpose Appends two models to form a composite model that retains the inputs and
outputs of the original models. In other words, for models in the MPC mod
format appmod replaces the append function of the Control Toolbox.
Syntax pmod = appmod(mod1,mod2)
Description The two input models combine as shown in the following block diagram:
mod1, mod2 and pmod are in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally
create mod1 and mod2 using either the tfd2mod, ss2mod,orth2mod function.
Restriction mod1 and mod2 must have been created with equal sampling periods.
See Also addmod, addmd, addumd, paramod, sermod
mod1
y
2
mod2
u
1
d
1
w
1
u
2
d
2
w
2
y
1
y
2
y
1
pmod
u
1
u
2
d
1
d
2
w
1
w
2
autosc, scal, rescal
4-11
1autosc, scal, rescal
Purpose Scales a matrix automatically or by specified mean and standard deviation.
Syntax [ax,mx,stdx] = autosc(x)
sx = scal(x,mx)
sx = scal(x,mx,stdx)
rx = rescal(x,mx)
rx = rescal(x,mx,stdx)
Description autosc scales an input matrix or vector x by its column means (mx)and
standard deviations (stdx) automatically and outputs mx and stdx as options.
By using scal, the input can also be scaled by some specified means and/or
standard deviations. rescal converts scaled data back to original data.
Output mx is a row vector containing the mean value for each column of x while
stdx is a row vector containing the standard deviation for each column.
Outputs ax and sx are obtained by dividing the difference of each column of x
and the mean for the column by the standard deviation for the column, i.e.,
ax(:,i) = (x :,i) – mx(i)/stdx(i).Outputrx is determined by multiplying
each column of x by the corresponding standard deviation and adding the
corresponding mean to that product.
If only two arguments are specified in scal or rescal, x is scaled by specified
means (mx)only.
Example See mlr for an example of the use of these functions.
See Also mlr, plsr, wrtreg
cmpc
1-12
1cmpc
Purpose Simulates closed-loop systems with hard bounds on manipulated variables
and/or outputs using models in the MPC step format. Solves the MPC
optimization problem by quadratic programming.
Syntax yp = cmpc(plant,model,ywt,uwt,M,P,tend,r)
[yp,u,ym]=cmpc(plant,model,ywt,uwt,M,P,tend,...
r,ulim,ylim,tfilter,dplant,dmodel,dstep)
Description
cmpc simulates the performance of the type of system shown in the above
diagram when there are bounds on the manipulated variables and/or outputs.
Measurement noise can be simulated by treating it as an unmeasured
disturbance.
The required input variables are as follows:
plant
Is a model in the MPC step format that represents the plant.
model
Is a model in the MPC step format that is to be used for state estimation in the
controller. In general, it can be different from plant if you want to simulate the
effect of plant/controller model mismatch.
y
Disturbance
Model
Controller
Disturbances
S
u
d
+
+
r
Plant
Plant
Disturbance
Setpoints
Plant
Outputs
d
y
cmpc
4-13
ywt
Is a matrix of weights that will be applied to the setpoint tracking errors. If
ywt=[ ], the default is equal (unity) weighting of all out- puts over the entire
prediction horizon. Ifywt|[ ],itmusthaven
y
columns, where n
y
is the number
of outputs. All weights must be ? 0.
You may vary the weights at each step in the prediction horizon by including
up to Prows inywt. Then the first row of n
y
values applies to the tracking errors
in the first step in the prediction horizon, the next row applies to the next step,
etc. See mpccon for details on the form of the optimization objective function.
If you supply only nrow rows, where 1 ￡ nrow < P, cmpc will use the last row to
fill in any remaining steps. Thus if you want the weighting to be the same for
all P steps, you need only specify a single row.
uwt
Same format as ywt,exceptthatuwt applies to the changes in the manipulated
variables. If you use uwt = [ ], the default is zero weighting. If uwt | [ ],it
must have n
u
columns, where n
u
is the number of manipulated variables.
M
Therearetwowaystospecifythisvariable:
If it is a scalar, cmpc interprets it as the input horizon (number of moves) as in
DMC.
If it is a row vector containing n
b
elements, each element of the vector indicates
the number of steps over which D u = 0 during the optimization and cmpc
interprets it as a set of n
b
blocking factors. There may be 1 ￡ n
b
￡ P blocking
factors, and their sum must be ￡ P
If you set M=[ ] and P | Inf, the default is M=P, which is equivalent to
M=ones(1,P). The default value for M is 1 if P=Inf.
P
The number of sampling periods in the prediction horizon. If P=Inf,the
prediction horizon is infinite.
tend
Is the desired duration of the simulation (in time units).
cmpc
1-14
r
Is a setpoint matrix consisting of N rows and n
y
columns, where n
y
is the
number of output variables, y:
where r
i
(k) is the setpoint for output j at time t=kT,andT is the sampling
period (as specified in the step format of plant and model). If tend > NT,the
setpoints vary for the first N periods in the simulation, as specified by r,and
are then held constant at the values given in the last row of r for the remainder
of the simulation.
In many simulations one wants the setpoints to be constant for the entire time,
in which case r need only contain a single row of n
y
values.
If you set r=[ ],thedefaultisarowofn
y
zeros.
The following input variables are optional. In general, setting one of them
equal to an empty matrix causes cmpc to use the default value, which is given
in the description.
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
cmpc
4-15
ulim
Is a matrix giving the limits on the manipulated variables. Its format is as
follows:
Note that it contains three matrices of N rows. In this case, the limits on N are
1 ￡ N ￡ n
b
,wheren
b
is the number of times the manipulated variables are to
change over the input horizon. If you supply fewer than n
b
rows, the last row
is repeated automatically.
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(2) is the lower bound for manipulated variable j for the
second move of the manipulated variables (where the first move is at the start
of the prediction horizon). If u
min
,
j
(k)=–inf, manipulated variable j will have
no lower bound for that move.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound for that move.
,,
,,
,,
,
()
,
()
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
……
,,
,,
,,
,
()
,
()
ulim
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
cmpc
1-16
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, cmpc will force|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The
limits on the rate of change must be nonnegative and finite.Ifyouwantittobe
unbounded, set the bound to a large number (but not too large — a value of 10
6
should work well in most cases).
The default is u
min
= –inf, u
max
=infand D u
max
=10
6
ylim
Same idea as for ulim, but for the lower and upper bounds of the outputs. The
first row applies to the first point in the prediction horizon. The default is y
min
= –inf,andy
max
=inf.
tfilter
Is a matrix of time constants for the noise filter and the unmeasured
disturbances entering at the plant output. The first row of n
y
elements gives
the noise filter time constants and the second row of n
y
elements gives the time
constants of the lags through which the unmeasured disturbance steps pass. If
tfilter only contains one row, the unmeasured disturbances are assumed to
be steps. If you set tfilter=[ ] or omit it, the default is no noise filtering and
steplike unmeasured disturbances.
dplant
Is a model in MPC step format representing all the disturbances (measured and
unmeasured) that affect plant in the above diagram. If dplant is provided,
then input dstep is also required. For output step disturbances, set
dplant=[ ]. The default is no disturbances.
dmodel
Is a model in MPC step format representing the measured disturbances. If
dmodel is provided, then input dstep is also required. If there are no measured
disturbances, set dmodel=[ ]. For output step disturbances, set dmodel=[ ].If
there are both measured and un- measured disturbances, set the columns of
dmodel corresponding to the unmeasured disturbances to zero. The default is no
measured disturbances.
cmpc
4-17
dstep
Is a matrix of disturbances to the plant. For output step disturbances
(dplant=[ ] and dmodel=[ ]), the format is the same as for r. For disturbances
through step-response models (dplant only or both dplant and dmodel
nonempty), the format is the same as for r, except that the number of columns
is n
d
rather than n
y
. The default is a row of zeros.
Notes ? You may use a different number of rows in the matrices r, ulim, ylim and
dstep, should that be appropriate for your simulation.
? The ulim constraints used here are fundamentally different from the usat
constraints used in the mpcsim function. The ulim constraints are defined
relative to the beginning of the prediction horizon, which moves as the
simulation progresses. Thus at each sampling period, k,theulim constraints
apply to a block of calculated moves that begin at sampling period k and
extend for the duration of the input horizon. The usat constraints, on the
otherhand,arerelativetothefixedpointt = 0, the start of the simulation.
The calculated outputs are as follows (all but yp are optional):
yp
Is a matrix containing M rows and n
y
columns, where M = max(fix(tend=T)+1,
2). The first row will contain the initial condition, and row k – 1 will give the
values of the plant outputs, y (see above diagram), at time t=kT.
u
Is a matrix containing the same number of rows asypand n
u
columns. The time
corresponding to each row is the same as for yp. The elements in each row are
the values of the manipulated variables, u (see above diagram).
ym
Is a matrix of the same structure as yp, containing the values of the predicted
output from the state estimator in the controller. These will, in general, differ
from those in yp if model|plant and/or there are unmeasured disturbances.
The prediction includes the effect of the most recent measurement, i.e., it is
.
For unconstrained problems, cmpc and mpcsim should give the same results.
The latter will be faster because it uses an analytical solution of the QP
problem, whereas cmpc solves it by iteration.
y? kk()
cmpc
1-18
Examples Consider the linear system:
The following statements build the model and set up the controller in the same
way as in the mpcsim example.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2; tfinal=90;
model=tfd2step(tfinal,delt,ny,g11,g21,g12,g22);
plant=model;
P=6; M=2; ywt=[ ]; uwt=[1 1];
tend=30; r=[0 1];
Here, however, we will demonstrate the effect of constraints. First we set a
limitof0.1ontherateofchangeofu
1
and a minimum of –0.15 for u
2
.
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[ ];
[y,u]=cmpc(plant,model,ywt,uwt,M,P,tend,r,ulim,ylim);
plotall(y,u,delt),pause
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
cmpc
4-19
Note that D u
2
has a large (but finite) limit. It never comes into play.
We next apply a lower bound of zero to both outputs:
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[0 0 inf inf];
[y,u]=cmpc(plant,model,ywt,uwt,M,P,tend,r,ulim,ylim);
plotall(y,u,delt),pause
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.2
?0.18
?0.16
?0.14
?0.12
?0.1
?0.08
?0.06
Manipulated Variables
Time
cmpc
1-20
The following results show that no constraints are violated.
Restriction Initial conditions of zero are used for all the variables. This simulates the
condition where all variables represent a deviation from a steady-state initial
condition.
Suggestion Problems with many inequality constraints can be very time consuming. You
can minimize the number of constraints by:
? Using small values for P and/or M.
? Leaving variables unconstrained (limits at ±inf) intermittently unless you
think the constraint is important.
See Also plotall, ploteach, mpccl, mpccon, mpcsim
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.2
?0.1
0
Manipulated Variables
Time
cp2dp
4-21
1cp2dp
Purpose Converts a single-input-single-output, continuous-time transfer function in
standard MATLAB polynomial form (including an optional time delay) to a
sampled-data transfer function.
Syntax [numd,dend] = cp2dp(num,den,delt)
[numd,dend] = cp2dp(num,den,delt,delay)
Description num and den are the numerator and denominator polynomials of the
continuous-time system (in the standard Control Toolbox polynomial format),
delt is the sampling period, and delay is the (optional) time delay (in time
units). If you omit delay, cp2dp assumes zero delay. The calculated results are
numd and dend, the numerator and denominator polynomials of the
corresponding discrete-time transfer function. cp2dp adds a zero-order hold at
the input of the continuous-time system during the conversion to discrete-time.
cp2dp accounts properly for the effect of a time delay that is a nonintegral
multiple of the sampling period. If delay=0, cp2dp is equivalent to the
MATLAB commands:
[a,b,c,d]=tf2ss(num,den);
[phi,gam]=c2dmp(a,b,delt);
[numd,dend]=ss2tf2(phi,gam,c,d,1);
Example See poly2tfd, poly format for an example of the use of this function.
Algorithm cp2dp first converts num and den to the equivalent discrete state-space form. It
then accounts for the fractional time delay (if any) using the formulas in
?str?m and Wittenmark (1984), pages 40–42. Finally, it converts the discrete
state-space model to a discrete transfer-function model, simultaneously
accounting for the whole periods of delay (if any).
Reference ?str?m,K.J.;Wittenmark,B.Computer Control Systems Theory and Design,
Prentice-Hall, Englewood Cliffs, N.J., 1984.
Restriction The order of num must be ￡ that of den.
See Also poly2tfd, poly format
dlqe2
1-22
1dlqe2
Purpose Solves the discrete Riccati equation by an iterative method to determine the
optimal steady-state gain (and optional covariance matrices) for a discrete
Kalman filter or state estimator.
Syntax k = dlqe2(phi,gamw,c,q,r)
[k,m,p] = dlqe2(phi,gamw,c,q,r)
Description Filter form:
Consider the state-space description:
where x is a vector of n state variables, u contains n
u
known inputs, is a
vector of n
y
measured outputs, y is the noise-free output, w is a vector of n
w
unmeasured disturbance inputs, z is a vector of n
y
measurement noise inputs,
and F , G
u
, G
w
,C and D are constant matrices. We assume that w and z are
stationary random-normal signals (white noise) with covariances
E{w(k)w
T
(k)} = Q
E{w(k)z
T
(k)} = R
12
=0
E{z(k)z
T
(k)} = R
The steady-state Kalman filter is
where is the estimate of x(k) based on the measurements available at
period k, is that based on the measurements available at period
xk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() zk()+=
Cx k() Du k() zk()++=
y
x? kk()x? kk 1–()Kyk() Cx? kk 1–()Du k()––[]+=
x? k 1+ k()Fx? kk()G
u
uk()+=
y? kk()Cx? kk()Du k()+=
x? kk()
x? kk 1–()
dlqe2
4-23
k – 1, etc. Note that is an estimate of the noise-free output, . The
steady-state Kalman gain, K,isthesolutionof
where M and P may be interpreted as the expected covariance of the errors in
the state estimates before and after the measurement update, respectively, i.e.,
where, by definition,
The dlqe2 function takes F , G
u
, C, R,andQ as inputs and calculates K, M,and
P. The last two output arguments are optional.
Note that the input and output arguments are identical to those for dlqe in the
Control Toolbox. The advantage of dlqe2 is that it can handle a singular
state-transition matrix (F ), e.g., for systems with time delay.
Predictor form:
y? y
KMC
T
RCMC
T
+[]
1–
=
M F PF
T
G
w
QG
w
T
+=
PMKCM–=
MEx? kk 1–()x? kk 1–()
T
{}=
PEx? kk()x? kk()
T
{}=
x? kk()xk() x? kk()–=
x? kk 1–()xk() x? kk 1–()–=
dlqe2
1-24
You can also use dlqe2 to calculate a state-estimator in the predictor form:
The relationship between K
p
, the estimator gain for the predictor form, and K
as calculated by dlqe2 is:
The matrix M calculated by dlqe2 is the expected covariance of the errors in
.
Algorithm dlqe2 callsdareiter
1
, which solves the discrete algebraic Riccati equation
using an iterative doubling algorithm.
Example Consider a system represented by the block diagram:
x? k 1+ k()Fx? kk 1–()G
u
uk() K
p
ek()++=
y? kk 1–()Cx? kk 1–()Du k()+=
ek() yk() y? kk 1–()–=
K
p
F K=
x? kk 1–()
1. We gratefully acknowledge Kjell Gustafsson, Department of Automatic Control, Lund
Institute of Technology, Lund, Sweden, who provided this function.
ySu
+
G
w
G
u
z
y
w
S
+
+
+
dlqe2
4-25
where G
u
and G
w
are first-order, discrete-time transfer functions.
and the statistics of the unmeasured inputs are Q =2,R =1.
We use the appropriate MPC Toolbox functions to build a model of the system,
then calculate the optimal gain:
delt=2; ny=1;
gu=poly2tfd(0.2,[1 -0.8],delt);
Gw=poly2tfd(0.3,[1 -0.95],delt);
[phi,gam,c,d]=mod2ss(tfd2mod(delt,ny,gu,Gw));
k=dlqe2(phi,gam(:,2),c,2,1)
The result is:
k = 0
1.0619
Note that the gain for the first state is zero since this corresponds to the state
of G
u
, which is unaffected by the disturbance, w.Alsonoticethatinthe
composite system, the second column of gam is G
w
. This is because of the order
in which gu and G
w
were specified as inputs to the tfd2mod function.
See Also smpcest
G
u
z()
0.20
10.8z
1–
–
--------------------------= G
w
z()
0.3
10.95z
1–
–
-----------------------------=
imp2step
1-26
1imp2step
Purpose Constructs a multi-input multi-output model in MPC step format from
multi-input single-output impulse response matrices.
Syntax plant = imp2step(delt,nout,theta1,theta2, ...,
theta25)
Description Given the impulse response coefficient matrices, theta1, theta2,etc.,amodel
in MPC step format is constructed. Each thetai is an n-by-n
u
matrix
corresponding to the impulse response coefficients for output i. n is the number
of the coefficients and n
u
is the number of inputs.
delt is the sampling interval used for obtaining the impulse response
coefficients. nout is the output stability indicator. For stable systems, this
argument is set equal to number of outputs, n
y
. For systems with one or more
integrating outputs, this argument is a column vector of length n
y
with
nout(i)=0 indicating an integrating output and nout(i)=1 indicating a stable
output.
Example See mlr and plsr for examples of the use of this function.
Restriction The limit on the number of impulse response matrices thetai is 25.
See Also mlr, plsr
mlr
4-27
1mlr
Purpose Determines impulse response coefficients for a multi-input single-output
system via Multivariable Least Squares Regression or Ridge Regression.
Syntax [theta, yres] = mlr(xreg,yreg,ninput)
[theta, yres] = mlr(xreg,yreg,ninput,plotopt,wtheta, ...
wdeltheta)
Description xreg and yreg are the input matrix and output vector produced by routines
such aswrtreg.ninputis number of inputs. Least Squares is used to determine
the impulse response coefficient matrix, theta. Columns of theta correspond
to impulse response coefficients from each input. Optional output yres is the
vector of residuals, the difference between the actual outputs and the predicted
outputs.
Optional inputs include plotopt, wtheta,andwdeltheta. No plot is produced
if plotopt is equal to 0 which is the default; a plot of the actual output and the
predicted output is produced if plotopt=1; two plots — plot of actual and
predicted output, and plot of residuals — are produced forplotopt=2.Penalties
on the squares of theta and the changes in theta can be specified through the
scalar weights wtheta and wdeltheta, respectively (defaults are 0). theta is
calculated as follows:
theta1 = (X
T
X)
–1
X
T
Y
where
X
xreg
wtheta I·
wdeltheta delI·
=
Y
yreg
0
0
=
…
mlr
1-28
where I is identity matrix of dimension n
*
n
u
dimension of delI is n
*
n
u
by n
*
n
u
.
then
Example Consider the following two-input single-output system:
Load the input and output data. The input and output data were generated
from the above transfer function and random zero-mean noise was added to the
output. Sampling time of 7 minutes was used.
load mlrdat;
Determine the standard deviations for input data using the function autosc.
[ax,mx,stdx] = autosc(x);
Scale the input data by their standard deviations only.
mx = [0,0];
sx = scal(x,mx,stdx);
Put the input and output data in a form such that they can be used to
determine the impulse response coefficients. 35 impulse response coefficients
(n)areused.
delI
1– 10… 0
01– 1 … 0
0 … 01– 1
0 … 00 1–
=
…
theta [theta11:n()theta1 n 1:2n+()…=
theta1 nu n 1–()1:nu n+ ]
**
ys()
5.72e
14s–
60s 1+
-------------------------
1.52e
15s–
25s 1+
-------------------------
u
1
s()
u
2
s()
=
mlr
4-29
n = 35;
[xreg,yreg] = wrtreg(sx,y,n);
Determine the impulse response coefficients via mlr.Byspecifyingplotopt=2,
two plots — plot of predicted output and actual output, and plot of the output
residual (or predicted error) — are produced.
ninput = 2;
plotopt = 2;
[theta,yres] = mlr(xreg,yreg,ninput,plotopt);
Scale theta based on the standard deviations used in scaling the input.
theta = scal(theta,mx,stdx);
Convert the impulse model to a step model to be used in MPC design. Recall
that a sampling time of 7 minutes was used in determining the impulse model.
Number of outputs (1 in this case) must be specified.
nout = 1;
delt = 7;
model = imp2step(delt,nout,theta);
0 20 40 60 80 100 120 140 160 180
?6
?4
?2
0
2
4
Actual value (o) versus Predicted Value (+)
Sample Number
Output
0 20 40 60 80 100 120 140 160 180
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
Output Residual or Prediction Error
Sample Number
Residual
mlr
1-30
Plot the step response coefficients.
plotstep(model)
See Also plsr, validmod, wrtreg
0 50 100 150 200 250
0
1
2
3
4
5
6
u1 step response : y1
TIME
0 50 100 150 200 250
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
u2 step response : y1
TIME
mod format
4-31
1mod format
Purpose The MPC mod format is a compact way to store the model of a linear system
for subsequent use in the MPC Toolbox functions.
Description
Consider the process shown in the above block diagram. Its discrete-time LTI
state-space representation is:
where x is a vector of n state variables, u represents the n
u
manipulated
variables, d represents n
d
measured but freely-varying inputs (i.e., measured
disturbances), w represents n
w
unmeasured disturbances, y is a vector of n
y
plant outputs, z is measurement noise, and F , G
u
, etc., are constant matrices of
appropriate size. The variable (k) represents the plant output before the
addition of measurement noise. Define:
G =[G
u
G
d
G
w
]
D =[D
u
D
ud
D
w
]
In some cases one would like to include n
ym
measured and n
yu
unmeasured
outputs in y, where n
ym
+ n
yu
= n
y
.Ifso,themod format assumes that the y
vector and the C and D matrices are arranged such that the measured outputs
come first, followed by the unmeasured outputs.
S
y y
Unmeasured
Disturbances
Manipulated
Variables
Measured
Disturbances
Measurement
Noise
w
u
d
++
z
Plant
Measured
Outputs
xk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() zk()+=
Cx k() D
u
uk() D
d
dk() D
w
wk() zk()+++ +=
y
mod format
1-32
The mod format is a single matrix that contains the F, G, C,andD matrices,
plus some additional information. Let M be the mod representation of the
above system. Its overall dimensions are:
? Number of rows = n + n
y
+1
? Number of columns = max(7, 1 + n + n
u
+ n
d
+ n
w
)
The minfo vector is the first seven elements of the first row in M. The elements
of minfo are:
The remainder of M contains the discrete state-space matrices:
Notes Since the minfo vector requires seven columns, this is the minimum possible
number of columns in the mod format, regardless of the dimensions of the
state-space matrices.
Also, the first column is reserved for other uses by MPC Toolbox routines. Thus
the state-space matrices start in column 2, as described above.
In order for the mpcinfo routine to recognize matrices in the MPC mod format,
the (2,1) element is set to NaN (Not-a-Number).
minfo (1) T, the sampling period used to create the model.
(2) n, the number of states.
(3) n
u
, the number of manipulated variable inputs.
(4) n
d
, the number of measured disturbances.
(5) n
w
, the number of unmeasured disturbances.
(6) n
ym
, the number of measured outputs.
(7) n
yu
, the number of unmeasured outputs.
F in rows 2 to n +1 columns2ton +1
G in rows 2 to n +1 columnsn +2ton + n
u
+n
d
+ n
w
+1
C in rows n +2ton + n
y
+1 columns2ton +1
D in rows n +2ton + n
y
+1 columnsn +2ton + n
u
+n
d
+ n
w
+1
mod format
4-33
Example See ss2mod for a mod format example.
See Also mod2ss, mod2step, step format, mpcinfo, ss2mod, step, tfd2mod, tf format,
th2mod, theta format
mod2frsp, varying format
1-34
1mod2frsp, varying format
Purpose Calculates the complex frequency response in varying format of a system in
MPC mod format.
Syntax frsp = mod2frsp(mod,freq)
[frsp,eyefrsp] = mod2frsp(mod,freq,out,in,balflg)
Description mod2frsp calculates the complex frequency response of a system (mod)inMPC
mod format. The desired frequencies are given by the input freq,arow vector
of 3 elements specifying the lower frequency as a power of 10, the upper
frequency as a power of 10, and the number of frequency points.
Optional inputs out and in are row vectors that specify the outputs and inputs
for which the frequency response is to be generated. If these variables are
omitted or empty, the default is to use all outputs and inputs.
Optional input balflg indicates whether the system’s F matrix should be
balanced (using the MATLAB balance command). If balflg is nonzero,
balancing is performed. Balancing improves the conditioning of the problem,
but may cause errors in the frequency response. If balflg=[ ] or is omitted, no
balancing is performed.
Outputfrsp is the frequency response matrix given in varying format. Let F(w )
denote a matrix-valued function of the independent variable w . Then the N
sampled values F(w
1
),...,F(w
N
) are contained in frsp as follows:
If the dimension of each submatrix F(w
i
)isn by m, then the dimensions offrsp
is n
.
N +1bym +1.
frsp
F w
1
() w
1
F w
i
() w
N
0
F w
N
()
0
0…0Ninf
=
… …
…
…
mod2frsp, varying format
4-35
Optional output eyefrsp is in varying format and represents I – F(w
i
)ateach
frequency. This output can only be specified for square submatrices and may
be useful in computing the frequency responses of both the sensitivity and
complementary sensitivity functions.
Example Consider the linear system:
See the mpccl example for the commands that build the a closed-loop model for
this process using a simple controller. However for this example, delt=6 and
tfinal=90 are used to reduce the number of step response coefficients.
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
mod2frsp, varying format
1-36
Now we will calculate and plot the frequency response of the sensitivity and
complementary sensitivity functions.
freq = [-3,0,30];
in = [1:ny]; % input is r for comp. sensitivity
out = [1:ny]; % output is yp for comp. sensitivity
[frsp,eyefrsp] = mod2frsp(clmod,freq,out,in);
plotfrsp(eyefrsp); % Sensitivity
pause;
plotfrsp(frsp); % Complementary Sensitivity pause;
10
?3
10
?2
10
?1
10
0
10
?4
10
?2
10
0
10
2
Frequency (radians/time)
Log Magnitude
BODE PLOTS
10
?3
10
?2
10
?1
10
0
?500
?400
?300
?200
?100
0
100
Frequency (radians/time)
Phase (degrees)
mod2frsp, varying format
4-37
10
?3
10
?2
10
?1
10
0
10
?4
10
?2
10
0
10
2
Frequency (radians/time)
Log Magnitude
BODE PLOTS
10
?3
10
?2
10
?1
10
0
?700
?600
?500
?400
?300
?200
?100
0
Frequency (radians/time)
Phase (degrees)
mod2frsp, varying format
1-38
Calculate and plot the singular values for the sensitivity function response.
[sigma,omega] = svdfrsp(eyefrsp);
clg;
semilogx(omega,sigma);
title(￠ Singular Values vs. Frequency￠ );
xlabel(￠ Frequency (radians/time)￠ );
ylabel(￠ Singular Values￠ );
Algorithm The algorithm to calculate the complex frequency response involves a matrix
inverse problem which is solved via a Hessenberg matrix.
Reference A.J. Laub, “Efficient Multivariable Frequency Response Computations,” IEEE
Transactions on Automatic Control, Vol. AC–26, No. 2, pp.407–408, April,
1981.
See Also mod, mpccl, plotfrsp, smpccl, svdfrsp
10
?3
10
?2
10
?1
10
0
0
0.5
1
1.5
2
2.5
Singular Values vs. Frequency
Frequency (radians/time)
Singular Values
mod2mod
4-39
1mod2mod
Purpose Changes the sampling period of a model in MPC mod format.
Syntax newmod = mod2mod(oldmod,delt2)
Description Input oldmod is the existing model in MPC mod format. Inputdelt2 is the new
sampling period for the model. mod2mod returnsnewmod, which is the system in
mod format converted to the new sampling time.
See Also mod, ss2mod
mod2ss
1-40
1mod2ss
Purpose Extracts the standard discrete-time state-space matrices and other
information from a model stored in the MPC mod format.
Syntax [phi,gam,c,d] = mod2ss(mod)
[phi,gam,c,d,minfo] = mod2ss(mod)
Description
Consider the process shown in the above block diagram. mod2ss assumes that
mod is a description of the above process in the MPC mod format (see modin the
online MATLAB Function Reference for more details). An equivalent
state-space representation is:
where x is a vector of n state variables, u represents the n
u
manipulated
variables, d represents n
d
measured but freely-varying inputs (i.e., measured
disturbances), w represents n
w
unmeasured disturbances, y is a vector of n
y
plant outputs, z is measurement noise, and F , G
u
, etc., are constant matrices of
appropriate size. The variable (k) represents the plant output before the
addition of measurement noise. Define:
G =[G
u
G
d
G
w
]
D =[D
u
D
d
D
w
]
S
y y
Unmeasured
Disturbances
Manipulated
Variables
Measured
Disturbances
Measurement
Noise
w
u
d
++
z
Plant
Measured
Outputs
xk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() zk()+=
Cx k() D
u
uk() D
d
dk() D
w
wk() zk()+++ +=
y
mod2ss
4-41
mod2ssextracts the F , G , C,andD matrices from the input variable, mod.Italso
extracts the vector minfo, which contains additional information about the
sampling period, number of each type of input and output, etc. see mod in the
online MATLAB Function Reference for more details on minfo.
Examples 1 See the example in the description of dlqe2.
2 Suppose you have a plant with the structure
where the inputs and outputs are all scalars, and you have constructed mod1
and mod2 using the commands:
phi1=diag([-0.7, 0.8]); gam1=[1, -1, 0; 0, 0, 1];
c1=[0.2 -0.4]; d1=zeros(1,3);
minfo1=[1 2 1 1 1 1 0];
mod1=ss2mod(phi1,gam1,c1,d1,minfo1);
phi2=-0.2; gam2=[1, -0.5, 0.2];
c2=3; d2=[0.2, 0, 0];
minfo2=[1 1 1 1 1 1 0];
mod2=ss2mod(phi2,gam2,c2,d2,minfo2);
pmod=addmod(mod1,mod2);
Now you want to calculate the response to a step change in d
2
,whichisthe
fourth input to the composite system, pmod.Onewaytodoitis:
[phi,gam,c,d,minfo]=mod2ss(pmod);
nstep=10;
ustep=[zeros(nstep,3) ones(nstep,1) zeros(nstep,2)];
% Define step in d2
y=dlsimm(phi,gam,c,d,ustep);
% simulate response to step input
plot([0:nstep-1],y)
mod2
S
y(k)
+
+
mod1
u
1
u
2
d
1
d
2
w
1
w
2
mod2ss
1-42
The results of the mod2ss call are:
phi=
-0.7000 0 3.0000
0 0.8000 0
0 0 -2.0000
gam=
1.0000 2.2000 -1.0000 0 0 0
0
0 0 0 0 1.0000 0
0 1.0000 0 -.5000 0 0.2000
c=
0.2000 -0.4000 0
d=
0 0 0 0 0 0
minfo=
1 3 2 2 2 1 0
And the step response is:
See Also mod, ss2mod
0 1 2 3 4 5 6 7 8 9
?0.35
?0.3
?0.25
?0.2
?0.15
?0.1
?0.05
0
mod2step, step format
4-43
1mod2step, step format
Purpose Uses a model in the mod format to calculate the step response of a SISO or
MIMO system in MPC step format.
Syntax plant = mod2step(mod,tfinal)
[plant,dplant] = mod2step(mod,tfinal,delt2,nout)
Description The input variable mod is assumed to be a model in the mod format (see mod in
the online MATLAB Function Reference for a description). You would normally
create it using ss2mod, tfd2mod,orth2mod. The input variable tfinal is the
timeatwhichyouwouldliketoendthestepresponse.
The optional input variable delt2 is the desired sampling period for the step
response. If you use delt2=[ ] or omit it, the default is equal to the sampling
period of mod (contained in the minfo vector of mod).
The optional input variable nout is the output stability indicator. For stable
systems, set nout equal to the number of outputs, n
y
.Forsystemswithoneor
more integrating outputs, nout is a column vector of length n
y
with nout(i)=0
indicating an integrating output and nout(i)=1 indicating a stable output. If
you use nout=[ ] or omit it, the default is nout=n
y
(only stable outputs).
plant and dplant are matrices in MPC step format containing the calculated
step responses.plantistheresponsetothemanipulatedvariables,anddplant
is the response to the disturbances (if any), both measured and unmeasured.
The overall dimensions of these matrices are:
plant n-by-n
y
+ n
y
+2rows,n
u
columns.
dplant n-by-n
y
+ n
y
+2rows,n
d
+ n
w
columns.
where n = round (tfinal/delt2)
It is assumed that stable step responses are nearly constant after n sampling
periods, while integrating responses increase with a constant slope after n–1
sampling periods.
Each column gives the step response with respect to the corresponding input
variable. Within each column, the first n
y
elements are the response for each
output at time t=T, the next n
y
elements give each output at time t=2T,etc.
mod2step, step format
1-44
The last n
y
+ 2 rows contain nout, n
y
and delt2, respectively (all in column 1
— any remaining elements in these rows are set to zero). In other words, for
plant thearrangementisasfollows:
where
S
kji
is the i
th
step response coefficient describing the effect of input j on output
k.
The arrangement of dplant is similar; the only difference is in the number of
columns.
plant
S
1
S
2
S
n
nout(1) 0 … 0
nout(2) 0 … 0
nout n
y
() 0 … 0
n
y
0 … 0
delt2 0 … 0
nn
y
n
y
2++ ()n
u
·
=
…
…
… …
S
i
S
11i,,
S
12i,,
… S
1 n
u
i,,
S
21i,,
S
22i,,
… S
2 n
u
i,,
S
n
y
1 i,,
S
n
y
2 i,,
… S
n
y
n
u
i,,
=
…
mod2step, step format
4-45
Example The following process has 3 inputs and 4 outputs:
phi=diag([0.3, 0.7, -0.7]);
gam=eye(3);
c=[1 0 0; 0 0 1; 0 1 1; 0 1 0];
d=[1 0 0; zeros(3,3)];
We first calculate its step response for 4 samples (including the initial
condition) with respect to each of the inputs using the Control Toolbox function,
dstep:
nstep=4; delt=1.5;
yu1=dstep(phi,gam,c,d,1,nstep)
yu2=dstep(phi,gam,c,d,2,nstep)
yu3=dstep(phi,gam,c,d,3,nstep)
The results are:
We then use mod2step to do the same job:
plant=mod2step(ss2mod(phi,gam,c,d,delt),(nstep-1)*delt)
Response to u
1
Response to u
2
Response to u
3
Time y
1
y
2
y
3
y
4
y
1
y
2
y
3
y
4
y
1
y
2
y
3
y
4
0100000000000
T200000110110
2T 2.3 0 0 0 0 0 1.7 1.7 0 0.3 0.3 0
3T 2.39 0 0 0 0 0 2.19 2.19 0 0.79 0.79 0
mod2step, step format
1-46
obtaining the results:
plant =
2.0000 0 0
0 0 1.0000
0 1.0000 1.0000
0 1.0000 0
2.3000 0 0
0 0 0.3000
0 1.7000 0.3000
0 1.7000 0
2.3900 0 0
0 0 0.7900
0 2.1900 0.7900
0 2.1900 0
1.0000 0 0
1.0000 0 0
1.0000 0 0
1.0000 0 0
4.0000 0 0
1.5000 0 0
See Also plotstep, ss2step, tfd2step
mpcaugss
4-47
1mpcaugss
Purpose Differences the states of a system and augments them with the output
variables. Mainly used as a utility function for setting up model predictive
controllers.
Syntax [phia,gama,ca] = mpcaugss(phi,gam,c)
[phia,gama,ca,da,na] = mpcaugss(phi,gam,c,d)
Description
Consider the process shown in the above block diagram. A state-space
representation is:
where x is a vector of n state variables, u is a vector of n
u
manipulated
variables, d is a vector of n
d
measured disturbances, w is a vector of n
w
unmeasured disturbances, y is a vector of n
y
plant outputs, z is measurement
noise, and F , G
u
, G
d
, G
w
, etc., are constant matrices of appropriate size. The
variable = Cx(k) represents the plant output before the addition of the
direct contribution of the inputs [D
u
u(k)+D
d
v(k)+D
w
w(k)] and the
measurement noise [z(k)]. (The variable is the output before addition of the
measurement noise). Define:
D u(k)=u(k)–u(k –1)
D x(k)=x(k)–x(k –1)
S
y y
Unmeasured
Disturbances
Manipulated
Variables
Measured
Disturbances
Measurement
Noise
w
u
d
++
z
Plant
Measured
Outputs
zk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() D
u
uk() D
d
dk() D
w
wk() zk()+++ +=
yk() zk()+=
y? k()
y
mpcaugss
1-48
etc. Then equations 4.28 and 4.29 can be converted to the form
x
a
(k+1) = F
a
x
a
(k)+G
ua
D u(k)+G
da
D (k)+G
wa
D w(k)
y(k)=C
a
x
a
(k)+D
u
u(k)+D
d
d(k)+D
w
w(k)+z(k)
where, by definition,
The mpcaugss function takes the matrices F , G (= [G
u
G
d
G
w
]), C as input, and
creates the augmented matrices F
a
, G
a
, C
a
and D
a
in the form shown above.
The D input matrix is optional. If you include it, mpcaugss assumes it has the
form D =[D
u
D
d
D
w
]. If you omit it, the default is zero. Note that all MPC
design routines require D
u
= D
d
=0.
The last output variable, na, is the order of the augmented system, i.e.,
n
a
= n + n
y
. It is optional.
Example The following system has 2 states, 3 inputs, and 2 outputs.
phi=diag([0.8, -0.2]);
gam=[1 -1 0;0 2 -0.5];
c=[0.4 0;0 1.5];
x
a
k()
D xk()
y
?
k()
=
F
a
F 0
CF I
G
a
G
ua
G
da
G
wa
==
G
ua
G
u
CG
u
G
da
G
d
CG
d
G
wa
G
w
CG
w
===
C
a
0 I
D
a
D
u
D
d
D
w
==
mpcaugss
4-49
Here is the augmentation command, followed by the calculated results:
[phia,gama,ca,da,na]=mpcaugss(phi,gam,c)
phia =
0.8000 0 0 0
0 -0.2000 0 0
0.3200 0 1.0000 0
0 -0.3000 0 1.0000
gama =
1.0000 -1.0000 0
0 2.0000 -0.5000
0.4000 -0.4000 0
0 3.0000 -0.7500
ca =
0 0 1 0
0 0 0 1
da =
0 0 0
0 0 0
na =
4
mpccl
1-50
1mpccl
Purpose Combines a plant model and a controller model in MPC step format, yielding a
closed-loop system model in the MPC mod format. This can be used for stability
analysis and linear simulations of closed-loop performance.
Syntax [clmod] = mpccl(plant,model,Kmpc)
[clmod,cmod] = mpccl(plant,model,Kmpc,tfilter,...
dplant, dmodel)
Description
plant
Is a model (in step format) representing the plant in the above diagram.
model
Is a model (in step format)thatistobeusedtodesigntheMPCcontrollerblock
shown in the diagram. It may be the same as plant (in which case there is no
“model error” in the controller design), or it may be different.
Kmpc
Is a controller gain matrix, which was calculated by the function mpccon.
Disturbance
Model
Controller
Disturbances
S
u
d
+
+
r
Plant
Plant
Disturbance
Setpoints
Noise-free
d
y
Plant
Outputs
S
Measured
Outputs
S
w
u
y
p
Measurement
Noise
+
+
+
+
z
mpccl
4-51
tfilter
Is a (optional) matrix of time constants for the noise filter and the unmeasured
disturbances entering at the plant output. If omitted or set to an empty matrix,
the default is no noise filtering and steplike unmeasured disturbances. See the
documentation for the function mpcsim for more details on the design and
proper format of tfilter.
dplant
Is a (optional) model (in step format) representing all the disturbances
(measured and unmeasured) that affect plant in the above diagram. If omitted
or set to an empty matrix, the default is that there are no disturbances.
dmodel
Is a (optional) model (in step format) representing the measured disturbances.
If omitted or set to an empty matrix, the default is that there are no measured
disturbances. See the documentation for the function mpcsim for more details
on how disturbances are handled when using step-response models.
mpccl
Calculates a model of the closed-loop system, clmod.Itisinthemod format and
can be used, for example, with analysis functions such as smpcgain and
smpcpole, and with simulation routines such as mod2step and dlsimm. mpccl
also calculates (as an option) a model of the controller element, cmod.
The closed-loop model, clmod, has the following state-space representation:
x
cl
(k +1)=F
cl
x
cl
(k)+G
cl
u
cl
(k)
y
cl
(k)=C
cl
x
cl
(k)+D
cl
u
cl
(k)
where x
cl
is a vector of n state variables, u
cl
is a vector of input variables, y
cl
is
a vector of outputs, and F
cl
, G
cl
, C
cl
,andD
cl
are matrices of appropriate size.
Theexpertusermaywanttoknowthesignificanceofthestatevariablesinx
cl
.
They are (in the following order):
? The n
p
states of the plant (as specified in plant),
? The n
i
state estimates (based on the model specified in model),
? n
d
integrators that operate on the D d signal to yield a d signal. If there are
no disturbances, these states are omitted.
mpccl
1-52
? n
u
integrators that operate on the D w
u
signal to yield a w
u
signal.
? n
u
integrators that operate on the D u signal produced by the standard MPC
formulation to yield a u signal that can be used as input to the plant and as
a closed-loop output.
The closed-loop input and output variables are:
where is the estimate of the noise-free plant output at sampling period
k based on information available at period k. This estimate is generated by the
controller element.
The state-space form of the controller model, cmod,canbewrittenas:
x
c
(k +1)=F
c
x
c
(k)+G
c
u
c
(k)
y
c
(k)=C
c
x
c
(k)+D
c
u
c
(k)
where
and the controller states are the same as those of the closed loop system except
that the n
p
plantstatesarenotincluded.
Example Consider the linear system:
u
cl
k()
rk()
zk()
w
u
k()
dk()
and y
cl
k()
y
p
k()
uk()
y? kk()
==
y? kk()
u
c
k()
rk()
yk()
dk 1–()
and y
c
k() uk 1–()==
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
mpccl
4-53
We build the step response model using the MPC Toolbox functions poly2tfd
and tfd2step.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2; tfinal = 60;
model=tfd2step(tfinal,delt,ny,g11,g21,g12,g22);
plant=model; % No plant/model mismatch
Now we design the controller. Since there is delay, we use M < P:Wespecify
the defaults for the other tuning parameters, uwt and ywt, then calculate the
controller gain:
P=6; % Prediction horizon.
M=2; % Number of moves (input horizon).
ywt=[ ]; % Output weights (default - unity
% on all outputs). uwt=[ ]; % Man. Var weights (default - zero
% on all man. vars).
Kmpc=mpccon(model,ywt,uwt,M,P);
Now we can calculate the model of the closed-loop system:
clmod=mpccl(plant,model,Kmpc);
You can use the closed-loop model to calculate and plot the step response with
respecttoalltheinputs.Theappropriatecommandsare:
tend=30;
clstep=mod2step(clmod, tend);
plotstep(clstep)
mpccl
1-54
Since the closed-loop system has m =6inputsandp = 6 outputs, only one of the
plots is reproduced here. It shows the response of the first 4 closed-loop outputs
to a unit step in the first closed-loop input, which is the setpoint for y
1
:
Closed-loop outputs y
1
and y
2
are the true plant outputs (noise-free). Output y
1
goes to the new setpoint quickly with a small overshoot. This causes a small,
short-term disturbance in y
2
. The plots for y
3
and y
4
show the required
variation in the manipulated variables.
Restriction model and plant must have been created using the same sampling period.
See Also cmpc, mod2step, step format, mpccon, mpcsim, smpcgain, smpcpole
0 10 20 30
0
0.2
0.4
0.6
0.8
1
1.2
1.4
u1 step response :y1
TIME
0 10 20 30
-0.2
-0.15
-0.1
-0.05
0
u1 step response :y2
TIME
0 10 20 30
0
0.05
0.1
0.15
0.2
0.25
u1 step response :y3
TIME
0 10 20 30
0
0.02
0.04
0.06
0.08
0.1
0.12
u1 step response :y4
TIME
mpccon
4-55
1mpccon
Purpose Calculates MPC controller gain using a model in MPC step format.
Syntax Kmpc = mpccon(model)
Kmpc = mpccon(model,ywt,uwt,M,P)
Description Combines the following variables (most of which are optional and have default
values) to calculate the MPC gain matrix, Kmpc.
model
is the model of the process to be used in the controller design (in the step
format).
The following input variables are optional:
ywt
Is a matrix of weights that will be applied to the setpoint tracking errors. If you
use ywt=[ ] or omit it, the default is equal (unity) weighting of all outputs over
theentirepredictionhorizon.Ifyouuseywt|[ ],itmusthaven
y
columns,
where n
y
is the number of outputs. All weights must be ? 0.
You may vary the weights at each step in the prediction horizon by including
up to Prows inywt. Then the first row of n
y
values applies to the tracking errors
in the first step in the prediction horizon, the next row applies to the next step,
etc.
If you supply only nrow rows, where 1 ￡ nrow < P, mpccon will use the last row
to fill in any remaining steps. Thus if you wish the weighting to be the same for
all P steps, you need only specify a single row.
uwt
Same format as ywt,exceptthatuwt applies to the changes in the manipulated
variables. If you use uwt = [ ] or omit it, the default is zero weighting. If uwt
|[],itmusthaven
u
columns, where n
u
is the number of manipulated
variables.
M
Therearetwowaystospecifythisvariable:
mpccon
1-56
? If it is a scalar, mpccon interprets it as the input horizon (number of moves)
as in DMC.
? If it is a row vector containing n
b
elements, each element of the vector
indicates the number of steps over which D u = 0 during the optimization and
cmpc interprets it as a set of n
b
blocking factors. There may be 1 ￡ n
b
￡ P
blocking factors, and their sum must be ￡ P.
If you set M=[ ]or omit it and P | Inf, the default is M=P, which is equivalent
to M=ones(1,P). The default value for M is 1 if P=Inf.
P
The number of sampling periods in the prediction horizon. If you set P=Inf or
omit it, the default is P=1.IfP=inf, the prediction horizon is infinite.
If you take the default values for all the optional variables, you get the “perfect
controller,” i.e., a model-inverse controller. This controller is not applicable
when one or more outputs can not respond to the manipulated variables within
one sampling period due to time delay. In this case, the plant-inverse controller
is unrealizable. For nonminimum phase discrete plants, this controller is
unstable. To counteract this you can penalize changes in the manipulated
variables (variable uwt), use blocking (variable M), and/or make P>>M.The
model-inverse controller is also relatively sensitive to model error and is best
used as a point of reference from which you can progress to a more robust
design.
Algorithm The controller gain is a component of the solution to the optimization problem:
with respect to (a series of current and future moves in the manipulated
variables), where (k + j) is a prediction of output i at a time j sampling
periods into the future (relative to the current time, k), which is a function of
Minimize Jk() ywt
i
j()r
i
kj+()y?
i
– kj+()[]
2
i 1=
n
y
Ge5
j 1=
p
Ge5
=
uwt
i
j()D u?
i
j()()
2
i 1=
n
u
Ge5
j 1=
n
b
Ge5
+
D u?
i
j()
y?
i
mpccon
4-57
, r
i
(k + j) is the corresponding future setpoint, and n
b
is the number of
blocks or moves of the manipulated variables.
Example Consider the linear system:
See the mpccl exampleforthecommandsthatbuildthemodelandasimple
controller for this process.
Here is a slightly more complex design with blocking and time-varying weights
on the manipulated and output variables:
P=6; M=[2 4];
uwt=[1 0; 0 1];
ywt=[1 0.1; 0.8 0.1; 0.1 0.1];
Kmpc=mpccon(model,ywt,uwt,M,P);
tend=30; r=[1 0];
[y,u]=mpcsim(plant,model,Kmpc,tend,r);
There is no particular rationale for using time varying weights in this case it is
only for illustration. The manipulated variables will make 2 moves during the
prediction horizon (see value of M,above).Theuwt selection gives u
1
a unity
weight and u
2
a zero weight for the first move, then switches the weights for
the second move. If there had been any additional moves they would have had
thesameweightingasthesecondmove.
D u?
i
j()
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
mpccon
1-58
The ywt value assigns a constant weight of 0.1 to y
2
, and a weight that
decreases over the first 3 periods to y
1
. The weights for periods 4 to 6 are the
same as for period 3. The resulting closed-loop (servo) response is:
See Also cmpc, mpccl, mpcsim
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
Manipulated Variables
Time
mpcinfo
4-59
1mpcinfo
Purpose Determines the type of a matrix and returns information about the matrix.
Syntax: mpcinfo(mat)
Description mpcinfo returns information about the type and size of the matrix, mat.The
information is determined from the matrix structure. The matrix types include
MPC step format, MPC mod format, varying format and constant. mpcinfo
returns text output to the screen.
If the matrix is in MPC step format, the output includes the sampling time used
to create the model, number of inputs, number of outputs and number of step
response coefficients; it also indicates which outputs are stable and which
outputs are integrating.
If the matrix is in MPC mod format, the output includes the sampling time
used to create the model, number of states, number of manipulated variable
inputs, number of measured disturbances, number of unmeasured
disturbances, number of measured outputs and number of unmeasured
outputs.
For a matrix in varying format, as formed in mod2frsp, the number of
independent variable values, and the number of rows and number of columns
of each submatrix are output.
For a constant matrix, the text output consists of the number of rows and
number of columns.
Examples 1 MPC step format: After running the mod2step example mpcinfo(plant)
returns:
This is a matrix in MPC Step format.
sampling time = 1.5
number of inputs = 3
number of outputs = 4
number of step response coefficients = 3
All outputs are stable.
mpcinfo
1-60
2 MPC mod format: After running the ss2mod example mpcinfo(pmod)
returns:
This is a matrix in MPC Mod format.
minfo = [2 3 1 1 1 1 0 ]
sampling time = 2
number of states = 3
number of manipulated variable inputs = 1
number of measured disturbances = 1
number of unmeasured disturbances = 1
number of measured outputs = 1
number of unmeasured outputs = 0
3 varying format: After running the mod2frsp example mpcinfo(eyefrsp)
returns:
varying: 30 pts 2 rows 2 cols
See Also mod, step, mod2frsp, varying format
mpcsim
4-61
1mpcsim
Purpose Simulates closed-loop systems with saturation constraints on the manipulated
variables using models in the MPC step format. Can also be used for open-loop
simulations.
Syntax yp = mpcsim(plant,model,Kmpc,tend,r)
[yp,u,ym] = mpcsim(plant,model,Kmpc,tend,r,usat,...
tfilter,dplant,dmodel,dstep)
Description
mpcsim provides a convenient way to simulate the performance of the type of
system shown in the above diagram. Measurement noise can be simulated by
treating it as an unmeasured disturbance. The required input variables are as
follows:
plant
Is a model in the MPC step format that is to represent the plant.
model
Is a model in the MPC step format that is to be used for state estimation in the
controller. In general, it can be different from plant if you want to simulate the
effect of plant/controller model mismatch. Note, however, that model should be
thesameasthatusedtocalculateKmpc.
Kmpc
Is the MPC controller gain matrix, usually calculated using the function
mpccon.
y
Disturbance
Model
Controller
Disturbances
S
u
d
+
+
r
Plant
Plant
Disturbance
Setpoints
Plant
Outputs
d
y
mpcsim
1-62
If you set Kmpc to an empty matrix, mpcsim will do an open-loop simulation.
Then the inputs to the plant will be r (which must be set to the vector of
manipulated variables in this case) and dstep.
tend
Is the desired duration of the simulation (in time units).
r
Is normally a setpoint matrix consisting of N rows and n
y
columns, where n
y
is
the number of output variables, y:
where r
i
(k) is the setpoint for output j at time t=kT,andT is the sampling
period (as specified in the step format of plant and model). If tend > NT,the
setpoints vary for the first N periods in the simulation, as specified by r,and
are then held constant at the values given in the last row of r for the remainder
of the simulation.
In many simulations one wants the setpoints to be constant for the entire time,
in which case r need only contain a single row of n
y
values.
If you set r=[ ],thedefaultisarowofn
y
zeros.
For open-loop simulations, r specifies the manipulated variables and must
contain n
u
columns.
The following input variables are optional. In general, setting one of them
equal to an empty matrix causes mpcsimto use the default value, which is given
in the description.
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
mpcsim
4-63
usat
Is a matrix giving the limits on the manipulated variables. Its format is as
follows:
Note that it contains three matrices of N rows. N may be different than that for
the setpoint matrix, r, but the idea is the same: the saturation limits will vary
for the first N sampling periods of the simulation, then be held constant at the
values given in the last row of usat for the remaining periods (if any).
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(k) is the lower bound for manipulated variable j at time t=
kT in the simulation. If u
min
,
j
(k)=–inf, manipulated variable j will have no
lower bound at t=kT.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound at t=kT.
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
,,
,,
,,
,
()
,
()
usat
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
,,
,,
,,
,
()
,
()
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
……
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
mpcsim
1-64
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, mpcsim will force
|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The limits on the rate of change must be
nonnegative.
The default is no saturation constraints, i.e., all the u
min
values will be set to –
inf, and all the u
max
and D u
max
values will be set to inf.
Note: Saturation constraints in mpcsim are enforced by simply clipping the
manipulated variable moves so that they satisfy all constraints. This is a
nonoptimal solution that, in general, will differ from the results you would get
using the ulim variable in cmpc.
tfilter
Is a matrix of time constants for the noise filter and the unmeasured
disturbances entering at the plant output. The first row of n
y
elements gives
the noise filter time constants and the second row of n
y
elements gives the time
constants of the lags through which the unmeasured disturbance steps pass. If
tfilter only contains one row, the unmeasured disturbances are assumed to
be steps. If you set tfilter=[ ] or omit it, the default is no noise filtering and
steplike unmeasured disturbances.
dplant
Is a model in MPC step format representing all the disturbances (measured and
unmeasured) that affect plant in the above diagram. If dplant is provided,
then input dstep is also required. For output step disturbances, set
dplant=[ ]. The default is no disturbances.
dmodel
Is a model in MPC step format representing the measured disturbances. If
dmodel is provided, then input dstep is also required. If there are no measured
disturbances, set dmodel=[ ]. For output step disturbances, set dmodel=[ ].If
there are both measured and un- measured disturbances, set the columns of
dmodel corresponding to the unmeasured disturbances to zero. The default is no
measured disturbances.
mpcsim
4-65
dstep
Is a matrix of disturbances to the plant. For output step disturbances
(dplant=[ ] and dmodel=[ ]), the format is the same as for r. For disturbances
through step-response models (dplant only or both dplant and dmodel
nonempty), the format is the same as for r, except that the number of columns
is n
d
rather than n
y
. The default is a row of zeros.
Note: You may use a different number of rows in the matrices r, usat and
dstep, should that be appropriate for your simulation.
The calculated outputs are as follows (all but yp are optional):
yp
Is a matrix containing M rows and n
y
columns, where M = round(tend/delt2)
+ 1 and delt2 is the sampling time. The first row will contain the initial
condition, and row k – 1 will give the values of the plant outputs, y (see above
diagram), at time t=kT.
u
Is a matrix containing the same number of rows as yp and n
u
columns. The
time corresponding to each row is the same as for yp. The elements in each row
are the values of the manipulated variables, u (see above diagram).
ym
Is a matrix of the same structure as yp, containing the values of the predicted
outputs from the state estimator in the controller. ym will, in general, differ
from yp if model|plant and/or there are unmeasured disturbances. The
prediction includes the effect of the most recent measurement, i.e., .y? kk()
mpcsim
1-66
Examples Consider the linear system:
The following statements build the model and calculate the MPC controller
gain:
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2; tfinal=90;
model=tfd2step(tfinal,delt,ny,g11,g21,g12,g22);
plant=model;
P=6; M=2;
ywt=[ ]; uwt=[1 1];
Kmpc=mpccon(imod,ywt,uwt,M,P);
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
mpcsim
4-67
Simulate and plot the closed-loop performance for a unit step in the setpoint for
y
2
, occurring at t =0.
tend=30; r=[0 1];
[y,u]=mpcsim(plant,model,Kmpc,tend,r);
plotall(y,u,delt),pause
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.24
?0.22
?0.2
?0.18
?0.16
?0.14
?0.12
?0.1
Manipulated Variables
Time
mpcsim
1-68
Try a pulse change in the disturbance that adds to u
1
:
r=[ ]; usat=[ ]; tfilter=[ ]; dmodel=[ ];
dplant=plant;
dstep=[ 1 0; 0 0];
[y,u]=mpcsim(plant,model,Kmpc,tend,r,usat,tfilter,...
dplant,dmodel,dstep);
plotall(y,u,delt),pause
For the same disturbance as in the previous case, limit the rates of change of
both manipulated variables.
usat=[-inf -inf inf inf 0.1 0.05];
[y,u]=mpcsim(plant,model,Kmpc,tend,r,usat,tfilter,...
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
1.5
Outputs
Time
0 5 10 15 20 25 30
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
0.1
Manipulated Variables
Time
mpcsim
4-69
dplant,dmodel,dstep);
plotall(y,u,delt),pause
Restriction Initial conditions of zero are used for all the variables. This simulates the
condition where all variables represent a deviation from a steady-state initial
condition.
See Also plotall, ploteach, cmpc, mpccl, mpccon
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
1.5
2
Outputs
Time
0 5 10 15 20 25 30
?0.3
?0.2
?0.1
0
0.1
Manipulated Variables
Time
nlcmpc
1-70
1nlcmpc
Purpose Model predictive controller for simulating closed-loop systems with hard
bounds on manipulated variables and/or controlled variables using linear
models in the MPC step format for nonlinear plants represented as Simulink
S-functions.
Description nlcmpc is a Simulink S-function block and can be invoked by typing nlmpclib
at the MATLAB prompt. Its usage is identical to other Simulink blocks. The
input to nlcmpc includes both the variables controlled by nlcmpc and measured
disturbances. The first n
y
elements of the input are treated as the controlled
variables while the rest is taken as the measured disturbances. The output
from nlcmpc are the values of the manipulated variables. Initial conditions for
the manipulated variables and the measured disturbances must be specified.
The controlled variables sent to nlcmpc and the manipulated variables
returned by nlcmpc are actual variables; they are not deviation variables.
Because of the limit on the number of masked variables that can be specified
for a Simulink block, model and dmodel are put together as “one” variable, r,
ywt,anduwt as “one” variable, and ylim and ulim as “one” variable. m and p
should be entered as one row vector. u0 and d0 should also be entered as one
row vector. The required input variables are as follows:
modelpd
Equals [model dmodel]. model is a linear model in the MPC step format that is
to be used for state estimation in the controller. In general, it is a linear
approximation of the nonlinear plant. dmodel is a model in MPC step format
representing the effect of the measured disturbances. The default is no
measured disturbances. Note that the truncation time for model and dmodel
must be the same and the number of outputs for model and dmodel must be the
same.
nlcmpc
4-71
ryuwt
Equals [r ywt uwt]. r is a setpoint matrix consisting of N rows and n
y
columns, where n
y
is the number of output variables, y:
Where r
i
(k) is the setpoint for output i at time t=kT,andT is the sampling
period (as specified in the step format of model). If the simulation time is larger
than NT, the setpoints vary for the first N periods in the simulation, as
specified by r, and are then held constant at the values given in the last row of
r for the remainder of the simulation. In many simulations one wants the
setpointstobeconstantfortheentiretime,inwhichcaser need only contain
asinglerowofn
y
values.
ywt
Must have n
y
columns, where n
y
is the number of outputs. All weights must be
? 0.
You may vary the weights at each step in the prediction horizon by including
up to Prows inywt. Then the first row of n
y
values applies to the tracking errors
in the first step in the prediction horizon, the next row applies to the next step,
etc. See mpccon for details on the form of the optimization objective function.
If you supply only nrow rows, where 1 ￡ nrow < P, nlcmpc will use the last row
to fill in any remaining steps. Thus if you wish the weighting to be the same for
all P steps, you need only specify a single row.
uwt
Has the same format as ywt,exceptthatuwt applies to the changes in the
manipulated variables. If uwt|[ ],itmusthaven
u
columns, where n
u
is the
number of manipulated variables.
Notice that the number of rows for r, ywt,anduwt should be the same. If not,
onecanenterthevariableasparpart(r, ywt, uwt). The function parpart
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
nlcmpc
1-72
appends extra rows to r, ywt,and/oruwt so that they have the same number
of rows. The default is r=y0, where y0 is the initial condition for the output,
equal (unity) weighting of all outputs over the entire prediction horizon and
zero weighting of all input.
mp
Equals [M P]. P equals the last element of MP. There are two ways to specify M:
If it is a scalar, nlcmpc interprets it as the input horizon (number of moves) as
in DMC; if it is a row vector containing n
b
elements, each element of the vector
indicates number of the steps over which D u(k) = 0 during the optimization and
nlcmpc interprets it as a set of n
b
blocking factors. There may be 1 ￡ n
b
￡ P
blocking factors, and their sum must be ￡ P.IfyousetM=[ ],thedefaultis
M = P, which is equivalent toM=ones(1,P). Pis the number of sampling periods
in the prediction horizon.
yulim
Equals [ylim ulim]. ulim is a matrix giving the limits on the manipulated
variables. Its format is as follows:
,,
,,
,,
,
()
,
()
,,
,,
,,
,
()
,
()
ulim
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
……
nlcmpc
4-73
Note that it contains three matrices of N rows. In this case, the limits on N are
1 ￡ N ￡ n
b
,wheren
b
is the number of times the manipulated variables are to
change over the input horizon. If you supply fewer than n
b
rows, the last row is
repeated automatically.
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(2) is the lower bound for manipulated variable j for the
second move of the manipulated variables (where the first move is at the start
of the prediction horizon). If u
min
,
j
(k)=–inf, manipulated variable j will have
no lower bound for that move.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound for that move.
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, cmpc will force|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The
limits on the rate of change must be nonnegative and finite.Ifyouwantittobe
unbounded, set the bound to a large number (but not too large — a value of 10
6
should work well in most cases).
ylim has the same format as ulim, but for the lower and upper bounds of the
outputs. The first row applies to the first point in the prediction horizon.
Note that the number of rows for ylim and ulim should be the same. If the
number of rows for ylim and ulim differs, one can use parpart(ylim, ulim).
The function parpart appends extra rows to ylim or ulim so that they have the
same number of rows. If you set yulim = [ ],thenu
min
= –inf, u
max
=inf, D u
max
=10
6
, y
min
= –inf and y
max
=inf.
tfilter
Is a matrix of time constants for the noise filter and the unmeasured
disturbances entering at the plant output. The first row of n
y
elements gives
the noise filter time constants and the second row of n
y
elements gives the time
constants of the lags through which the unmeasured disturbance steps pass. If
tfilter only contains one row, the unmeasured disturbances are assumed to
be steps. If you set tfilter= [ ], no noise filtering and steplike unmeasured
disturbances are assumed.
nlcmpc
1-74
ud0
Equals [u0 d0]. u0 are initial values of the manipulated variables arranged in
a row vector having n
u
elements; n
u
is the number of the manipulated variables
computed by nlcmpc. d0 are initial values of the measured disturbances
arranged in a row vector having n
d
elements; n
d
is the number of the measured
disturbances. The default is u0 = 0 and d0 = 0.
Notes ? Initial conditions for the manipulated variables that are calculated by
nlcmpc are specified through nlcmpc while initial conditions for the
controlled variables are specified through the S-function for the nonlinear
plant.
? You may use a different number of rows in the matrices r, ulim and ylim,
should that be appropriate for your simulation.
? The ulim constraints used here are fundamentally different from the usat
constraints used in the nlmpcsim block. The ulim constraints are defined
relative to the beginning of the prediction horizon, which moves as the
simulation progresses. Thus at each sampling period, k,theulim constraints
apply to a block of calculated moves that begin at sampling period k and
extend for the duration of the input horizon. The usat constraints, on the
otherhand,arerelativetothefixedpointt = 0, the start of the simulation.
? For unconstrained problems, nlcmpc and nlmpcsim should give the same
results. The latter will be faster because it uses an analytical solution of the
QP problem, whereas nlcmpc solves it by iteration.
Example See the examples for nlmpcsim with one modification: replace the block
nlmpcsim with nlcmpc. Clearly, additional variables should be defined
appropriately.
See Also cmpc, nlmpcsim
nlmpcsim
4-75
1nlmpcsim
Purpose Model predictive controller for simulating closed-loop systems with saturation
constraints on the manipulated variables using linear models in the MPC step
format for nonlinear plants represented as Simulink S-functions.
Description nlmpcsim is a Simulink S-function block and can be invoked by typing
nlmpclib at the MATLAB prompt. Its usage is identical to other Simulink
blocks. The input to nlmpcsim includes both the variables controlled by
nlmpcsim and measured disturbances. The first n
y
elements of the input are
treated as the controlled variables while the rest is taken as the measured
disturbances. The output from nlmpcsim arethevaluesofthemanipulated
variables. Initial conditions for the manipulated variables and the measured
disturbances must be specified. Both the controlled variables sent to nlmpcsim
and the manipulated variables returned by nlmpcsim are the actual variables;
they are not deviation variables.
Because of the limit on the number of masked variables that can be specified
for a Simulink block, model and dmodel are put together as one variable. u0 and
d0 should be entered as one row vector. The required input variables are as
follows:
modelpd
Equals [model dmodel]. model is a linear model in the MPC step format that is
to be used for state estimation in the controller. In general, it is a linear
approximation for the nonlinear plant. Note, however, thatmodelshould be the
same as that used to calculate Kmpc. dmodel is a model in MPC step format
representing the measured disturbances. If dmodel = [ ], the default is no
measured disturbances. Note that the truncation time for model and dmodel
should be the same and the number of outputs for model and dmodel should be
the same.
Kmpc
Is the MPC controller gain matrix, usually calculated using the function
mpccon.
nlmpcsim
1-76
r
Is a setpoint matrix consisting of N rows and n
y
columns, where n
y
is the
number of controlled variables, y:
Where r
i
(k) is the setpoint for output i at time t=kT,andT is the sampling
period (as specified in the step format of model). If the simulation time is larger
than NT, the setpoints vary for the first N periods in the simulation, as
specified by r, and are then held constant at the values given in the last row of
r fortheremainderofthesimulation.
In many simulations one wants the setpoints to be constant for the entire time,
in which case r need only contain a single row of n
y
values.
Note that r is the actual setpoint. If you set r=[ ], the default is y0.
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
nlmpcsim
4-77
usat
Is a matrix giving the saturation limits on the manipulated variables. Its
format is as follows:
Note that it contains three matrices of N rows. N may be different from that for
the setpoint matrix, r, but the idea is the same: the saturation limits will vary
for the first N sampling periods of the simulation, then be held constant at the
values given in the last row of usat for the remaining periods (if any).
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(k) is the lower bound for manipulated variable j at time t=
kT in the simulation. If u
min
,
j
(k)=–inf, manipulated variable j will have no
lower bound at t=kT.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound at t=kT.
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
,,
,,
,,
,
()
,
()
,,
,,
,,
,
()
,
()
usat
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
……
nlmpcsim
1-78
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, mpcsim will force|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The
limits on the rate of change must be nonnegative.
If usat= [ ],thenalltheu
min
values will be set to –inf, and all the u
max
and
u
max
values will be set to inf.
Note: Saturation constraints are enforced by simply clipping the manipulated
variable moves so that they satisfy all constraints. This is a nonoptimal
solution that, in general, will differ from the results you would get using the
ulim variable in cmpc or nlcmpc.
tfilter
Is a matrix of time constants for the noise filter and the unmeasured
disturbances entering at the plant output. The first row of n
y
elements gives
the noise filter time constants and the second row of n
y
elements gives the time
constants of the lags through which the unmeasured disturbance steps pass. If
tfilter only contains one row, the unmeasured disturbances are assumed to
be steps. If you set tfilter= [ ], no noise filtering and steplike unmeasured
disturbances are assumed.
ud0
Equals [u0 d0]. u0 are initial values of the manipulated variables arranged in
a row vector having n
u
elements; n
u
is the number of the manipulated variables
computed by nlmpcsim. d0 are initial values of the measured disturbances
arranged in a row vector having n
d
elements; n
d
is the number of the measured
disturbances. The default is u0 = 0 and d0 = 0.
Note: You may use a different number of rows in the matrices r and usat,
should that be appropriate for your simulation.
Examples Let us now demonstrate the use of the controller nlmpcsim. Since the plant
used in Example 1 is linear, using mpcsim would be much faster. The point,
however, is to show how masked variables are specified for nlmpcsim.
nlmpcsim
4-79
1 The plant is linear with two inputs and two outputs. It is represented by
The Simulink S-function for this plant is in mpcplant.m. The nominal
steady-state operating condition is y0 = [58.3 1.5] and u0 = [100 1].The
Simulink block to simulate this plant using nlmpcsim is in nlmpcdm1.m and
showninFigure1-1.
Figure 1-1 Simulink Block for Example 1
The following statements build the step response model and specify the
parameter values. Note that model does not equal the plant model stored in
dx
dt
------
1.2– 0
0
1
1.5
--------–
x
0.2
1
u
50
0
++=
yx=
y1, y2
nlmpcsim Demo #1
(Double click on the "?" for more info)
nlmpcsim
nlmpcsim
t
Save Time
y
Save Outputs
u
Save Inputs
mpcplant
Sfunction
Plot
Results
?
Load
Data
Clock
nlmpcsim
1-80
mpcplant.m. The important thing to notice is that both r and usat are actual
variables. They are not deviation variables.
g11=poly(0.4,[1 2]);
g21=poly2tfd(0,1);
g12=poly2tfd(0,1);
g22=poly2tfd(1,[1 1]);
tfinal=8;
delt=0.2;
nout=2;
model=tfd2step(tfinal,delt,nout,g11,g21,g12,g22);
ywt=[1 1];
uwt=[0 0];
M=4;
P=10;
r=[68.3 2];
usat=[100 1 200 3 200 200];
tfilter=[ ];
Kmpc = mpccon(model,ywt,uwt,M,P);
dmodel = [ ];
Therearetwowaystosimulatetheclosedloopsystem.Wecansetthe
simulation parameters and click on Start under Simulation or via the
following statements.
plant=￠ nlmpcdm1￠ ; y0=[58.3 1.5];
u0=[100 1];
tfsim = 2;
tol=[1e-3];
minstep=[ ];
maxstep=[ ];
[t,yu]=gear(plant,tfsim,[y0 u0],[tol,minstep,maxstep]);
Figure1-2showstheresponseforthesetpointchange.
nlmpcsim
4-81
Figure 1-2 Output responses for a setpoint change for Example 1
2 The plant is the paper machine headbox discussed in the section,
“Application: Paper Machine Headbox Control” in Chapter 3. The nonlinear
plant model is represented as a Simulink S-function and is in pap_mach.m.
The plant has two inputs, three outputs, four states, one measured
disturbance, and one unmeasured disturbance. All these variables are zero
at the nominal steady-state. Since the model for nlmpcsim must be linear,
we linearize the nonlinear plant at the nominal steady-state to obtain a
linear model. Since the model is simple, we can linearize it analytically to
obtain A, B, C, and D.
The Simulink block to simulate this nonlinear plant using nlmpcsim is in
nlmpcdm2.m and shown in Figure 1-3.
nlmpcsim
1-82
Figure 1-3 Simulink Block for Example 2
The following statements build the step response model and specify the
parameter values.
A=[-1.93 0 0 0; .394 -.426 0 0; 0 0 -.63 0; .82 -.784
.413 -.426];
B=[1.274 1.274 0; 0 0 0; 1.34 -.65 .203; 0 0 0];
C=[0 1 0 0; 0 0 1 0; 0 0 0 1];
D=zeros(3,3);
% Discretize the linear model and save in MOD form.
dt=2;
[PHI,GAM]=c2dmp(A,B,dt);
minfo=[dt,4,2,1,0,3,0];
imod=ss2mod(PHI,GAM,C,D,minfo);
% Store plant model and measured disturbance model in MPC
% step format
Gp, Gw
Np
Nw
nlmpcsim Demo #2
(Double click on the "?" for more info)
nlmpcsim
nlmpcsim
Step Fcn1
Step Fcn
t
Save Time
y
Save Outputs
u
Save Inputs
papmach
Sfunction
Plot
Results
Mux
Mux2
Mux
Mux1
?
Load
Data
Clock
nlmpcsim
4-83
[model,dmodel]=mod2step(imod,20);
m=5;
p=20;
ywt=[1 0 5]; % unequal weighting of y1 and y3, no control
% of y2
uwt=[1 1]; % Equal weighting of u1 and u2
ulim=[-10 -10 10 10 2 2]; % Constraints on u
ylim=[ ]; % No constraints on y
usat=ulim;
tfilter=[ ];
y0=[0 0 0];
u0=[0 0];
r=[0 0 0];
Kmpc=mpccon(model,ywt,uwt,M,P);
Figure 1-4 shows the output responses for a unit-step measured disturbance
Np = 1 and a step unmeasured disturbance with Nw =5.
Figure 1-4 Output responses for a unit-step measured disturbance Np =1 and
a step unmeasured disturbance Nw =5
See Also mpcsim, nlcmpc
-6
-4
2
0
-2
4
5 1015202503
Ou
tpu
t
Time
Solid: y1 (H2)
Dashed: y2 (N1)
Dotted: y3 (N2)
y1
y2
y3
paramod
1-84
1paramod
Purpose Puts two models in parallel by connecting their outputs. Mimics the utility
function mpcparal,exceptthatparamod worksonmodelsintheMPCmod
format.
Syntax pmod = paramod(mod1,mod2)
Description
mod1 and mod2 are models in the MPC mod format (see mod in the online
MATLAB Function Reference format section for a detailed description). You
would normally create them using either the tfd2mod, ss2mod or th2mod
functions.
paramod combines them to form a composite system, pmod,asshowninthe
above diagram. It is also in the mod format.Notehowtheinputstomod1 and
mod2 are ordered in pmod.
Restriction mod1 and mod2 must have been created with equal sampling periods and they
must have the same number of measured and unmeasured outputs.
See Also addmd, addmod, addumd, appmod, sermod
mod1
y
2
mod2
u
1
d
1
w
1
u
2
d
2
w
2
y
1
y
pmod
u
1
u
2
d
1
d
2
w
1
w
2
y
S
+
+
plotall
4-85
1plotall
Purpose Plots outputs and manipulated variables from a simulation, all on one “page.”
Syntax plotall(y,u) plotall(y,u,t)
Description Input variables y and u are matrices of outputs and manipulated variables,
respectively. Each row represents a sample at a particular time. Each column
shows how a particular output (or manipulated) variable changes with time.
Input variable t is optional. If you supply it as a scalar, plotall interprets is
as the sampling period, and calculates the time axis for the plots accordingly.
It can also be a column vector, in which case it must have the same number of
rows as y and u and is interpreted as the times at which the samples of y and
u were taken. If you do not supply t, plotall uses a sampling period of 1 by
default.
plotall plots all the outputs on a single graph. If there are multiple outputs
that have very different numerical scales, this may be unsatisfactory. In that
case, use ploteach.
plotall plots all the manipulated variables in “stairstep” form (i.e., assuming
a zero-order hold) on a single graph. Again, ploteach may be the better choice
if scales are very different.
plotall
1-86
Example output: (mpccon example)
See Also ploteach, plotstep, plotfrsp
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
Manipulated Variables
Time
ploteach
4-87
1ploteach
Purpose Plots outputs and manipulated variables from a simulation on separate
graphs, up to four per page.
Syntax ploteach(y)
ploteach(y,u)
ploteach([ ],u)
ploteach(y,[ ],t)
ploteach([ ],u,t)
ploteach(y,u,t)
Description Input variables y and u are matrices of outputs and manipulated variables,
respectively. Each row represents a sample at a particular time. Each column
shows how a particular output (or manipulated) variable changes with time. As
shown above, you may supply both y and u,oromiteitheroneofthem.
Input variable t is optional. If you supply it as a scalar, ploteach interprets is
as the sampling period, and calculates the time axis for the plots accordingly.
It can also be a column vector, in which case it must have the same number of
rows as y and u and is interpreted as the times at which the samples of y and
u were taken. If you do not supply t, ploteach uses a sampling period of 1 by
default.
ploteach plots the manipulated variables in “stairstep” form (i.e., assuming a
zero-order hold).
ploteach
1-88
Example output: (mpccon example)
See Also plotall, plotfrsp, plotstep
plotfrsp
4-89
1plotfrsp
Purpose Plots the frequency response generated by mod2frsp as a Bode plot.
Syntax plotfrsp(vmat)
plotfrsp(vmat,out,in)
Description vmat is a varying matrix which contains the data to be plotted.
Let F(w ) denote the matrix (whose entries are functions of the independent
variable w )whosesampledvaluesF(w
1
), ... , F(w
N
) are contained in vmat.
plotfrsp(vmat) will generate Bode plots of all elements of F(w ).
Optional inputs out and in are row vectors which specify the row and column
indices respectively of a submatrix of F(w ). plotfrsp will then generate Bode
plots of the elements of the specified submatrix of F(w ).
Example Output: (mod2frsp example)
See Also mod2frsp, varying format
10
?3
10
?2
10
?1
10
0
10
?4
10
?2
10
0
10
2
Frequency (radians/time)
Log Magnitude
BODE PLOTS
10
?3
10
?2
10
?1
10
0
?500
?400
?300
?200
?100
0
100
Frequency (radians/time)
Phase (degrees)
plotstep
1-90
1plotstep
Purpose Plots multiple step responses as calculated by mod2step, ss2step or tfd2step.
Syntax plotstep(plant)
plotstep(plant,opt)
Description plant is a step-response matrix in the MPC step format created by mod2step,
ss2step or tfd2step.
opt is an optional scalar or row vector that allows you to select the outputs to
be plotted. If you omit opt, plotstep plots every output. If, for example, plant
contains results for 6 outputs, setting opt=[1,4,5] would cause only y
1
, y
4
and
y
5
to be plotted.
Example output: (tfd2step example)
0 10 20 30 40 50 60 70 80 90
0
2
4
6
8
10
12
14
u1 step response : y1
TIME
0 10 20 30 40 50 60 70 80 90
0
1
2
3
4
5
6
7
u1 step response : y2
TIME
plotstep
4-91
See Also imp2step, mod2step, step format, plotall, ploteach, plotfrsp, ss2step,
tfd2step
0 10 20 30 40 50 60 70 80 90
?20
?15
?10
?5
0
u2 step response : y1
TIME
0 10 20 30 40 50 60 70 80 90
?20
?15
?10
?5
0
u2 step response : y2
TIME
plsr
1-92
1plsr
Purpose Determine the impulse response coefficients for a multi-input single-output
system via Partial Least Squares (PLS).
Syntax [theta,w,cw,ssqdif,yres] = plsr(xreg,yreg,ninput,lv)
[theta,w,cw,ssqdif,yres] = plsr(xreg,yreg,ninput,lv, plotopt)
Description Given a set of regression data, xreg and yreg, the impulse response coefficient
matrix, theta, is determined via PLS. Column i of theta corresponds to the
impulse response coefficients for input i. Only a single output is allowed. The
number of inputs, ninput, and the number of latent variables, lv,mustbe
specified.
Optional output w is a matrix of dimension n (number of impulse response
coefficients) by lv consisting of orthogonal column vectors maximizing the
cross variance between input and output. Column vector cw (optional) contains
the coefficients associated with each orthogonal vector for calculating theta
(theta=w
*
cw).
Optional output ssqdif is an lv-by-2 matrix containing the percent variances
captured by PLS. The first column contains information for the input; the
second column for the output. Row i of ssqdif gives a measure of the variance
captured by using the first i latent variables.
The output residual or prediction error (yres) is also returned (optional).
No plot is produced if plotopt is equal to 0, which is the default; a plot of the
actual output and the predicted output is produced if plotopt=1;twoplots—
plot of actual and predicted output, and plot of output residual — are produced
for plotopt=2.
Example Consider the following two-input single-output system:
Load the input and output data. The input and output data were generated
from the above transfer function and random zero-mean noise was added to the
output. Sampling time of 7 minutes was used.
load plsrdat;
ys()
5.72e
14s–
60s 1+
-------------------------
1.52e
15s–
25s 1+
-------------------------
u
1
s()
u
2
s()
=
plsr
4-93
Put the input and output data in a form such that they can be used to
determine the impulse response coefficients. 30 impulse response coefficients
(n)areused.
n = 30;
[xreg,yreg] = wrtreg(x,y,n);
Determine the impulse response coefficients via plsrusing 10 latent variables.
By specifying plotopt=2, two plots — plot of predicted output and actual
output, and plot of the output residual (or predicted error) — are produced.
ninput = 2;
lv = 10;
plotopt = 2;
theta = plsr(xreg,yreg,ninput,lv,plotopt);
Use a new set of data to validate the impulse model.
[newxreg,newyreg] = wrtreg(newx,newy,n);
yres = validmod(newxreg,newyreg,theta,plotopt);
0 20 40 60 80 100 120
?3
?2
?1
0
1
2
3
Actual value (o) versus Predicted Value (+)
Sample Number
Output
0 20 40 60 80 100 120
?0.4
?0.3
?0.2
?0.1
0
0.1
0.2
0.3
Output Residual or Prediction Error
Sample Number
Residual
plsr
1-94
Convert the impulse model to a step model to be used in MPC design. Sampling
time of 7 minutes was used in determining the impulse model. Number of
outputs (1 in this case) must be specified.
nout = 1;
delt = 7;
model = imp2step(delt,nout,theta);
0 5 10 15 20 25 30
?6
?4
?2
0
2
4
Actual value (o) versus Predicted Value (+)
Sample Number
Output
0 5 10 15 20 25 30
?0.6
?0.4
?0.2
0
0.2
0.4
Output Residual or Prediction Error
Sample Number
Residual
plsr
4-95
Plot the step response coefficients.
plotstep(model)
See Also mlr, validmod, wrtreg
0 50 100 150 200 250
?1
0
1
2
3
4
5
6
u1 step response : y1
TIME
0 50 100 150 200 250
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
u2 step response : y1
TIME
poly2tfd, poly format
1-96
1poly2tfd, poly format
Purpose poly2tfd converts a transfer function (continuous or discrete) from the
standard MATLAB poly format into the MPC tf format.
Syntax g = poly2tfd(num,den)
g = poly2tfd(num,den,delt,delay)
Description Consider a continuous-time (Laplace domain) transfer function such as
or a discrete-time transfer function such as
where z is the forward-shift operator. Using the MATLAB poly format, you
would represent either of these as a numerator polynomial and a denominator
polynomial, giving the coefficients of the highest-order terms first:
num = [b
0
b
1
... b
n
];
den = [a
0
a
1
... a
n
];
If the numerator contains leading zeros, they may be omitted, i.e., the number
of elements in num can be ￡ the number of elements in den.
poly2tfd uses num and den as input to build a transfer function, g,intheMPC
tf format (see tf section for details). Optional variables you can include are:
delt
The sampling period. If this is zero or you omit it, poly2tfd assumes that you
are supplying a continuous-time transfer function. If you are supplying a
discrete-time transfer function you must specify delt.Otherwiseg will be
misinterpreted when you use it later in the MPC Toolbox functions.
delay
The time delay. For a continuous-time transfer function, delay should be in
time units. For a discrete-time transfer function, delay should be specified as
Gs()
b
0
s
n
b
1
s
n 1–
… b
n
+++
a
0
s
n
a
1
s
n 1–
… a
n
+++
----------------------------------------------------------------=
Gz()
b
0
b
1
z
1–
… b
n
z
n–
+++
a
0
a
1
z
1–
… a
n
z
n–
+++
-------------------------------------------------------------=
poly2tfd, poly format
4-97
the integer number of sampling periods of time delay. If you omit it, poly2tfd
assumes a delay of zero.
Examples Consider the continuous-time transfer function:
G(s) = 0.5 .
It has no delay. The following command creates the MPC tf format:
g=poly2tfd(0.5*[3 -1],[5 2 1]);
Now suppose there were a delay of 2.5 time units:
G(s) = 0.5 . You could use:
g=poly2tfd(0.5*[3 -1],[5 2 1],0,2.5);
Next let’s get the equivalent transfer function in discrete form. An easy way is
to get the correct poly form using cp2dp,thenusepoly2tfd to get it in the tf
form. Here are the commands to do it using a sampling period of 0.75 time
units:
delt=0.75;
[numd,dend]=cp2dp(0.5
*
[3 -1],[5 2 1],delt,rem(2.5,delt));
g=poly2tfd(numd,dend,delt,fix(2.5/delt));
Note that cp2dp is used to handle the fractional time delay and the integer
number of sampling periods of time delay is an input to poly2tfd.Theresults
are:
numd =
0 0.1232 -0.1106 -0.0607
dend =
1.0000 -1.6445 0.7408 0
g =
0 0.1232 -0.1106 -0.0607
1.0000 -1.6445 0.7408 0
0.7500 3.0000 0 0
See Also cp2dp, tf, th2mod, tfd2step
3s 1–
5s
2
2s 1++
--------------------------------
3s 1–
5s
2
2s 1++
--------------------------------e
2.5s–
scmpc
1-98
1scmpc
Purpose Simulates closed-loop systems with hard bounds on manipulated variables
and/or outputs using models in the MPC mod format. Solves the MPC
optimization problem by quadratic programming.
Syntax yp = scmpc(pmod,imod,ywt,uwt,M,P,tend,r)
[yp,u,ym] = scmpc(pmod,imod,ywt,uwt,M,P,tend, ...
r,ulim,ylim,Kest,z,d,w,wu)
Description
scmpc simulates the performance of the type of system shown in the above
diagram when there are bounds on the manipulated variables and/or outputs.
The required input variables are as follows:
pmod
Is a model in the MPC mod format that is to represent the plant.
imod
Is a model in the MPC mod format that is to be used for state estimation in the
controller. In general, it can be different from pmod if you want to simulate the
effect of plant/controller model mismatch.
Disturbances
u
d
r Plant
Setpoints
Noise-free
d
y
Plant
Outputs
S
Measured
Outputs
S
w
u
y
p
Measurement
Noise
+
+
+
+
z
Measured
Disturbances
Unmeasured
Controller
w
d
scmpc
4-99
ywt
Is a matrix of weights that will be applied to the setpoint tracking errors
(optional). If ywt=[ ], the default is equal (unity) weighting of all outputs over
theentirepredictionhorizon.Ifywt|[ ],itmusthaven
y
columns, where n
y
is
the number of outputs. All weights must be ? 0.
You may vary the weights at each step in the prediction horizon by including
up to Prows inywt. Then the first row of n
y
values applies to the tracking errors
in the first step in the prediction horizon, the next row applies to the next step,
etc. See smpccon for details on the form of the optimization objective function.
If you supply only nrow rows, where 1 ￡ nrow < P, scmpc will use the last row
to fill in any remaining steps. Thus if you wish the weighting to be the same for
all P steps, you need only specify a single row.
uwt
Is as for ywt,exceptthatuwt applies to the changes in the manipulated
variables. If you use uwt=[ ], the default is zero weighting. If uwt|[ ],itmust
have n
u
columns, where n
u
is the number of manipulated variables.
M
Therearetwowaystospecifythisvariable:
If it is a scalar, scmpc interprets it as the input horizon (number of moves) as
in DMC.
If it is a row vector containing n
b
elements, each element of the vector indicates
the number of steps over which D u = 0 during the optimization and scmpc
interprets it as a set of n
b
blocking factors. There may be 1 ￡ n
b
￡ P blocking
factors, and their sum must be ￡ P.
If you set M=[ ], the default is M=P, which is equivalent to M=ones(1,P).
P
The number of sampling periods in the prediction horizon.
tend
Is the desired duration of the simulation (in time units).
scmpc
1-100
r
Is a setpoint matrix consisting of N rows and n
y
columns, where n
y
is the
number of output variables, y:
where r
i
(k) is the setpoint for output j at time t=kT,andT is the sampling
period (as specified by the minfo vector in the mod format of pmod and imod). If
tend > NT, the setpoints vary for the first N periods in the simulation, as
specified by r, and are then held constant at the values given in the last row of
r fortheremainderofthesimulation.
In many simulations one wants the setpoints to be constant for the entire time,
in which case r need only contain a single row of n
y
values.
If you set r=[ ],thedefaultisarowofn
y
zeros.
The following input variables are optional. In general, setting one of them
equal to an empty matrix causes scmpc to use the default value, which is given
in the description.
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
scmpc
4-101
ulim
Is a matrix giving the limits on the manipulated variables. Its format is as
follows:
Note that it contains three matrices of N rows. In this case, the limits on N are
1 ￡ N ￡ n
b
,wheren
b
is the number of times the manipulated variables are to
change over the input horizon. If you supply fewer than n
b
rows, the last row
is repeated automatically.
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(2) is the lower bound for manipulated variable j for the
second move of the manipulated variables (where the first move is at the start
of the prediction horizon). If u
min
,
j
(k)=–inf, manipulated variable j will have
no lower bound for that move.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound for that move.
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
,,
,,
,,
,
()
,
()
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
……
,,
,,
,,
,
()
,
()
ulim
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
scmpc
1-102
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, cmpc will force|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The
limits on the rate of change must be nonnegative and finite.Ifyouwantittobe
unbounded, set the bound to a large number (but not too large — a value of 10
6
should work well in most cases).
The default is u
min
= –inf, u
max
=infand D u
max
=10
6
ylim
Same format as for ulim, but for the lower and upper bounds of the outputs.
The first row applies to the first point in the prediction horizon. The default is
y
min
= –inf,andy
max
=inf.
Kest
Is the estimator gain matrix. The default is the DMC estimator. See smpcest
for more details.
z
Is measurement noise that will be added to the outputs (see above diagram).
The format is the same as for r. The default is a row of n
y
zeros.
d
Is a matrix of measured disturbances (see above diagram). The format is the
same as for r, except that the number of columns is n
d
rather than n
y
.The
default is a row of n
d
zeros.s
w
Is a matrix of unmeasured disturbances (see above diagram). The format is the
same as for r, except that the number of columns is n
w
rather than n
y
.The
default is a row of n
w
zeros.
wu
Is a matrix of unmeasured disturbances that are added to the manipulated
variables (see above diagram). The format is the same as for r, except that the
number of columns is n
u
rather than n
y
.Thedefaultisarowofn
u
zeros.
Notes ? You may use a different number of rows in the matrices r, z, d, w and wu,
should that be appropriate for your simulation.
scmpc
4-103
? The ulim constraints used here are fundamentally different from the usat
constraints used in the smpcsim function. The ulim constraints are defined
relative to the beginning of the prediction horizon, which moves as the
simulation progresses. Thus at each sampling period, k,theulim constraints
apply to a block of calculated moves that begin at sampling period k and
extend for the duration of the input horizon. The usat constraints, on the
otherhand,arerelativetothefixedpointt = 0, the start of the simulation.
The calculated outputs are as follows (all but yp are optional):
yp
Is a matrix containing M rows and n
y
columns, where M = max(fix(tend/T)+1,
2). The first row will contain the initial condition, and row k – 1 will give the
values of the noise-free plant outputs, (see above diagram), at time t=kT.
u
Is a matrix containing the same number of rows asypand n
u
columns. The time
corresponding to each row is the same as for yp. The elements in each row are
the values of the manipulated variables, u (see above diagram).
Note Theuvaluesarethosecomingfromthecontrollerbeforetheadditionof
the unmeasured disturbance, w
u
.
ym
Is a matrix of the same structure as yp, containing the values of the predicted
output from the state estimator in the controller. These will, in general, differ
from those in yp if imod|pmod and/or there are unmeasured disturbances. The
prediction includes the effect of the most recent measurement, i.e., it is .
For unconstrained problems, scmpc and smpcsim should give the same results.
The latter will be faster because it uses an analytical solution of the QP
problem, whereas scmpc solves it by iteration.
y
p
y? kk()
scmpc
1-104
Examples Consider the linear system:
The following statements build the model and set up the controller in the same
way as in the smpcsim example.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2;
imod=tfd2mod(delt,ny,g11,g21,g12,g22);
pmod=imod; P=6; M=2; ywt=[ ]; uwt=[1 1];
tend=30; r=[0 1];
Here, however, we will demonstrate the effect of constraints. First we set a
limitof0.1ontherateofchangeofu
1
and a minimum of –0.15 for u
2
.
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[ ];
[y,u]=scmpc(pmod,imod,ywt,uwt,M,P,tend,r,ulim,ylim);
plotall(y,u,delt), pause
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
scmpc
4-105
Note that D u
2
has a large (but finite) limit. It never comes into play.
We next apply a lower bound of zero to both outputs:
ulim=[-inf -0.15 inf inf 0.1 100];
ylim=[0 0 inf inf];
[y,u]=scmpc(pmod,imod,ywt,uwt,M,P,tend,r,ulim,ylim);
plotall(y,u,delt),
pause
The following results show that no constraints are violated.
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.2
?0.18
?0.16
?0.14
?0.12
?0.1
?0.08
?0.06
Manipulated Variables
Time
scmpc
1-106
Restrictions ? Initial conditions of zero are used for all states in imod and pmod.This
simulates the condition where all variables represent a deviation from a
steady-state initial condition.
? The first n
u
+ n
d
columns of the D matrices in imod and pmod must be zero.
In other words, neither u nor d may have an immediate effect on the outputs.
Suggestions Problems with many inequality constraints can be very time consuming. You
can minimize the number of constraints by:
? Using small values for P and/or M.
? Leaving variables unconstrained (limits at ±inf) intermittently unless you
think the constraint is important.
See Also plotall, ploteach, smpccl, smpccon, smpcest, smpcsim
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.2
?0.1
0
Manipulated Variables
Time
sermod
4-107
1sermod
Purpose Puts two models in series by connecting the output of one to the input of the
other. Mimics the series function of the Control System Toolbox, except that
sermod worksonmodelsintheMPCmod format.
Syntax pmod = sermod(mod1,mod2)
Description
mod1 and mod2 are models in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally
create them using either the tfd2mod, ss2mod or th2mod functions.
sermodcombines them to form a composite system, pmod,asshownintheabove
diagram. It is also in the mod format. Note how the inputs to mod1 and mod2
are ordered in pmod.
Restrictions ?mod1 and mod2 must have been created with equal sampling periods.
? The number of measured outputs in mod1 must equal the number of
manipulated variables in mod2.
See Also addmd, addmod, addumd, appmod, paramod
mod2
pmod
y
1
=u
2
u
1
y
mod1
u
2
d
1
d
2
w
1
w
2
u
1
d
1
d
2
w
1
w
2
y
smpccl
1-108
1smpccl
Purpose Combines a plant model and a controller model in the MPC mod format,
yielding a closed-loop system model in the MPC format. This can be used for
stability analysis and linear simulations of closed-loop performance.
Syntax [clmod,cmod] = smpccl(pmod,imod,Ks)
[clmod,cmod] = smpccl(pmod,imod,Ks,Kest)
Description
pmod
Is a model (in the mod format) representing the plant in the above diagram.
imod
Is a model (in the same format) that is to be used to design the MPC controller
blockshowninthediagram.Itmaybethesameaspmod (in which case there is
no model error in the controller design), or it may be different.
Ks
Is a controller gain matrix, which must have been calculated by the function
smpccon.
Disturbances
u
d
r Plant
Setpoints
Noise-free
d
y
Plant
Outputs
S
Measured
Outputs
S
w
u
y
p
Measurement
Noise
+
+
+
+
z
Measured
Disturbances
Unmeasured
Controller
w
d
smpccl
4-109
Kest
Is an (optional) estimator gain matrix. If omitted or set to an empty matrix, the
default is to use the DMC estimator index DMC estimator. See the
documentation for the function smpcest for more details on the design and
proper format of Kest.
smpccl
Calculates a model of the closed-loop system, clmod.Itisinthemod format and
can be used, for example, with analysis functions such as smpcgain and
smpcpole, and with simulation routines such as mod2step and dlsimm. smpccl
also calculates a model of the controller element, cmod.
The closed-loop model, clmod, has the following state-space representation:
x
cl
(k +1)=F
cl
x
cl
(k)+G
cl
u
cl
(k)
y
cl
(k)=C
cl
x
cl
(k)+D
cl
u
cl
(k)
where x
cl
is a vector of n state variables, u
cl
is a vector of input variables, y
cl
is
a vector of outputs, and F
cl
, G
cl
, C
cl
,andD
cl
are matrices of appropriate size.
Theexpertusermaywanttoknowthesignificanceofthestatevariablesinx
cl
.
They are (in the following order):
? The n
p
states of the plant (as specified in pmod),
? The n
i
changes in the state estimates (based on the model specified in imod
and the estimator gain, Kest),
? The n
y
estimates of the noise-free plant output (from the state
estimator),
? n
u
integrators that operate on the D u signal produced by the standard MPC
formulation to yield a u signal that can be used as input to the plant and as
a closed-loop output, and
? n
d
differencing elements that operate on the d signal to produce the D d signal
required in the standard MPC formulation. If there are no measured
disturbances, these states are omitted.
y? kk 1–()=
smpccl
1-110
The closed-loop input and output variables are:
where istheestimateofthenoise-freeplantoutputatsampling
period k based on information available at period k.Thisestimateisgenerated
by the controller element.
Note that u
cl
will include d and/or w automatically whenever pmod includes
measured disturbances and/or unmeasured disturbances. Thus the length of
the u
cl
vector will depend on the inputs you have defined in pmod and imod.
Similarly, ycl depends on the number of outputs and manipulated variables.
Let m and p be the lengths of u
cl
and y
cl
,respectively.Then
m=2n
y
+n
u
+n
d
+n
w
p=2n
y
+n
u
The state-space form of the controller model, cmod,canbewrittenas:
x
c
(k +1)=F
c
x
c
(k)+G
cl
u
c
(k)
y
c
(k)=C
c
x
c
(k)+D
c
u
c
(k)
where
and the controller states are the same as those of the closed loop system except
that the n
p
plantstatesarenotincluded.
u
cl
k()
rk()
zk()
w
u
k()
dk()
wk()
and y
cl
k()
y
p
k()
uk()
y? kk()
==
y? kk()=
u
c
k()
rk()
yk()
dk()
and y
c
k() uk()==
smpccl
4-111
Examples Consider the linear system:
We build this model using the MPC Toolbox functions poly2tfd and tfd2mod.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2;
imod=tfd2mod(delt,ny,g11,g21,g12,g22);
pmod=imod; % No plant/model mismatch
Now we design the controller. Since there is delay, we use M < P:Wespecify
the defaults for the other tuning parameters, uwt and ywt, then calculate the
controller gain:
P=6; % Prediction horizon.
M=2; % Number of moves (input horizon).
ywt=[ ]; % Output weights (default - unity on
% all outputs).
uwt=[ ]; % Man. Var weights (default - zero on
% all man. vars).
Ks=smpccon(imod,ywt,uwt,M,P);
Now we can calculate the model of the closed-loop system and check its poles
for stability:
clmod=smpccl(pmod,imod,Ks);
maxpole=max(abs(smpcpole(clmod)))
The result is:
maxpole = 0.8869
Since this is less than 1, the plant and controller combination will be
closed-loop stable. (The closed-loop system has 20 states in this example).
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
smpccl
1-112
You can also use the closed-loop model to calculate and plot the step response
with respect to all theinputs.Theappropriatecommandsare:
tend=30;
clstep=mod2step(clmod,tend);
plotstep(clstep)
Since the closed-loop system has m =6inputsandp = 6 outputs, only one of the
plots is reproduced here. It shows the response of the first 4 closed-loop outputs
to a step in the first closed-loop input, which is the setpoint for y
1
:
Closed-loop outputs y
1
and y
2
are the true plant outputs (noise-free). Output y
1
goes to the new setpoint quickly with a small overshoot. This causes a small,
short-term disturbance in y
2
. The plots for y
3
and y
4
show the required
variation in the manipulated variables.
0 10 20 30
0
0.2
0.4
0.6
0.8
1
1.2
1.4
u1 step response : y1
TIME
0 10 20 30
?0.2
?0.15
?0.1
?0.05
0
u1 step response : y2
TIME
0 10 20 30
0
0.05
0.1
0.15
0.2
0.25
u1 step response : y3
TIME
0 10 20 30
0
0.02
0.04
0.06
0.08
0.1
0.12
u1 step response : y4
TIME
smpccl
4-113
The following commands show how you could use dlsimm to calculate the
response of the closed-loop system to a step in the setpoint for y
1
,withadded
random measurement noise.
r=[ones(11,1) zeros(11,1)];
z=0.1*rand(11,2);
wu=zeros(11,2);
d=[ ];
w=[ ];
ucl=[r z wu d w];
[phicl,gamcl,ccl,dcl]=mod2ss(clmods);
ycl=dlsimm(phicl,gamcl,ccl,dcl,ucl);
y=ycl(:,1:2); u=ycl(:,3:4); ym=ycl(:,5:6);
Restrictions ?imod and pmod must have been created using the same sampling period, and
an equal number of outputs, measured disturbances, and manipulated
variables.
? Both imod and pmod must be strictly proper, i.e., the D matrices in their
state-space descriptions must be zero. Exception: the last n
w
columns of the
D matrices may be nonzero, i.e., the unmeasured disturbance may have an
immediate effect on the outputs.
See Also mod2step, scmpc, smpccon, smpcest, smpcgain, smpcpole, smpcsim
smpccon
1-114
1smpccon
Purpose Calculates MPC controller gain using a model in MPC mod format.
Syntax Ks = smpccon(imod)
Ks = smpccon(imod,ywt,uwt,M,P)
Description Combines the following variables (most of which are optional and have default
values) to calculate the state-space MPC gain matrix, Ks.
imod is the model of the process to be used in the controller design (in the
mod format).
The following input variables are optional:
ywt
Is a matrix of weights that will be applied to the setpoint tracking errors. If you
use ywt=[ ] or omit it, the default is equal (unity) weighting of all outputs over
theentirepredictionhorizon.Ifywt|[ ],itmusthaven
y
columns, where n
y
is
the number of outputs. All weights must be ? 0.
You may vary the weights at each step in the prediction horizon by including
up to Prows inywt. Then the first row of n
y
values applies to the tracking errors
in the first step in the prediction horizon, the next row applies to the next step,
etc.
If you supply only nrow rows, where 1 ￡ nrow < P, smpccon will use the last row
to fill in any remaining steps. Thus if you wish the weighting to be the same for
all P steps, you need only specify a single row.
uwt
Same format as ywt,exceptthatuwt applies to the changes in the
manipulated variables. If you use uwt=[ ] or omit it, the default is zero
weighting. If uwt|[ ],itmusthaven
u
columns, where n
u
is the number of
manipulated variables.
M
Therearetwowaystospecifythisvariable:
If it is a scalar, smpccon interprets it as the input horizon (number of moves) as
in DMC.
smpccon
4-115
If it is a row vector containing n
b
elements, each element of the vector indicates
the number of steps over which D u = 0 during the optimization and smpccon
interprets it as a set of n
b
blocking factors. There may be 1 ￡ n
b
￡ P blocking
factors, and their sum must be ￡ P.
If you set M=[ ] or omit it, the default is M=P, which is equivalent to
M=ones(1,P).
P
The number of sampling periods in the prediction horizon. If you set P=[ ] or
omit it, the default is P=1.
If you take the default values for all the optional variables, you get the “perfect
controller,” i.e., a model-inverse controller. This controller is not applicable in
the following situations:
? When one or more outputs cannot respond to the manipulated variables
within 1 sampling period due to time delay, the plant-inverse controller is
unrealizable. To counteract this you can penalize changes in the
manipulated variables (variable uwt), use blocking (variable M), and/or make
P>>M.
? When imod contains transmission zeros outside the unit circle the
plant-inverse controller will be unstable. To counteract this, you can use
blocking (variable M), restrict the input horizon (variable M), and/or penalize
changes in the manipulated variables (variable uwt).
The model-inverse controller is also relatively sensitive to model error and is
best used as a point of reference from which you can progress to a more robust
design.
Algorithm The controller gain is a component of the solution to the optimization problem:
Minimize
with respect to (a series of current and future moves in the manipulated
variables), where (k + j) is a prediction of output i at a time j sampling
periods into the future (relative to the current time, k), which is a function of
Jk() ywt
i
j()r
i
kj+()– y?
i
kj+()[]
2
i 1=
n
y
Ge5
j 1=
p
Ge5
=
uwt
i
j()D u?
i
j()()
2
i 1=
n
u
Ge5
j 1=
n
b
Ge5
+
D u?
i
j()
y?
i
smpccon
1-116
(j), r
i
(k + j) is the corresponding future setpoint, and n
b
is the number of
blocks or moves of the manipulated variables.
References Ricker, N. L. “Use of Quadratic Programming for Constrained Internal Model
Control,” Ind. Eng. Chem. Process Des. Dev., 1985, 24, 925–936.
Ricker, N. L. “Model-predictive control with state estimation,” I & EC Res.,
1990, 29, 374.
Example Consider the linear system:
See the smpccl example for the commands that build the model and a simple
controller for this process.
Here is a slightly more complex design with blocking and time-varying weights
on the manipulated and output variables:
P=6; M=[2 4];
uwt=[1 0; 0 1];
ywt=[1 0.1; 0.8 0.1; 0.1 0.1];
Ks=smpccon(imod,ywt,uwt,M,P);
tend=30; r=[1 0];
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
There is no particular rationale for using time varying weights in this case —
it is only for illustration. The manipulated variables will make 2 moves during
the prediction horizon (see value of M,above).Theuwt selection gives u
1
a unity
weight and u
2
a zero weight for the first move, then switches the weights for
the second move. If there had been any additional moves they would have had
the same weighting as the second move.
D u?
i
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
smpccon
4-117
The ywt value assigns a constant weight of 0.1 to y
2
, and a weight that
decreases over the first 3 periods to y
1
. The weights for periods 4 to 6 are the
same as for period 3. The resulting closed-loop (servo) response is:
See Also scmpc, smpccl, smpcsim
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
Outputs
Time
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
Manipulated Variables
Time
smpcest
1-118
1smpcest
Purpose Sets up a state-estimator gain matrix for use with MPC controller design and
simulation routines using models in MPC mod format. Can use either a
disturbance/noise model that you specify, or a simplified form in which each
output is affected by an independent disturbance (plus measurement noise).
Syntax For the general case:
[Kest] = smpcest(imod,Q,R)
For simplified disturbance modeling:
[Kest,newmod] = smpcest(imod)
[Kest,newmod] = smpcest(imod,tau,signoise)
Description
In the above block diagram, u is a vector of n
u
manipulated variables (n
u
? 1),
d is a vector of n
d
measured disturbances (n
d
? 0), w is a vector of unmeasured
disturbances, z is measurement noise, y is a vector of outputs, and represents
these outputs before the addition of measurement noise. The objective of the
state estimator in MPC is to estimate the present and future values of ,
rejecting as much of the measurement noise as possible. The inputs u and d are
assumed perfectly measurable, whereas w and z are unknown and must be
inferred from the measurements. G
w
is a transfer function matrix representing
the effect of each element of w on each output in y.
General Case imod
Is the model (in mod format) to be used as the basis for the state estimator. It
should be the same as that used to calculate the controller gain (see smpccon).
It must include a model of the disturbances, i.e., the G
w
element in the above
yS
u
+
G
w
Plant
z
y
w
S
+
+
+
d
y
y
smpcest
4-119
diagram. You could, for example, use addumd to combine a plant and
disturbance model, yielding a composite model in the proper form.
Q
Is a symmetric, positive semi-definite matrix giving the covariances of the
disturbances in w.Itmustben
w
by n
w
,wheren
w
(? 1) is the number of
unmeasured disturbances in imod (i.e., the length of w).
R
Is a symmetric, positive-definite matrix giving the covariances of the
measurement noise, z.Itmustben
ym
by n
ym
,wheren
ym
(? 1) is the number of
measured outputs in imod.
The calculated output variable is:
Kest
The estimator gain matrix. It will contain n + n
y
rows and n
ym
columns, where
n is the number of states in imod,andn
y
is the total number of outputs
(measured plus unmeasured).
Simplified
disturbance
modeling
For the simplified disturbance/noise model we make the following
assumptions:
? The vectors w, z, y and are all length n
y
.
? G
w
is diagonal. Thus each element of w affects one (and only one) element of
y. Diagonal element G
wi
has the discrete (sampled-data) form:
where a
i
= e
–T/
t
i
,0￡t
i
￡￥ ,andT is the sampling period.
As , Gwi(q) approaches a unity gain, while as , G
wi
becomes an
integrator.
? Element i of D w is a stationary white-noise signal with zero mean and
standard deviation s
wi
(where w
i
(k)=w
i
(k)–w
i
(k –1)).
? Element i of z is a stationary white-noise signal with zero mean and standard
deviation s
zi
.
y
G
wi
q()
1
qa
i
–
--------------=
t
i
0fit
i
￥fi
smpcest
1-120
The input variables are then as follows:
imod
Is the model (in mod format) to be used as the basis for the state estimator. It
should be the same as that used to calculate the controller gain (see smpccon).
tau
Is a row vector,lengthn
y
, giving the values of t
i
to be used in eq. 1. Each element
must satisfy: 0 ￡t
i
￡￥ .Ifyouusetau=[ ], smpcest uses the default, which is
n
y
zeros.
signoise
Is a row vector,lengthn
y
, giving the signal-to-noise ratio for the each
disturbance, defined as g
i
= s
wi
=s
zi
. Each element must be nonnegative. If
omitted, smpcsim uses an infinite signal-to-noise ratio for each output.
The calculated output variables are:
Kest
The estimator gain matrix.
newmod
The modified version of imod, which must be used in place of imod in any
simulation/analysis functions that require Kest (e.g., smpccl, smpcsim, scmpc).
If imod contains n states, and there are n
1
outputs for which t
i
>0,then
newmod will have n + n
1
states. The optimal gain matrix, Kest, will have n +
n
1
+ n
y
rows and n
ym
columns. The first n rows will be zero, the next n
1
rows
will have the gains for the estimates of the n
1
added states (if any), and the last
n
y
rows will have the gains for estimating the noise-free outputs, .
Examples Consider the linear system:
y
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
-----------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
-----------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
3.8e
8– s
14.9s 1+
------------------------
4.9e
3s–
13.2s 1+
------------------------
ws()+=
smpcest
4-121
The following statements build two models: pmod, which contains the model of
the disturbance, w,andimod, which does not.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=1; ny=2;
imod=tfd2mod(delt,ny,g11,g21,g12,g22);
gw1=poly2tfd(3.8,[14.9 1],0,8);
gw2=poly2tfd(4.9,[13.2 1],0,3);
pmod=addumd(imod,tfd2mod(delt,ny,gw1,gw2));
Calculate the gain for a typical MPC (unconstrained) controller
P=6; M=2;
ywt=[ ]; uwt=[1 1];
Ks=smpccon(imod,ywt,uwt,M,P);
Next design an estimator using the G
w
model in pmod.ThechoicesofQandR
are arbitrary. R was made relatively small (since measurement noise will be
negligible in the simulations).
Kest1=smpcest(pmod,1,0.001*eye(ny));
Ks1=smpccon(pmod,ywt,uwt,M,P);
Now design another estimator using a simplified disturbance model in which
each output is affected by a disturbance with a first-order time constant of 10
and a signal-to-noise ratio of 3.
tau=[10 10]; signoise=[3 3];
[Kest2,newmod]=smpcest(imod,tau,signoise);
Ks2=smpccon(newmod,ywt,uwt,M,P);
Compare the performance of these two estimators to the default (DMC)
estimator when there is a unit step in w:
r=[ ]; ulim=[ ]; z=[ ]; d=[ ]; w=[1]; wu=[ ]; tend=30;
[y
1
,u
1
]=smpcsim(pmod,pmod,Ks1,tend,r,ulim,Kest1, z,d,w,wu);
[y
2
,u
2
]=smpcsim(pmod,newmod,Ks2,tend,r,ulim,Kest2, z,d,w,wu);
[y
3
,u3]=smpcsim(pmod,imod,Ks,tend,r,ulim,[ ],z,d,w,wu);
smpcest
1-122
The solid lines in the following plots are for y
1
(or u
1
) and the dashed lines are
for y
2
(or u
2
). Both outputs have setpoints at zero. You can see that the default
estimator is much more sluggish than the others in counteracting this type of
disturbance. The simplified disturbance design does nearly as well as that
using the exact model of the disturbances. The main difference is that it allows
more error in y
1
following the disturbance in y
2
.
The first 14 states in both imod and pmod arefortheresponseoftheoutputsto
u. Since the unmeasured disturbance has no effect on them, their gains are
2
1
0
-1
0102030
Time
Outputs using general estimator
2
1
0
-1
0 102030
Time
Outputs using simplified estimator
2
1
0
-1
0102030
Time
Outputs using default estimator
1
0.5
0
0102030
Time
Man.vars. using general estimator
smpcest
4-123
zero. pmod contains 10 additional disturbance states and there are 2 outputs,
so the last 12 rows of Kest1 are nonzero:
Kest1(15:26,:)=
-0.0556 8.8659
-0.0594 7.1499
-0.0635 5.1314
-0.0679 2.7748
-0.0725 0.0411
-0.0781 -0.0182
-0.0915 -0.0008
-0.0520 0.0001
1.2663 0.0000
0.0281 -0.0000
0.3137 0.0000
0.0000 0.9925
and the last 4 rows of Kest2 are nonzero:
Kest2(15:18,:) =
0.7274 0
0 0.7274
0.9261 0
0 0.9261
Algorithm In the general case, smpcest uses dlqe2 to calculate the optimal estimator
gain, Kest. In the simplified case, it uses an analytical solution of the discrete
Riccati equation (which is possible to obtain in this case because the
disturbances are independent with low-order dynamics).
The number of rows in Kest is larger than that in newmod because the MPC
analysis and simulation functions augment the model states with the outputs
(see mpcaugss), and Kest must be set up to account for this.
If all t
i
=0andallg
i
= ￥ ,wegettheDMCestimator, which has n rows of zeros
followed by an identity matrix of dimension n
y
. This is the default for all of the
MPC analysis and simulation routines that require an estimator gain as input.
smpcest
1-124
Important note: smpcest decides whether you are using the general case or
the simplified approach by checking the number of output arguments you
have supplied. If there is only one, it assumes you want the general case.
Otherwise, it proceeds as for the simplified case. It checks the dimensions of
your input arguments to make sure they are consistent with this decision.
If you get unexpected results or an error message, make sure you have specified
the correct number of output arguments.
See Also scmpc, smpccl, smpccon, smpcsim
smpcgain, smpcpole
4-125
1smpcgain, smpcpole
Purpose Calculates steady-state gain matrix or poles for a system in the MPC mod
format.
Syntax g = smpcgain(mod)
poles = smpcpole(mod)
Description mod is a dynamic model in the MPC mod format. smpcgain and smpcpole
convert it to its equivalent state-space form:
x(k +1)=F x(k)+G v (k)
y(k)=Cx(k)+Dv(k)
where n includes all of the inputs in mod. smpcgain then calculates the gain
matrix:
G = C(I – F )
1
G + D
which contains n
y
rows, corresponding to each of the outputs in mod,andn
u
+
n
d
+ n
w
columns, corresponding to each of the inputs.
smpcpole calculates the poles, i.e., the eigenvalues of the F matrix.
Example See smpccl for an example of the use of smpcpole.
Restriction If mod is not asymptotically stable, smpcgain terminates with an error message.
See Also mod
smpcsim
1-126
1smpcsim
Purpose Simulates closed-loop systems with saturation constraints on the manipulated
variables using models in the MPC mod format. Can also be used for open-loop
simulations.
Syntax yp = smpcsim(pmod,imod,Ks,tend,r)
[yp,u,ym] = smpcsim(pmod,imod,Ks,tend,r,usat,...
Kest,z,d,w,wu)
Description
smpcsim provides a convenient way to simulate the performance of the type of
system shown in the above diagram. The required input variables are as
follows:
pmod
Is a model in the MPC mod format that is to represent the plant.
imod
Is a model in the MPC mod format that is to be used for state estimation in the
controller. In general, it can be different from pmod if you wish to simulate the
effect of plant/controller model mismatch. Note, however, that imod should be
thesameasthatusedtocalculateKs.
Disturbances
u
d
r Plant
Setpoints
Noise-free
d
y
Plant
Outputs
S
Measured
Outputs
S
w
u
y
p
Measurement
Noise
+
+
+
+
z
Measured
Disturbances
Unmeasured
Controller
w
d
smpcsim
4-127
Ks
Is the MPC controller gain matrix, usually calculated using the function
smpccon.
If you set Ksto an empty matrix, smpcsimwill do an open-loop simulation. Then
the inputs to the plant will be r (which must be set to the vector of manipulated
variablesinthiscase),d, w,andwu. The measurement noise input, z,willbe
ignored.
tend
Is the desired duration of the simulation (in time units).
r
Is normally asetpointmatrixconsistingofN rows and n
y
columns, where n
y
is
the number of output variables, y:
where r
i
(k) is the setpoint for output j at time t=kT,andT is the sampling
period (as specified by the minfo vector in the mod format of pmod and imod). If
tend > NT, the setpoints vary for the first N periods in the simulation, as
specified by r, and are then held constant at the values given in the last row of
r fortheremainderofthesimulation.
In many simulations one wants the setpoints to be constant for the entire time,
in which case r need only contain a single row of n
y
values.
If you set r=[ ], the default is a row of n
y
zeros.
For open-loop simulations, r specifies the manipulated variables and must
contain n
u
columns.
The following input variables are optional. In general, setting one of them
equal to an empty matrix causes smpcsim to use the default value, which is
given in the description.
r
r
1
1() r
2
1() … r
n
y
1()
r
1
2() r
2
2() … r
n
y
2()
…
r
1
N()r
2
N() … r
n
y
N()
=
… … …
smpcsim
1-128
usat
Is a matrix giving the saturation limits on the manipulated variables. Its
format is as follows:
Note that it contains three matrices of N rows. N may be different than that for
the setpoint matrix, r, but the idea is the same: the saturation limits will vary
for the first N sampling periods of the simulation, then be held constant at the
values given in the last row of usat for the remaining periods (if any).
The first matrix specifies the lower bounds on the n
u
manipulated variables.
For example, u
min
,
j
(k) is the lower bound for manipulated variable j at time t=
kT in the simulation. If u
min
,
j
(k)=–inf, manipulated variable j will have no
lower bound at t=kT.
The second matrix gives the upper bounds on the manipulated variables. If
u
max,j
(k) =inf, manipulated variable j will have no upper bound at t=kT.
The lower and upper bounds may be either positive or negative (or zero) as long
as u
min
,
j
(k) ￡ u
max,j
(k).
,,
,,
,,
,
()
,
()
,,
,,
,,
,
()
,
()
usat
u
min 1,
1() … u
min n
u
,
1()
u
min 1,
2() … u
min n
u
,
2()
…
u
min 1,
N() … u
min n
u
,
N()
=
D u
max 1,
1() … Du
max n
u
,
1()
D u
max 1,
2() … Du
max n
u
,
2()
…
D u
max 1,
N() … Du
max n
u
,
N()
u
max 1,
1() … u
max n
u
,
1()
u
max 1,
2() … u
max n
u
,
2()
…
u
max 1,
N() … u
max n
u
,
N()
… …
……
……
smpcsim
4-129
The third matrix gives the limits on the rate of change of the manipulated
variables. In other words, smpcsim will force|u
j
(k)–u
j
(k –1)|￡D u
max,j
(k). The
limits on the rate of change must be nonnegative.
The default is no saturation constraints, i.e., all the umin values will be set to
–inf, and all the u
max
and D u
max
values will be set to inf.
Note: Saturation constraints are enforced by simply “clipping” the
manipulated variable moves so that they satisfy all constraints. This is a
nonoptimal solution that, in general, will differ from the results you would get
using the ulim variable in scmpc.
Kest
Is the estimator gain matrix. The default is the DMC estimator. See smpcest
for more details.
z
Is measurement noise that will be added to the outputs (see above diagram).
The format is the same as for r.Thedefaultisarowofn
y
zeros.
d
Is a matrix of measured disturbances (see above diagram). The format is the
same as for r, except that the number of columns is n
d
rather than n
y
.The
default is a row of n
d
zeros.
w
Is a matrix of unmeasured disturbances (see above diagram). The format is the
same as for r, except that the number of columns is n
w
rather than n
y
.The
default is a row of n
w
zeros.I
wu
Is a matrix of unmeasured disturbances that are added to the manipulated
variables (see above diagram). The format is the same as for r,exceptthatthe
number of columns is n
u
rather than n
y
.Thedefaultisarowofn
u
zeros.
smpcsim
1-130
Note: You may use a different number of rows in the matrices r, usat, z, d, w
and wu, should that be appropriate for your simulation.
The calculated outputs are as follows (all but yp are optional):
yp
Is a matrix containing M rows and n
y
columns, where
M = max(fix(tend=T) + 1, 2). The first row will contain the initial condition, and
row k – 1 will give the values of the plant outputs, y (see above diagram), at
time t=kT.
u
Is a matrix containing the same number of rows asypand n
u
columns. The time
corresponding to each row is the same as for yp. The elements in each row are
the values of the manipulated variables, u (see above diagram).
Note: The u values are those coming from the controller before the addition
of the unmeasured disturbance, w
u
.
ym
Is a matrix of the same structure as yp, containing the values of the predicted
output from the state estimator in the controller. These will, in general, differ
from those in yp if imod|pmod and/or there are unmeasured disturbances. The
prediction includes the effect of the most recent measurement, i.e., it is .
Examples Consider the linear system:
y? kk()
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
smpcsim
4-131
The following statements build the model and calculate the MPC controller
gain:
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2;
imod=tfd2mod(delt,ny,g11,g21,g12,g22);
pmod=imod;
P=6; M=2;
ywt=[ ]; uwt=[1 1];
Ks=smpccon(imod,ywt,uwt,M,P);
Simulate and plot the closed-loop performance for a unit step in the setpoint for
y
2
, occurring at t =0.
tend=30; r=[0 1];
[y,u]=smpcsim(pmod,imod,Ks,tend,r);
plotall(y,u,delt),pause
0 5 10 15 20 25 30
?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Outputs
Time
0 5 10 15 20 25 30
?0.24
?0.22
?0.2
?0.18
?0.16
?0.14
?0.12
?0.1
Manipulated Variables
Time
smpcsim
1-132
Try a pulse change in the disturbance that adds to u
1
:
r=[ ]; usat=[ ]; Kest=[ ]; z=[ ]; d=[ ]; w=[ ];
wu=[ 1 0; 0 0];
[y,u]=smpcsim(pmod,imod,Ks,tend,r,usat,Kest,z,d,w,wu);
plotall(y,u,delt),pause
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
1.5
Outputs
Time
0 5 10 15 20 25 30
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
0.1
Manipulated Variables
Time
smpcsim
4-133
For the same disturbance as in the previous case, limit the rates of change of
both manipulated variables.
usat=[-inf -inf inf inf 0.1 0.05];
[y,u]=smpcsim(pmod,imod,Ks,tend,r,usat,Kest,z,d,w,wu);
plotall(y,u,delt),pause
Restrictions ? Initial conditions of zero are used for all states in imod and pmod.This
simulates the condition where all variables represent a deviation from a
steady-state initial condition.
? The first n
u
+ n
d
columns of the D matrices in pmod and imod must be zero.
In other words, neither u nor d may have an immediate effect on the outputs.
See Also plotall, ploteach, scmpc, smpccl, smpccon, smpcest
0 5 10 15 20 25 30
?1
?0.5
0
0.5
1
1.5
2
Outputs
Time
0 5 10 15 20 25 30
?0.3
?0.2
?0.1
0
0.1
Manipulated Variables
Time
ss2mod
1-134
1ss2mod
Purpose Converts a discrete-time state-space system model into the MPC mod format.
Syntax pmod = ss2mod(phi,gam,c,d)
pmod = ss2mod(phi,gam,c,d,minfo)
Description
Consider the process shown in the above block diagram. ss2mod assumes the
following state-space representation:
where x is a vector of n state variables, u represents n
u
manipulated variables,
d represents n
d
measured disturbances, w represents n
w
unmeasured
disturbances, y is a vector of n
y
plant outputs, z is measurement noise, and F ,
G
u
, etc., are constant matrices of appropriate size. The variable (k)represents
the plant output before the addition of measurement noise. We further define:
D =[D
u
D
d
D
w
]
G =[G
u
G
d
G
w
]
ss2mod uses the F, G , C, and D matrices you supply to build a model, pmod,in
the MPC mod format. See the mod section for more details.
S
y y
Unmeasured
Disturbances
Manipulated
Variables
Measured
Disturbances
Measurement
Noise
w
u
d
++
z
Plant
Measured
Outputs
xk 1+()Fxk() G
u
uk() G
d
dk() G
w
wk()+++=
yk() yk() zk()+=
Cx k() D
u
uk() D
d
dk() D
w
wk() zk()+++ +=
y
ss2mod
4-135
You can also divide the outputs into n
ym
measured outputs and n
yu
unmeasured outputs, where n
ym
+ n
yu
= n
y
. Then the first n
ym
elements in y
and the first n
ym
rows in C and D are assumed to be for the measured outputs,
and the rest are for the unmeasured outputs.
minfo is an optional variable that allows you to specify certain characteristics
of the system. The general form is a row vector with 7 elements, the
interpretation of which is as follows:
If you specify minfo as a scalar, ss2mod takesitasthesamplingperiodandsets
the remaining elements of minfo as follows:
minfo(2) = # rows in phi, minfo(3) = # columns in gam,
minfo(4) = minfo(5) = 0, minfo(6) = # rows in c, minfo(7) = 0.
In other words, the default is to assume that all inputs are manipulated
variables and all outputs are measured. If you omit minfo, ss2mod sets the
sampling period to 1 and uses the defaults for the remaining elements.
minfo (1) T, the sampling period used to create the model.
(2) n, the number of states.
(3) n
u
, the number of manipulated variable inputs.
(4) n
d
, the number of measured disturbances.
(5) n
w
, the number of unmeasured disturbances.
(6) n
ym
, the number of measured outputs.
(7) n
yu
, the number of unmeasured outputs.
ss2mod
1-136
Example
Suppose you have the situation shown in the above diagram where u, d, w,and
y are scalar signals, and the three transfer functions are first-order:
The sampling period is T =2.
One way to build the model of the complete system is to convert these to
state-space form and use ss2mod:
[phiu,gamu,cu,du]=tf2ss(0.7,[1 -0.9]);
[phid,gamd,cd,dd]=tf2ss(-1.5,[1 -0.85]);
[phiw,gamw,cw,dw]=tf2ss(1,[1 0.6]);
[phi,gam,c,d]=mpcparal(phiu,gamu,cu,du,phid,gamd,cd,dd);
[phi,gam,c,d]=mpcparal(phi,gam,c,d,phiw,gamw,cw,dw);
delt=2;
minfo=[delt 3 1 1 1 1 0];
pmod=ss2mod(phi,gam,c,d,minfo)
You must be careful to build up the parallel structure in the correct order. For
example, the columns corresponding to G
u
must always come first in G .
yS
u
+
G
d
G
u
y
d
S
+
+
+
G
w
w
G
d
z()
1.5–
10.85z
1–
–
-----------------------------=G
u
z()
0.7
10.9z
1–
–
--------------------------=
G
w
z()
1
10.6z
1–
+
--------------------------=
ss2mod
4-137
Another, more foolproof way is to use the addmd and addumd functions:
ny=1;
gu=poly2tfd(0.7,[1 -0.9],delt);
gd=poly2tfd(-1.5,[1 -0.85],delt);
gw=poly2tfd(1,[1 0.6],delt);
pmod=tfd2mod(delt,ny,gu);
pmod=addmd(pmod,tfd2mod(delt,ny,gd));
pmod=addumd(pmod,tfd2mod(delt,ny,gw))
Using either approach, the result is:
pmod =
2.0000 3.0000 1.0000 1.0000 1.0000 1.0000 0
NaN 0.9000 0 0 1.0000 0 0
0 0 0.8500 0 0 1.0000 0
0 0 0 -0.6000 0 0 1.0000
0 0.7000 -1.5000 1.0000 0 0 0
See Also mod format, mod2ss
ss2step
1-138
1ss2step
Purpose Uses a model in state-space format to calculate the step response of a SISO or
MIMO system, in MPC step format.
Syntax plant = ss2step(phi,gam,c,d,tfinal)
plant = ss2step(phi,gam,c,d,tfinal,delt1,delt2,nout)
Description The input variables phi, gam, c,andd are assumed to be a state-space model
of a process. The model can be either continuous time:
or discrete time:
x(k +1)=F x(k)+G u(k)
y(k)=Cx(k)+Du(k)
where x is a vector of n state variables, u is a vector of n
u
inputs (usually but
not necessarily manipulated variables), y is a vector of n
y
plant outputs, and F ,
G , etc., are constant matrices of appropriate size. The ss2step function
calculates the step responses of all the outputs of this process with respect to
all the inputs in u, and puts this information into the variable plant in MPC
step format. The section for mod2step describes the step format in detail.
The input variable tfinal isthetimeatwhichyouwouldliketoendthestep
response calculation, and delt1 is the sampling period. For continuous
systems, use delt1=0. If you do not specify delt1, the default is delt1=0.
The optional input variable delt2 is the desired sampling period for the step
response model. If you use delt2=[ ] or omit it, the default is delt2=delt1 if
delt1 is specified and delt1 neq 0;otherwise,thedefaultisdelt2=1.
The optional input variable nout is the output stability indicator. For stable
systems, set nout equal to the number of outputs, n
y
. For systems with one or
more integrating outputs, nout is a column vector of length n
y
with nout(i)=0
indicating an integrating output and nout(i)=1 indicating a stable output. If
you use nout=[ ] or omit it, the default is nout=n
y
(only stable outputs).
x
·
t() F xt() G ut()+=
yt() Cx t() Du t()+=
ss2step
4-139
Example The following process has 3 inputs and 4 outputs (and is the same one used for
the example in the mod2step section):
phi=diag([0.3, 0.7, -0.7]);
gam=eye(3);
c=[1 0 0; 0 0 1; 0 1 1; 0 1 0];
d=[1 0 0; zeros(3,3)];
The following command duplicates the results obtained with mod2step:
delt1=1.5; tfinal=3
*
1.5;
plant=ss2step(phi,gam,c,d,tfinal,delt1)
See Also plotstep, mod2step, tfd2step
svdfrsp
1-140
1svdfrsp
Purpose Calculates the singular values of a varying matrix, for example, the frequency
response generated by mod2frsp.
Syntax [sigma,omega] = svdfrsp(vmat)
Description vmat is a varying matrix which contains the sampled values F(w
1
), ... , F(w
N
)of
a matrix function F(w ).
If the smaller dimension of F(w
i
)ism,andifs
1
(w
i
), . . . , s
m
(w
i
)arethesingular
values of F(w
i
), in decreasing magnitude, then the output sigma is a matrix of
singular values arranged as follows:
The output omega is a column vector containing the frequencies w
1
,...,w
N
.
Example See mod2frsp, varying format for an example of the use of this function.
See Also mod2frsp
sigma
s
1
w
1
()s
2
w
1
() … s
m
w
1
()
s
1
w
2
()s
2
w
2
() … s
m
w
2
()
s
1
w
N
()s
2
w
N
() … s
m
w
N
()
=
… ……
tfd2mod, tf format
4-141
1tfd2mod, tf format
Purpose tfd2mod converts a transfer function (continuous or discrete) from the MPC tf
format into the MPC mod format, converting to discrete time if necessary.
Syntax model = tfd2mod(delt2,n
y
,g1,g2,g3,...,g25)
Description Consider a transfer function such as
or
The MPC tf format is a matrix consisting of three rows:
The tf matrix will always have at least two columns, since that is the minimum
width of the third row.
The input arguments for tfd2mod are:
row 1 The n coefficients of the numerator polynomial, b
0
to b
n
.
row 2 The n coefficients of the denominator polynomial, a
0
to a
n
.
row 3 column 1: The sampling period. This must be zero if the
coefficients in the above rows are for a continuous system. It
must be positive otherwise.
column 2: The time delay. For a continuous-time transfer
function, it is in time units. For a discrete-time transfer
function, it is the integer number of sampling periods of time
delay.
Gs()
b
0
s
n
b
1
s
n 1–
… b
n
+++
a
0
s
n
a
1
s
n 1–
… a
n
+++
----------------------------------------------------------------=
Gz()
b
0
b
1
z
1–
… b
n
z
n–
+++
a
0
a
1
z
1–
… a
n
z
n–
+++
-------------------------------------------------------------=
tfd2mod, tf format
1-142
delt2
The sampling period for the system. If any of the transfer functions g1, ...,
gN are continuous-time or discrete-time with sampling period not equal to
delt2, tfd2mod will convert them to discrete-time with this sampling period.
ny
The number of output variables in the plant you are modeling.
g1, g2,...gN
A sequence of N transfer functions in the tf format described above, where N ?
1. These are assumed to be the individual elements of a transfer-function
matrix:
Thus it should be clear that N must be an integer multiple (n
u
) of the number
of outputs, n
y
.
Also, tfd2mod assumes that you are supplying the transfer functions in a
column-wise order. In other words, you should first give the n
y
transfer
functions for input 1 (g
1,1
to g
ny
, 1), then the n
y
transfer functions for input 2
(g
1,2
to g
ny
,2),etc.
tfd2mod converts the transfer functions to discrete-time, if necessary, and
combines them to form the output variable, model, which is a composite system
in the MPC mod form.
Example Consider the linear system:
g
11,
g
12,
… g
1 n
u
,
g
21,
g
22,
… g
2 n
u
,
g
n
y
1,
g
n
y
2,
… g
n
y
n
u
,
………
.
.
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
3.8e
8– s
14.9s 1+
------------------------
4.9e
3s–
13.2s 1+
------------------------
ws()+=
tfd2mod, tf format
4-143
The following commands build separate models of the response to the
manipulated variables, u, and the unmeasured disturbance, w,allfora
sampling period T = 3 then combines them using addumd to get a model of the
entire system (the pmod variable):
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2;
umod=tfd2mod(delt,ny,g11,g21,g12,g22);
gw1=poly2tfd(3.8,[14.9 1],0,8);
gw2=poly2tfd(4.9,[13.2 1],0,3);
wmod=tfd2mod(delt,ny,gw1,gw2);
pmod=addumd(umod,wmod);
Restriction The current limit on the number of input transfer functions is N = 25.
See Also mod, poly2tfd, tfd2step
tfd2step
1-144
1tfd2step
Purpose CalculatestheMIMOstepresponseofamodelintheMPCtf format. The
resulting step response is in the MPC step format.
Syntax plant = tfd2step(tfinal,delt2,nout,g1)
plant = tfd2step(tfinal,delt2,nout,g1,...,g25)
Description The input variables are as follows:
tfinal
Truncation time for step response.
delt2
Desired sampling period for step response.
nout
Output stability indicator. For stable systems, this argument is set equal to the
number of outputs, n
y
. For systems with one or more integrating outputs, this
argument is a column vector of length n
y
with nout(i)=0 indicating an
integrating output and nout(i)=1 indicating a stable output.
g1, g2,...gN
AsequenceofN transfer functions in the tf format (see tf format section),
where N ? 1. These are assumed to be the individual elements of a
transfer-function matrix:
Thus it should be clear that N must be an integer multiple (n
u
) of the number
of outputs, n
y
.
tfd2step assumes that you are supplying the transfer functions in a
column-wise order. In other words, you should first give the n
y
transfer
functions for input 1 (g
1,1
to g
ny
, 1), then the n
y
transfer functions for input 2
(g
1,2
to g
ny
,2),etc.
g
11,
g
12,
… g
1 n
u
,
g
21,
g
22,
… g
2 n
u
,
g
n
y
1,
g
n
y
2,
… g
n
y
n
u
,
………
.
.
tfd2step
4-145
The output variableplant isthecalculatedstepresponseofthen
y
outputs with
respect to all inputs. The format is as described in the step section.
Example Consider the linear system:
which is the same as that considered in the mpcsim example. We build the
individual tf format models, then calculate and plot the MIMO step response.
g11=poly2tfd(12.8,[16.7 1],0,1);
g21=poly2tfd(6.6,[10.9 1],0,7);
g12=poly2tfd(-18.9,[21.0 1],0,3);
g22=poly2tfd(-19.4,[14.4 1],0,3);
delt=3; ny=2; tfinal=90;
plant=tfd2step(tfinal,delt,ny,g11,g21,g12,g22,gw1,gw2);
plotstep(plant)
The plots should match the example output in the plotstep description.
Restriction The current limit on the number of input transfer functions is N = 25.
See Also mod2step, plotstep, ss2step
y
1
s()
y
2
s()
12.8e
s–
16.7s 1+
------------------------
18.9e–
3s–
21.0s 1+
-------------------------
6.6e
7s–
10.9s 1+
------------------------
19.4e–
3s–
14.4s 1+
-------------------------
u
1
s()
u
2
s()
=
th2mod, theta format
1-146
1th2mod, theta format
Purpose Converts a SISO or MISO model from the theta format (as used in the System
Identification Toolbox) to one in the MPC mod format. Can also combine such
models to form a MIMO system.
Syntax umod = th2mod(th)
[umod,emod] = th2mod(th1,th2,...,thN)
Description The System Identification Toolbox allows you to identify single-input,
single-output (SISO) and multi-input, single-output (MISO) transfer functions
from data. The MISO form relating an output, y,tom inputs, u
1
to u
m
,anda
noise input, e,is:
where A, B
i
,C,D,andF
i
are polynomials in the forward-shift operator, z.
The System Identification Toolbox automatically stores such models in a
special format, the theta format. See the System Identification Toolbox User’s
Guide for details.
th2mod converts one or more MISO theta models into the MPC mod format,
which you can then use with the MPC Toolbox functions. If you supply a single
input argument, th, and a single output argument, umod,thenumod will model
theresponseofasingleoutput,y,tom inputs, u
1
to u
m
,wherem ? 1. The value
of m depends on the number of inputs included in the input model, th.Note
that umod will reflect the values of the A(z), B(z), and F(z) polynomials in eq. 1.
If you supply a second output argument, emod, it will model the response of y to
the noise, e, i.e., the A(z), C(z) and D(z) polynomials in eq. 1.
If you supply p input models (1 ￡ p ￡ 8), tfd2mod assumes that they define a
MIMO system in the following form:
Az()yk()
B
1
z()
F
1
z()
---------------u
1
k()
B
2
z()
F
2
z()
---------------u
2
k() …
B
m
z()
F
m
z()
---------------- u
m
k()
Cz()
Dz()
------------ ek()+++ +=
A
1
z()y
1
k()
B
11
z()
F
11
z()
----------------- u
1
k() …
B
1m
z()
F
1m
z()
------------------- u
m
k()
C
1
z()
D
1
z()
--------------- e
1
k()++ +=
A
p
z()y
p
k()
B
p1
z()
F
p1
z()
------------------u
1
k() …
B
pm
z()
F
pm
z()
------------------- u
m
k()
C
p
z()
D
p
z()
--------------- e
p
k()++ +=
… … … …
th2mod, theta format
4-147
The p output variables have independent noise inputs. In this case, each of the
p input models must include the same number of inputs, m.Thep outputs are
arranged in parallel in the resulting umod output model (and the emod model, if
used).
If the th models are auto-regressive (i.e., m =0),thenumod will be set to an
emptymatrixandonlyemod will be nonempty.
Example The following commands create three SISO theta models using the mktheta
command (System Identification Toolbox), then converts them to the
equivalent mod form:
th1=mktheta([1 0 -.2],[0 0 -1]);
th2=mktheta([1 -.8 .1],[0 -.5 .3],1,1,1);
th3=mktheta([1 -.2],[0 1],[1 2 0],[1 -1.2 .3], 1);
[umod,emod]=th2mod(th1,th2,th3)
th2mod, theta format
1-148
The results are:
umod =
1.0000 5.0000 1.0000 0 0 3.0000 0
NaN 0 0.2000 0 0 0 1.0000
0 1.0000 0 0 0 0 0
0 0 0 0.8000 -0.1000 0 1.0000
0 0 0 1.0000 0 0 0
0 0 0 0 0 2.0000 1.0000
0 0 -1.0000 0 0 0 0
0 0 0 -0.5000 0.3000 0 0
0 0 0 0 0 1.0000
emod =
Columns 1 through 7
1.0000 7.0000 3.0000 0 0 3.0000 0
NaN 0 0.2000 0 0 0 0
0 1.0000 0 0 0 0 0
0 0 0 0.8000 -0.1000 0 0
0 0 0 1.0000 0 0 0
0 0 0 0 0 1.4000 -0.5400
0 0 0 0 0 1.0000 0
0 0 0 0 0 0 1.0000
0 0 0.2000 0 0 0 0
0 0 0 0.8000 -0.1000 0 0
0 0 0 0 0 3.4000 -0.5400
Columns 8 through 11
0 0 0 0
0 1.0000 0 0
0 0 0 0
0 0 1.0000 0
0 0 0 0
0.0600 0 0 1.0000
0 0 0 0
0 0 0 0
0 1.0000 0 0
0 0 1.0000 0
0.0600 0 0 1.0000
th2mod, theta format
4-149
Restriction The System Identification Toolbox must be installed to use this function.
See Also mod
validmod
1-150
1validmod
Purpose Validates an impulse response model for a new set of data.
Syntax yres = validmod(xreg,yreg,theta)
yres = validmod(xreg,yreg,theta,plotopt)
Description Model validation is a very important part of building a model. For a new set of
data, xreg and yreg, the impulse response model is tested by calculating the
output residual, yres. theta consists of impulse response coefficients as
determined by routines such as plsr or mlr.
Optional input, plotopt, can be supplied to produce various plots. No plot is
produced ifplotoptis equal to 0 which is the default; a plot of the actual output
andthepredictedoutputisproducedifplotopt=1; two plots — plot of actual
and predicted output, and plot of output residual — are produced for
plotopt=2.
Example See plsr for an example of the use of this function.
See Also mlr, plsr
wrtreg
4-151
1wrtreg
Purpose Writes input and output data matrices for a multi-input single-output system
such that they can be used in regression routines mlr and pls for determining
impulse response coefficients.
Syntax [xreg,yreg] = wrtreg(x,y,n)
Description x is the input data of dimension N by n
u
where N is number of data points and
n
u
is number of inputs. y is the output of dimension N by 1. n is number of
impulse response coefficients for all inputs. x is rearranged to produce xreg of
dimension (N–n–1) by n
*
n
u
while yreg is produced by deleting the first n
rows of y. This operation is illustrated as follows:
then
A single sampling delay is assumed for all inputs. y must be a column vector,
i.e., only one output can be specified.
x
x
1
1() … x
n
u
1()
x
1
2() … x
n
u
2()
x
1
N() … x
n
u
N()
=
…
…
y
y 1()
yN()
=
…
xreg
x
1
n() … x
1
1() … x
n
u
n() … x
n
u
1()
x
1
n 1+()…x
1
2() …x
n
u
n 1+()…x
n
u
2()
x
1
N 1–()…x
1
Nn–()…x
n
u
N 1()– … x
n
u
Nn–()
=
… … … …
yreg
yn 1+()
yN()
=
…
wrtreg
1-152
Example See mlr and plsr for examples of the use of this function.
See Also mlr, plsr
I-1
Index
A
abcdchk 4-6
abcdchkm 4-6
addmd 4-4, 4-9, 4-149
addmod 4-4, 4-10
addumd 4-4, 4-11, 4-131, 4-149, 4-155
append 4-13
appmod 4-4, 4-13
associated variable 2-36
autosc 2-7, 2-8, 4-2, 4-14, 4-31
B
balance 4-38
bilinear 3-26
blocking 3-15, 4-60
blocking factor 3-15, 4-76, 4-109
Bode plot 4-97
bound, see constraint
output 4-120
bound, see constraint
manipulated variable 4-97, 4-107, 4-108, 4-110,
4-112
output 4-113, 4-117, 4-120, 4-123, 4-125, 4-128
C
c2d 4-3
c2dmp 4-3, 4-24, 4-89
closed-loop
model in mod format 4-107
system model 4-56
cmpc 2-21, 2-24, 2-29, 4-4, 4-15
complementary sensitivity function 2-18, 4-39
constraint, see bound 2-20, 2-22, 3-12, 3-20
hard 2-20
continuous 3-9
controller gain 4-59, 4-66, 4-81, 4-137
for model in mod format 4-124
covariance 4-25
covariance matrix 4-129
cp2dp 4-3, 4-24
D
d2c 4-3
d2cmp 4-3
dantzgmp 4-6, 4-15, 4-107, 4-126
dareiter 4-6, 4-27
demo file 1-3
dimpulse 4-6
dimpulsm 4-6
discrete 3-9
disturbance
measured 3-10, 3-17, 3-26, 3-33, 4-35
time constant 4-17, 4-66, 4-75, 4-79
unmeasured 3-13, 3-17, 3-20, 3-23, 4-35
disturbance model 4-128, 4-129
disturbance time constant 2-12, 4-131
dlqe2 4-6, 4-25–4-28, 4-133
dlsimm 4-6, 4-119
DMC estimator 4-132, 4-133
DMC, 1, see dynamic matrix control 2-12, 3-18,
3-20, 4-60, 4-76, 4-109, 4-125
dynamic matrix control, 1, see DMC 2-12
E
estimator 4-70
estimator gain 4-112, 4-120, 4-129, 4-130, 4-141
for models in mod format 4-128
Index
I-2
F
fcc_demo.m 2-34
feedforward compensation 4-9
feedforward control 3-33
filter 4-25
fluid catalytic cracking 2-31–2-38
forward-shift operator 3-6
frequency resoponse 2-18
frequency response 4-38, 4-41
plot 4-97
G
gain matrix
for model in mod format 4-135
H
Hessenberg matrix 4-41
horizon 2-11, 2-12
control 3-20, 3-30
moving 2-11
prediction 3-22, 3-30
reciding 2-11
I
IDCOM 1-2
identification 2-6
idle speed control 2-22–2-30
idlectr.m 2-24
imp2step 4-2, 4-29, 4-102
impulse response coefficient 2-6, 4-29
infeasibility 3-23
initial condition 4-79
input horizon 4-60, 4-78, 4-109, 4-125
integrating process 2-7, 2-34
integrating processes 2-2
inverse response 3-29
L
least squares regression 4-30
linear quadratic optimal control 1-2
M
manipulated variable
bound 4-18, 4-68
rate limit 4-68, 4-84, 4-140
saturation limit 4-68, 4-82, 4-140
matrix type 4-63
mean 4-14
measurement noise 2-13, 4-15, 4-25, 4-35, 4-51
minfo 3-29, 4-44, 4-147
minfo vector 4-37
mktheta 4-159
mlr 4-2, 4-14, 4-30, 4-35, 4-37, 4-162, 4-163
mod format 4-35, 4-38, 4-63, 4-107
from discrete-time state-space model 4-146
from model in tf format 4-153
from model in theta format 4-158
matrix type infomation 4-63
mod2frsp 2-18, 3-12, 4-5, 4-38–??, 4-63, 4-152
varying format 4-38
mod2mod 4-3, 4-42
mod2ss 3-10, 4-3, 4-10, 4-28, 4-43, 4-43–??
mod2step 3-29, 4-3, 4-47–??, 4-90
step format 4-47
model algorithmic control 1-2
model-inverse controller 4-60, 4-125
mpcaugss 4-6, 4-51, 4-51–4-53, 4-133
mpccl 4-4, 4-54, 4-54–4-58, 4-61
Index
I-3
mpccon 2-13, 2-15, 2-29, 2-36, 4-4, 4-54, 4-57, 4-59,
4-59–4-62
mpcinfo 4-2, 4-37, 4-63
mpcparal 4-6, 4-9, 4-11, 4-92, 4-148
mpcsim 2-14, 2-29, 2-36, 4-4, 4-65
mpcstair 4-6
mpctut.m 2-3
mpctutid 2-7
mpctutss.m 3-12, 3-20
N
nargchk 4-6
nargchkm 4-6
nlcmpc 4-4, 4-74, 4-74–??
nlmpcdm1.m 4-86
nlmpclib 4-74, 4-81
nlmpcsim 4-4, 4-81
noise filter
time constant 2-12, 4-19, 4-69, 4-79, 4-85
noise model 4-128
nonlinear plant 4-74
O
operating conditions
drive positions 2-22
transmission in neutral 2-22
output
bound 4-19, 4-78, 4-107, 4-108
measured 2-11, 4-25, 4-36, 4-43
unmeasured 2-12, 4-36
output stability indicator 4-29, 4-47, 4-151, 4-156
P
pap_mach.m 3-37, 4-88
paper machine headbox control 3-26
paramod 3-10, 4-4, 4-9, 4-92
parpart4-78
partial least squares 4-100
perfect controller 3-13, 4-60, 4-125
plotall 2-14, 2-15, 3-12, 4-2, 4-93
ploteach 3-12, 4-2, 4-93, 4-95
plotfrsp 2-18, 3-12, 4-2, 4-97
plotstep 2-4, 2-9, 4-2, 4-98
PLS 4-100
pls 4-163
plsr 2-7, 4-2, 4-100, 4-162
pm_lin.m 3-27
pm_nonl.m 3-37
pole
for model in mod format 2-18, 4-135
poly format
conversion to tf format 4-104
poly2tfd2-3, 2-13, 3-6, 3-7, 3-13, 3-20,4-3,
4-21, 4-28, 4-57, 4-104, 4-121, 4-131,
4-142, 4-155, 4-157
prediction horizon 4-59, 4-75, 4-108, 4-125
predictor 4-27
Q
QP 4-21, 4-79, 4-114
quadratic program 4-6, 4-15, 4-107, 4-126
Index
I-4
R
ramp 2-13
rate limit 4-68, 4-82, 4-138
reference value 2-11
regression 4-163
least squares 4-30
partial least squares 4-100
ridge 4-30
rescal 4-2, 4-14
ridge regression 4-30
ringing 3-14
robust 2-22
robustness 2-12
S
sampling period
change in mod format 4-42
saturation constraint 4-65, 4-136
saturation limit 4-68, 4-84, 4-138
scal 2-9, 4-2, 4-14, 4-31
scaling 4-14
scmpc 3-21, 3-24, 3-30, 3-31, 3-32, 3-34, 4-5,
4-107,4-140
scmpcnl 3-37
sensitivity function 2-18, 4-41
sermod 3-10, 4-4, 4-117
setpoint 2-11
signal-to-noise ratio 3-37, 4-130
Simulink 1-3, 3-37, 4-4, 4-81, 4-85, 4-88, 4-89
singular value 2-18, 4-41
of varying matirix 4-152
smpccl 3-12, 3-14, 4-5, 4-118, 4-126
smpccon 3-12, 3-13, 4-5, 4-108, 4-118, 4-124,
4-129, 4-131, 4-134, 4-137, 4-142
smpcest 3-12, 3-18, 4-5, 4-119, 4-128, 4-141
smpcgain 3-12, 4-5, 4-135
smpcpole 3-12, 4-5, 4-121, 4-135
smpcsim 3-12, 3-14, 3-16, 3-18, 3-21, 4-5, 4-113
ss2mod 3-9, 3-29, 4-37, 4-146
ss2moda 4-3
ss2step 2-5, 4-3, 4-150
ss2tf 3-11, 4-3
ss2tf2 4-3
stability 2-12
stability analysis 4-54
stairstep 4-6, 4-93, 4-95
standard deviation 4-14
state estimation 4-15, 4-126
state estimator 3-35, 4-25
state space 2-18, 4-25
step format 4-15, 4-17, 4-29, 4-47–??, 4-55, 4-65
matrix type information 4-63
step response
from mod format 4-47
from state-space model 4-150
from tf format 4-156
plot 4-98
step response coefficient 2-2
step response model 2-2
svdfrsp 3-12, 4-5, 4-152
System Identification Toolbox 3-9, 4-158
system requirement 1-3
T
tf format 4-153
tf2ss 2-5, 4-3
tf2ssm 4-3
tfd2mod 3-4, 3-5, 3-6, 3-7, 3-13, 3-20, 4-3, 4-28,
4-121, 4-142, 4-153, 4-153, 4-159
Index
I-5
tfd2step 2-4, 2-13, 2-24, 4-3, 4-21, 4-57, 4-156
th2mod 3-9, 4-3, 4-9
theta4-158
theta format 3-9, 4-158
time constant
noise filter 4-19, 4-79
unmeasured disturbance 4-69, 4-85, 4-90,
4-113, 4-123
tutorial 1-3
U
usage display 1-3
V
validation 4-162
validmod 4-2, 4-101, 4-162
varying format
matrix type information 4-63, 4-97
varying matrix 4-152
vec2mat 4-6
W
weight 4-16, 4-59, 4-75, 4-108, 4-124
time varying 4-61, 4-126
weighting matrix 2-11, 2-12
white noise 4-25, 4-130
wrtreg 2-7, 2-8, 4-2, 4-163
Index
I-6
| 课件名称： | 自动控制原理（双语教学） |
| 课件分类： | 电气与自动化 |
| 课件类型： | 电子教案 |
| 文件大小： | 46.6MB |
| 下载次数： | 2 |
| 评论次数： | 3 |
| 用户评分： | 5.3 |
- 5. 自动控制原理（双语教学）：05.simulink reference
- 6. 自动控制原理（双语教学）：06.creating GUI
- 7. 自动控制原理（双语教学）：07.Graphics
- 8. 自动控制原理（双语教学）：08.optimization toolbox
- 9. 自动控制原理（双语教学）：09.s function reference
- 10. 自动控制原理（双语教学）：10.S-function
- 11. 自动控制原理（双语教学）：11.MAT-file Format
- 12. 自动控制原理（双语教学）：12.ident
- 13. 自动控制原理（双语教学）：13.mpc simulink library
- 14. 自动控制原理（双语教学）：14.mpc
- 15. 自动控制原理（双语教学）：1.4.1
- 16. 自动控制原理（双语教学）：1.4.2
- 17. 自动控制原理（双语教学）：1.4.3
- 18. 自动控制原理（双语教学）：1.4.4
- 19. 自动控制原理（双语教学）：1.4.5
- 20. 自动控制原理（双语教学）：2.9.1.1
- 21. 自动控制原理（双语教学）：2.9.1.10
- 22. 自动控制原理（双语教学）：2.9.1.2
- 23. 自动控制原理（双语教学）：2.9.1.3
