ModelonDemand Model Predictive Control Toolbox for Simulink
A. Stenman,^{*} M. W. Braun,^{**}
and D. E. Rivera^{**
*}Division of Automatic Control, Department of Electrical Engineering
Linköping University, SE581 83 Linköping, Sweden
^{**}Department of Chemical and Materials
Engineering
Control Systems Engineering Laboratory,
Arizona State University, Tempe, Arizona 852876006
This file is best viewed with Internet Explorer or Netscape Communicator 4.7.
1 Introduction
The ModelonDemand Model Predictive Control (MoDMPC) Toolbox
for Simulink provides researchers and industrial practitioners an easytouse
interface to try out MoD and MoDMPC on both SISO and MIMO nonlinear process
system simulations. This help section is organized as follows:

Introduction

ModelonDemand Estimation and Control: A Brief Review

The MoD Visualization and Validation Toolbox

The MoDMPC Simulink Toolbox
2 ModelonDemand Estimation and Control: A Brief Review
The ModelonDemand concept is very useful from the standpoint
that the user decisions required to determine meaningful parameters are
well guided and simple in comparison to other nonlinear estimation techniques.
Since the ModelonDemand technique is datadriven, the user can spend
less time making decisions regarding model structure and more time in designing
informative dynamical experiments and data analysis. The user should expect
the data length and quality requirements to be much greater in comparison
to linear modeling. Keep in mind however, that other blackbox modeling
techniques such as neural networks also require more data in comparison
to linear modeling. In this section, we will briefly touch on the theory
behind MoD and MoDMPC so the user can understand the effect the userselected
parameters have on the estimate. This helps the user understand the basis
for the suggestions in the context sensitive help and encourages the user
not to blindly "fiddle" with the parameters.
The MoD modeling formulation used in this toolbox follows from the
approach discussed in [9, 10].
The user may refer to published examples of the utility of this modeling
paradigm such as [1, 2]
and others available at
http://www.control.isy.liu.se/
and
http://www.eas.asu.edu/~csel.
2.1 Global vs. Local Modeling
To understand "ModelonDemand" estimation, one must first
understand the concepts of local and global modeling. Suppose we are presented
with a noisy dataset (Y_{i},X_{i})_{i=1}^{N}
and we are looking for some function m to construct a relationship
Y_{i}=m(X_{i})+e_{i},
i=1,...,N,
(1)
where X_{i} represents an independent
variable, Y_{i} represents a response
variable, and e_{i}
depicts a set of independent, identically distributed random variables
with zero means and known variances s_{i}^{2}.
m
is some nonlinear function that maps X_{i}
to Y_{i}. N represents the number
of observations in the database.
One way to estimate the true value y at the current operating
point x is to take the average of the response variables produced
by a small neighborhood of the independent variables around the current
operating point. To try to improve the estimate, it is also possible to
give more weight to data that are closer to the operating point and less
weight to data that are farther away. The weights can be chosen according
to a kernel function, which explicitly defines the shape of the
weights much like the idea of windowing in digital signal processing. What
shape these windows have and how much data should reside in the neighborhood
are issues which dramatically effect the performance of the estimate. The
estimation of the response variable from a small portion of the dataset
is known as local modeling.
In contrast, global modeling first assumes a structure to be used for
the entire dataset. This structure can be based on linear or nonlinear
assumptions, or on firstprinciples knowledge of the system. For example,
linear regression assumes that the function m is of the form
m(X_{i},q_{1},q_{2})=q_{1}+X_{i}^{T}q_{2},
(2)
and that the parameters q_{1}
and q_{2} have
been fit by minimizing the sum of square errors over the entire dataset,



(Y_{i}m(X_{i},q_{1},q_{2}))
(3) 
Such dynamical models as the AutoRegressive with Exogenous input (ARX),
which has a linear structure

A(q^{1})y(k)=q 

B(q^{1})u(k),
(4) 
where

A(q^{1})=1+a_{1}q^{1}+...+a 

q 



B(q^{1})=b_{0}+b_{1}q^{1}+...+b 

q 

, 

are considered global models. Note that q^{1}
is the backward timeshift operator. The idea behind global modeling is
to compress all the information contained within the data into a few parameters.
This produces a computationally simple predictor, at the expense of a biased
estimate, and extensive computational load up front.
2.2 ModelonDemand Concept
Recalling the description of local modeling in the previous
section, the main issues are deciding on how much data to include when
a local estimate is made and what weighting to give the data in that subset.
The MoD concept addresses these concerns and extends the local modeling
idea by taking advantage of modern computer and database technologies.
The MoD estimator performs the data selection and weighting online at
each time instant. Therefore, as the process moves away from the current
operating point, a subset of the database is selected, and a new dynamical
model is constructed. The ModelonDemand estimator is conceptually very
similar to the Artificial Intelligence notion of Lazy Learning. MoD has
the following advantages:

the predictions made by the estimator are datadriven rather than structure
dependent

the predictions are less likely to get trapped by local minima in the optimization
step as compared to global modeling

new data can readily be added to and/or old data can be discarded from
the database
The disadvantages of this approach include:

bruteforce database search for local subsets of data

the estimator suffers from the "curse of dimensionality"

the estimate is dependent on online computational resources

suboptimal estimates of the system can arise when the current operating
condition falls outside of data support
2.3 ModelonDemand Formulation
Potentially, there are many different variants of the MoD
concept that could be used. The exact formulation chosen and the the corresponding
reasoning are the subject of this section. Remaining userdecisions are
discussed as well. In choosing a local model structure, it is desirable
to use one that allows for a closed form, computationally efficient solution
to the local estimation regression problem. By choosing a locally linear
or locally quadratic polynomial, we retain the simple leastsquares solution.
Note that we are now assuming that x and X_{i}
can take on a vector form and they have the same length. Using an l_{2}
norm, we can write
b=arg 



(Y_{i}B^{T}(X_{i}x)b)^{2}w_{i}(x),
(5) 
where
b=(b_{0}b_{1}
...
b_{p})^{T}
(6)
b represents the parameter vector, where p
is dependent on d the length of the X_{i}
vector. For the linear case, p=d. For the quadratic case
p=d+d(d+1)/2.
B(X_{i}x)
represents the vector of regressors for the polynomial model. w_{i}
represents the corresponding smoothing weight. W_{k}(x)
denotes the selected regressor neighborhood. Hence, the linear and quadratic
local models, respectively, can be represented as

B(X_{i}x) 
= 
[1 (X_{i}x)^{T}]^{T} 
(7) 

and 

B(X_{i}x) 
= 
[1 (X_{i}x)^{T}vech^{T}((X_{i}x)(X_{i}x)^{T})]^{T} 

where vech(A) is the operation of obtaining a vector by eliminating
the abovediagonal entries of matrix A and stacking the remaining
entries columnwise.
The MoD formulation allows a natural extension to MIMO prediction.
Additional inputs, outputs, or measured disturbances can be included into
B
to improve the output prediction. Two MoD estimators can be used in parallel
to produce output predictions for two outputs, three for three outputs,
and so on.
Tuning of the weights w_{i} takes
place locally according to a multivariable kernel function of the form
w_{i}(x)=K 
æ
ç
ç
è 

ö
÷
÷
ø 
. (8) 
where K is a kernel function, d is a distance function and
h
is the bandwidth. The specific formulation in this work makes use of a
tricube kernel function. This function is computationally efficient and
has a continuous derivative. The function smoothly descends to zero as
shown in Figure 1. This property reduces high frequency
noise leakage in the estimate as compared to the boxcar and other windows
which do not smoothly descend to zero.
The distance function d is designed to determine how far away the
data in W(k) are from the current operating
point. This distance decides how the data will be weighted by the kernel
function. This distance function can be described as
d(X_{i},x)=(X_{i}x)^{T}M(X_{i}x)
(9)
Defining
M to be equal to the inverse covariance of the regressors
suits the formulation for this paper, since the inverse covariance of the
regressors is readily available and provides an automatic scaling of the
regressors to account for unit and value differences.
The bandwidth h that provides the best tradeoff between bias
and variance errors in the estimate is not known a priori. The MoD
formulation adopted begins by selecting enough data in the local neighborhood
W_{k},
such that the least squares problem is not illconditioned. The neighborhood
size is then expanded until a reasonable tradeoff between bias and variance
errors is reached. Ideally, one seeks to minimize the pointwise Mean Square
Error (MSE)
MSE(m(x,h)) 

E((m(x,h)m(x))^{2},
(10) 
which is expressed on the next line in terms of a squared bias term and
a variance term, respectively.
MSE(m(x,h))=E(m(x,h)m(x))^{2}+Var(m(x,h))
(11)
Recall that m is the true nonlinear relationship between the independent
and dependent variables, while m is the estimated relationship.
Observing the relationship between these two terms and the bandwidth h,
one finds that the variance error of the estimate decreases, while the
bias error of the estimate increases with increasing bandwidth. If the
true function m was known, however, there would be no need to estimate
it. Therefore, localized automatic bandwidth selectors are used to estimate
a goodnessoffit. These selectors are the workhorses in MoD, striking
a balance between bias and variance error.
First, some notation will be defined so that the bandwidth selectors
can be compactly presented. Let
y_{k} 

(Y 

... Y 

)^{T},
(12) 
and define the local design matrix of basis functions as
X_{k} 

æ
ç
ç
ç
ç
ç
ç
è 

ö
÷
÷
÷
÷
÷
÷
ø 

. (13) 
The local weight matrix can be written
W_{k} 

diag(w 

... w 

). (14) 
The solution to the local regression problem and the corresponding estimate
can now be written as

b=(X_{k}^{T}W_{k}X_{k})^{1}X_{k}^{T}W_{k}y_{k} 

m(x,k)=[1 0 ...
0]^{T}b. 
For notational rigor, also note
m(X_{i},k) 

B^{T}(X_{i}x)b.
(15) 
Localized AIC The localized Akaike Information Criterion (AIC)
is written below, where a represents the user
adjustable penalty on the variance term. The AIC criteria has been a popular
method for choosing the number of parameters in linear identification.

AIC(x,k) 



w_{i}(x)(Y_{i}m(X_{i},k))^{2} 


tr(W_{k}) 



x exp( 
a tr((X_{k}^{T}W_{k}X_{k})^{1}(X_{k}^{T}W_{k}^{2}X_{k})) 

tr(W_{k}) 

) 

Akaike's FPE A localized version of Akaike's Final Prediction Error
(FPE) is also available as

FPE(x,k) 



w_{i}(x)(Y_{i}m(X_{i},k))^{2} 


tr(W_{k}) 



x 
2tr(W_{k})+a((X_{x}^{T}W_{k}X_{k})^{1}(X_{k}^{T}W_{k}^{2}X_{k})) 

2tr(W_{k})a((X_{x}^{T}W_{k}X_{k})^{1}(X_{k}^{T}W_{k}^{2}X_{k})) 

. 

Other bandwidth selectors are also included in the toolbox. Now
that the bandwidth selectors are reviewed, selection of candidate values
for bandwidth must be addressed. To bruteforce through the entire dataset,
computing a model and evaluating the goodnessoffit for the successive
addition of one datapoint to the database would require an immense amount
of time. Instead, an exponential search scheme is used as shown in the
following
procedure.

Fit at a very small bandwidth h_{0},
for which the estimation problem is well posed.

Increase bandwidth exponentially according to
h_{i}=C_{h}·h_{i}_{1},
C_{h}=1+ 

(19) 
where C_{h}>1, and d is the
dimension of the regressor space. Compute the corresponding fits. Proceed
until the goodnessoffit falls below a predefined significance level.

Use the model of the bandwidth with lowest goodnessoffit cost in 2.
In actuality, the user can set a range for the number of datapoints used
by the algorithm, k_{min} to
k_{max},
in which the exponential update will search. Part of the validation of
the MoD estimator involves picking appropriate values for this range. The
user can plot the values of k_{opt}
used for each data point in the prediction and inspect the estimated output
for instabilities due to poor data support. Based on this, the user has
a good estimate of how far to increase k_{min}.
In practice, there is no real limit on k_{max}other
than computation time and database length. k_{max}
can be set beyond the highest value for k_{opt}
observed for the validation dataset. Choosing k_{min}
and k_{max} properly and determining
a reasonable regressor structure are the two critical decisions the user
has to make. Regressor structure selection is given it's own treatment
in the next section.
2.4 NARX Regressor Selection
One particularly critical user decision associated with
MoD estimation (or any nonlinear modeling approach for that matter), is
the selection of lagged inputs and outputs to be used as regressors for
the model. This subsection discusses this issue as it relates to MoD, and
supports the use of the best fit global linear ARX model structure as can
be determined in the System Identification Toolbox for MATLAB^{TM}.
2.4.1 Linear ARX
In practice, often the regressor space that provides the
best fit for the global linear ARX modeling approach provides the best
fit for the Model on Demand approach as well. Intuitively, this is attractive
because the best regressor space for the linear case is readily available
by using the parametric order selection option in the System Identification
Toolbox for MATLAB^{TM} for both SISO
and MISO identification. Data can easily be imported into the toolbox.
Pretreatment operations such as filtering, subtraction of means, and estimation/validation
range selection can be performed. When the order selection option is invoked,
the user can visually inspect the model fits and choose between parameterizations
preferred by the Akaike Information Criterion (AIC), Best Fit, or Rissanen's
Minimum Description Length (MDL) approach [6].
The global linear ARX model is in fact a special case of the MoD predictor.
For example, by using a boxcar window instead of a tricube window and forcing
the MoD predictor to use the entire estimation database and only the entire
estimation database, a global linear ARX model would be created at each
step. Therefore, regressors which work well for the global linear ARX case
are likely to be a good starting point for the more general case of MoD
estimation.
3 Model Predictive Control Framework
MPC refers to a class of discretetime control systems which
make explicit use of a model to predict the process output at future time
instants. A receding horizon strategy is used so that at each instant the
horizon is displaced towards the future (see Figure 2).
At each instant the controller calculates a corresponding control sequence
that minimizes a given objective function. This objective function typically
penalizes excessive use of control energy and control error. Only the first
control move calculated is fed to the plant, new information becomes available,
and then the horizon is advanced and the control problem is resolved.
To incorporate the MoD estimator into an MPC control framework, we start
with the MPC objective function given in Meadows and Rawlings [4].
Given the model description and knowledge of the current system state,
we seek a control that minimizes the objective function

J= 

Q_{e}(k)(r(t+k+1)y_{est}(t+k+1))^{2} 
+Q_{u}(k)u^{2}(t+k)+Q 

(k)Du^{2}(t+k)
(20a) 

where Q_{e}(k), Q_{u}(k)
and Q_{Du}(k)
represent penalty weights on the control error, control signal, and control
move size, respectively. r represents the setpoint trajectory, u
the control signal, and y_{est} the
estimated output.
Of the N_{u} future control actions
that minimize J, only the first one is used. Often only 4 or 5 control
increments are solved for, with the implicit assumption that the control
action will be held constant for the remaining moves up to the prediction
horizon N_{p}. This is advantageous
since optimization of Equation 20a can be computationally
intensive for large prediction horizons N_{p}.
The practice of using smaller control horizons N_{u}
has the effect of producing less aggressive controllers and providing stable
control for nonminimum phase systems [4].
A special feature of the formulation is the presence of the move size,
Du(t+k)=u(t+k)u(t+k1),
(20b)
in the objective. In some process control applications, the move size is
restricted due to actuator constraints. The main advantage of MPC is that
hard constraints can be specified by the user on control signal magnitude
and the move size, i.e.,

u_{min}£u(t+k)£u_{min}, 

Du_{min}£Du(t+k)£Du_{min}, 
Most MPC algorithms also make use of a simple disturbance estimate
d(t),
which is defined as the difference of the measured output at the next time
step y(t) and the estimated output y_{est}(t)
based on information at time t1.
d(t)=y(t)y_{est}(t)
the output prediction is then updated with d(t) such that
the estimate over the prediction horizon is
for k=0...N_{p},
y_{est}(t+k)=d(t)+y_{est}(t+k)
More complex output and disturbance models can be included in the MPC framework.
Other common adaptations include the use of the Extended Kalman filter
for online parameter estimation, measured disturbance models, and setpoint
trajectory filters.
3.1 ModelonDemand Model Predictive Control
To adopt the Generalized Predictive Control approach, this
work makes use of a Control ARIMA model (CARIMA), shown in Equation 24,
where D = 1q^{1}.
This forces the controller to have integral action, since implicitly, the
disturbance model is assumed to be nonstationary.
A(q^{1})Dy(t)=B(q^{1})Du(t)
(24)
To express the output prediction at time t+k as a function
of future controls, introduce the following identity

1=A(q^{1})F_{k}(q^{1})D+q^{k}G_{k}(q^{1}) 
i.e.,

y(t+k)=B(q^{1})F_{k}(q^{1})Du(t+k1) 

+G_{k}(q^{1})y(t). 
By partitioning B(q^{1})F_{k}(q^{1})
as
B(q^{1})F_{k}(q^{1})=S_{k}(q^{1})+q^{k}S_{k}(q^{1}),
(25)
where deg S_{k}(q^{1})=k1
and deg
S_{k}(q^{1})=n_{b}2,
the output prediction above can be rewritten as
y(t+k)=S_{k}(q^{1})Du(t+k1)+y(t+k).
(26)
The first term depends on future control actions whereas the remaining
terms depend on measured quantities only. By introducing

y 

(y(t+1) ... y(t+N))^{T}, 


Du 

Du(t) ... Du(t+N1))^{T}, 


y 

(y(t+1) ... y(t+N))^{T}, 


S 


æ
ç
ç
ç
ç
ç
è 
s_{0} 
0 
... 
0 
s_{1} 
s_{0} 
... 
0 
·
·
· 

·
·
· 
·
·
· 
s_{N}_{1} 
s_{N}_{2} 
... 
s_{0} 

ö
÷
÷
÷
÷
÷
ø 

, 

where s_{i} are the coefficients of
S_{k}(q^{1})
and
y=y+SDu
(27)
Therefore
y=y+SDu
(28)
where
Du 

(Du(t) ... Du(t+N_{u}1))^{T},
(29) 
and
The control moves can also be expressed in vector form
u=TDu+u
(31)
where T is a matrix with 1's in the lower triangle
and diagonal, and u is a vector containing only u(t1)'s.
The objective can then be simplified as
where r denotes the desired reference trajectory. For the unconstrained
case, the minimizing control sequence is obtained explicitly by ordinary
least squares theory. For the constrained case, the constraints
can be reformulated in matrix/vector form and the problem is solved efficiently
using standard numerical optimization algorithms.
The full MPC formulation proposed in this toolbox takes the following
form


ry SDu 
 

+ TDu+u 
 

+Du 
 

+ 


(r(t+N_{p})y_{est}(t+N_{p}))^{T}Q_{w}(r(t+N_{p})y_{est}(t+N_{p}))
+e 
 

+Q_{y}_{'}e 

s.t.

I(TDu+u) 
£ 
u_{max} 
(32) 

I(TDu+u) 
£ 
u_{min} 
(33) 

D(TDu+u) 
£ 
u_{max}+u 
(34) 

D(TDu+u) 
£ 
u_{min}u 
(35) 

y+SDu 
£ 
y_{max}+e_{max} 
(36) 

ySDu 
£ 
y_{min}+e_{min} 
(37) 
where
e=[e_{max}^{T}e_{min}^{T}]^{T}
(38)
and
e³0.
(39)
This formulation makes use of a weighted endpoint condition for ensuring
stability, per similar arguments as those of Demircioglu and Clarke [3].
The number of endpoint conditions (conditions appended to time t+N_{p}
) can be adjusted to cover the order of the system, to essentially force
the states of the system to zero, while meeting setpoint at the end of
the prediction horizon. While it is possible to derive a theoretical
value for the penalty Q_{w}, Demircioglu
and Clarke have shown that this penalty may require some adjustment by
the user to minimize overshoot for certain systems. Thus it is left as
a tunable parameter for the user. Additional necessary conditions
for the stability of the closedloop system include; the prediction horizon
should be 2 greater than the order of the system, and the number of endpoint
constraints should be 1 greater than the order of the system.
The soft constraints approach developed by Scokaert and Rawlings
[8] for linear systems, is adopted to provide
intuitive infeasibility handling for nonlinear systems which exhibit nonminimum
phase characteristics (ie. inverse response and/or delay). Two penalties
associated with this approach appear in the MoDMPC cost function. The first
penalty Q_{y} on the 2norm
of the output constraint slack variables, is generally chosen to be an
order of magnitude or more larger than the 1norm penalty, Q_{y}_{'}.
The 1norm exists in the cost function to keep the qp optimization
problem well conditioned. Increasing these penalties has the effect of
demanding a smaller output constraint violation, at the expense of allowing
the violation to persist.
4 The MoD Visualization and Validation Toolbox
The MoD Visualization and Validation Toolbox GUI has been
designed to provide researchers a userfriendly interface for executing
the following:

visualizing data support in the regressor spaces

examining power spectra of the input and output for both the estimation
and validation data

visually inspecting noise and nonlinearities apparent in the ETFE

iterative refinement of the MoD estimator user decisions
The toolbox has been designed to handle the full MIMO problem. Thus
the user can naturally include measured disturbances (additional inputs),
and associated variables (additional outputs), and include this information
in the estimation/validation of the MoD estimator and subsequently in a
MoDMPC Simulink block. The only limitation at this time, is the user cannot
operate the MoDMPC Simulink blocks with anything other than the locally
linear p=1 model structure, although the user can validate other
local models as discussed below. The following section takes the user stepbystep
through each section of the GUI. The V&V Toolbox has context sensitive
help buttons throughout the display which will bring the user back to the
appropriate paragraph in this section.
4.1 Importing Estimation and Validation Data
The MoD Visualization and Validation Toolbox only
accepts data in the System Identification Toolbox for MATLAB^{TM
}Z
format [6]. Therefore, the user must place
both the estimation and validation data in the following format
Z= 
é
ê
ê
ê
ê
ë 

y_{1}(k) 
... 
y_{ny}(k) 
u_{1}(k) 
... 
u_{nu}(k) 
·
·
· 

·
·
· 
·
·
· 

·
·
· 
y_{1}(k+N1) 
... 
y_{ny}(k+N1) 
u_{1}(k+N1) 
... 
u_{nu}(k+N1) 


ù
ú
ú
ú
ú
û 
, (40) 
where n_{y} is the number of outputs
from the system and n_{u} is the number
of inputs to the system. The user must save the estimation and validation
data to separate files in any desired directory. When saving the data in *.mat format, the file name and variable name must be the same.
Having done this, the user can import the data into the GUI, by clicking
on the Load Estimation Data and Load Validation Data buttons
in the GUI and selecting the appropriate paths and filenames.
4.2 Plotting Data in the Time Domain
The user may now specify the number of inputs, outputs
and the sampling time of the data in the GUI. The GUI will check the consistency
of the entered values with the actual sizes of the estimation and validation
datasets, and will stop execution if they are not consistent. The Plot
Time Domain button can be pressed to verify the identity of the data
and compare the estimation and validation datasets in the time domain.
The user can decide to include the validation in all visualization plots
by checking the include validation data in ALL plots option. All
data plotted in the time domain will have the appropriate values for the
time axis. The user is free to modify all plots generated by the GUI for
presentations, manuscript preparation, etc., by using the figure editing
tools.
4.3 Plotting the Data in the Space Domains
The most critical factor in the performance of MoD, is
data support in the regressor spaces. Specifically, if data is not available
in the input, output, and dynamical input/output spaces where the process
is expected to operate, the MoD algorithm can only base it's prediction
on the closest neighborhood of data. For linear and weakly nonlinear systems,
the user may still find MoD prediction surprisingly good. For highly nonlinear,
highly interacting, and/or illconditioned systems, MoD prediction can
be quite poor if the data does not support the operating regimes. As a
result, MoDMPC control can exhibit poor setpoint tracking, excessive manipulated
variable movement, and unstable control of the system.
The user may look at the input, output and input/output spaces by checking
the appropriate option, and pressing Plot Space Domains. The user
may look at differenced data plotted in the spaces selected by pressing
Plot
Dif. Space Domains. This allows the user to visualize the 2^{nd}
order properties of the data (note by 2^{nd}
order, we're referring to state information, not order of nonlinearity).
Normally, the user has little control over the spaces covered by differenced
data in the design of excitation signals. However, different classes of
inputs may cover these spaces differently.
4.4 Plotting Data in the Frequency Domains
The input and output periodograms (power spectra) can
be viewed by selecting the appropriate option. The periodograms are based
on the element by element squaring of the absolute values of the fft of
the data sequence. No windowing, or zero padding is used to reduce variance
or account for nonperiodic data. These plots can be used to verify the
input signal designs, compare input signal designs against validation data
or operational data, or observe the frequency response of the system. The
user who has not removed means from the input/output data can disregard
the dc components of the spectra such that the automatic axis settings
are reasonable for viewing all of the data. The ETFE is simply the elementbyelement
multiplication of the fft of the output multiplied by the elementbyelement
inverse of the fft of the input [7]. The
ETFE provides the user with a simple assessment of the behavior of the
system in the frequency domain. Future versions of this toolbox may include
Coherency analysis to determine whether the process is nonlinear.
4.5 ModelonDemand Crossvalidation
In the MoD Visualization and Validation GUI, the
user may crossvalidate the ModelonDemand estimator, which is using the
estimation dataset for it's database, on validation data. This gives the
user a sense for how well the MoD estimator can ``match'' the output from
the process on new data. This part of the GUI also allows the user to improve
their choice of MoD decisions.
4.6 NARX Structure Matrix
The NARX structure matrix tells the MoD algorithm how
many lagged outputs n_{a}, lagged inputs
n_{b},
and how long to delay the input
n_{k}
in the estimation of the plant output. Since the MoD algorithm will handle
MIMO prediction, this information must be present for each output in the
estimation/validation datasets such that the NARX structure has the same
format as the n_{i} below (note the
y_{i }and u_{i
}are
not part of the matrix passed to the toolbox).


y_{1} 
··· 

u_{1} 
··· 

u_{1} 
··· 

y_{1} 

··· 


··· 


··· 

·
·
· 
·
·
· 

·
·
· 
·
·
· 

·
·
· 
·
·
· 

·
·
· 


··· 


··· 


··· 


(41) 
As an example, consider a 2 input, 2 output MIMO process. If each output
can be described as a function of past inputs and outputs, ie.

y_{1}(k) 
= 
f(y_{1}(k1),u_{1}(k2),u_{1}(k3),u_{2}(k1)) 
(42) 

y_{2}(k) 
= 
f(y_{1}(k1),y_{2}(k1),y_{2}(k2),u_{1}(k1),u_{2}(k1)), 

the appropriate NARX structure matrix for this system would be
Once created, the user may save the NARX structure matrix and call it into
the GUI by pressing the Load NARX Matrix button.
4.7 MoD User Decisions
4.7.1 Prediction Horizon
The user may decide on how far into the future the MoD algorithm
can predict before using actual output data to estimate the next output.
The user may select one step ahead prediction, or any number of steps up
to the length of the validation data, or an infinite prediction horizon
(pure simulation). One step ahead prediction is interesting for the user
considering immediate initial time performance of the algorithm.
The user may also choose a horizon up to the 95% settling time of the system.
Simulation is preferred for those considering the ability of MoD to replicate
steadystate outputs and low frequency characteristics in the data. In
future versions of this toolbox, the user will have the option to specify
a numeric value for the prediction horizon.
4.7.2 kmin and kmax
The user can select the minimum number of data k_{min}
and the maximum number of data k_{max}
to be considered by the MoD algorithm. The user should pick k_{min},
such that the performance of MoD does not ``burst'' due to an illconditioned
local model. The user will be able to pick an appropriate k_{min}
value by inspecting the lower plot in the validation figure.
Initially, k_{max} can be made
the same value as the length of the estimation data. k_{max}
can be reduced to a value above the maximum value observed in the lower
plot in the validation figure. This will reduce the computational load
of the algorithm, while maintaining performance.
4.7.3 Polynomial Order
Local model polynomial orders of p=0, p=1,
and p=2 are supported in the GUI. It is also possible to select
the ``Opt'' option, which will use the MoD optimization approach discussed
in Section 4.5.4 of [10]. Note that the
MoDMPC Simulink blocks can only make use of a MoD choice of p=1
at this time.
4.7.4 Goodnessoffit Measure
The goodnessoffit measures noted in the table below are
supported in this toolbox. Some measures such as the Improved Akaike's
Information Criteria, are relatively new in the statistical literature,
and have not yet received the extensive use as the AIC or FPE criteria
have. However, there are theoretical arguments for why one should try them
[5].
Goodnessoffit Measure

Abbreviation

Specify a?




crossvalidation

CV

no

generalized crossvalidation

GCV

no

Akaike's Final Prediction Error

FPE

yes

Akaike's Information Criteria per Loader

AICN

yes

Akaike's Information Criteria

AIC

yes

Improved Akaike's Information Criteria

AICC

no

Along with the goodnessoffit measure, the user must choose the value
a
which is used to increase the penalty on the variance of the estimate.
a
must be greater than or equal to 2 and typical values are 2 or 3.
4.7.5 Goodnessoffit Optimization
This option determines how the optimal bandwidth at each
time step will be determined. The global option finds the best goodnessoffit
over all bandwidths considered. The polynomial option fits a polynomial
to the values of goodnessoffit and estimates an optimal bandwidth from
the polynomial function.
4.7.6 Include 95% Confidence Intervals in Plots
By checking this box, the user may include 95% confidence
interval bounds in the validation plots.
4.8 The Validate Button
The Validate button uses the current MoD specifications
to compute an estimate of the validation output. A wait bar will appear
which will give the user a sense for when the computations will finish.
Once finished, the GUI will plot the outputs and the corresponding number
of data used at each time step during the prediction. In the plot title,
the user can view the rms and maximum errors between the estimated output
of the MoD estimator and the actual output from the process. The number
of data used by the MoD algorithm at each time step is useful for the user
to decide on reasonable datasize limits k_{min}
and k_{max} for the MoD estimator.
4.9 Export Parameters
The MoD V&V GUI can save all of your MoD parameters,
the filenames and paths of your estimation and validation data, and the
corresponding estimation rms and maximum errors in one structured variable.
This is NOT done automatically, the user must type or select the
filename/variable name in the save as window brought up by pressing
the Save Parameters button.
5 The MoDMPC Simulink Library
The MoDMPC Simulink Library is designed to use the parameters the user
has found to work well in the MoD Visualization and Validation GUI.
The user may then select the weighting and constraints for the control
problem, initialize the controller, and evaluate control performance.
The MoDMPC Simulink blocks can be incorporated in any Simulink file to
control a simulation. Multiple MoDMPC controllers are permissible,
provided that the parameter names are unique. In this section, the
individual MoDMPC controllers in the MoDMPC Simulink Library are discussed.
Each subsection provides detail regarding the structure of the specific
information required to run that particular block. The first subsection
discusses fields common to most blocks.
5.1 Generic MoDMPC Simulink Block Fields
Simulation Time is the final time at which the simulation will end (a.k.a
stop time in Simulink Parameters dialogue). This specification allows
multiple MoDMPC blocks in the same model, since each block will be able
to know a priori how long to make the discrete states vector. The
Parameters Variable Name entered should be the one corresponding to the
parameters structure created with the MoD V&V that the user would like
to use to run the MoDMPC controller. The variable should be available
in the workspace. Both Simulation Time and Parameters Variable Name
must be entered for the block to function.
The penalties for the MoDMPC objective function must be defined appropriately
for the number of inputs and number of outputs the controller is expected
to handle. MoDMPC blocks 1 through 4 can only handle square systems.
The user need not expand the matrices to cover the prediction or control
horizon. So for example, suppose the user has a twoinput/twooutput
system, and would like to penalize the first control error by 1, the second
control error by 2, the first input move suppression by 0.1, and the second
input move suppression by 0.2, the user would enter in [[1 0;0 2] zeros(2)
[0.1 0;0 0.2]].
The output horizon entry N_{p} is
a scalar value that determines the output horizon for all output predictions.
This entry should be in units of integer multiples of the sampling time
of the estimation database used by the MoD estimator. In the
traditional practice of MPC, this value has been selected to cover the
95% sampling time. This is not a necessary condition for good performance
of MoDMPC, and when setpoint trajectories are preprogrammed, performance
may be enhanced by choosing a smaller value for N_{p}.
In fact, if an endpoint condition or endpoint hard constraint is used in
conjunction with a preprogrammed setpoint trajectory, the controller may
make premature changes in the control action if long prediction horizons
are used. Keep in mind, large control horizons will have the effect
of producing a computationally burdensome optimization problem.
The control horizon entry N_{u}
is also a scalar value and remains constant for all inputs. Some
literature suggests that small control horizons have the effect of producing
less aggressive control for nonminimum phase systems. The GPC literature
suggests this value be one larger than the number of unstable modes in
the system. From inspection of the formulation, large control horizons
will also have the effect of producing a computationally burdensome optimization
problem.
Output limits for each output can be specified in the Output Limits
field. The first column of entries corresponds to the lower output
limits for each output. The second column corresponds to the upper
output limits for each output. For example, if the user would like
to constraint two outputs such that 0 £ y_{1
}£
1, and 1 £y_{2}£
2, the user would enter [[0;1] [1;2]].
Input limits are defined in the same way as output limits. As
an example, if the user would like to constraint three inputs, such that
0.5 £ u_{1 }£
2, 8 £ u_{2
}£
12, and 100 £
u_{3
}£
250, the user would enter [[0.5;8;100] [2;12;250]].
Increment limits are defined in the same way as input and output limits.
Rarely will a user require moves in the negative direction to be bounded
differently than in the positive direction, but this option is available
nonetheless.
The reference filter time constant allows the user to filter the setpoint
so as to produce a controller with two degreeoffreedom performance (i.e.
less aggressive setpoint tracking is available, while disturbance rejection
is handled with more aggressive control moves). The reference filter
is defined below.
Thus the reference filter time constant entered will correspond to
t.
Initial state (or initial input/output) vectors are defined as [[y_{1}
; y_{2} ;... ] [u_{1
};
u_{2
};...
]].
The neighborhood size entry [k_{min} k_{max}
] corresponds to the limits on the number of datapoints the localized bandwidth
selector can choose. If no entry is made, the MoDMPC controller will
use the values found in the structure specified in the Parameters variable
name.
Lastly, the user has the option of allowing new inputs and outputs to
be transformed into the appropriate NARX regressor structure and included
in the database. It is possible that by allowing additional data
into the database, the MoD estimate will improve in accuracy. For
systems with unmeasured disturbances, allowing additional "corrupted" data
into the database may produce a less accurate MoD estimate.
5.2 MoDMPC Block 1
MoDMPC Block 1 uses a solitary hard constraint at the end of the prediction
horizon, to force the predicted output y_{est},
to the setpoint r. N_{p}
is usually selected beyond the 95% settling time for this endpoint constraint
to be sufficiently stabilizing. To use the endpoint constraint, select
this option by clicking the check box to create a check mark. This
block can only handle square systems.
For additional field descriptions, see the Generic
MoDMPC Simulink Block Fields section.
5.3 MoDMPC Block 2
MoDMPC Block 2 replaces the endpoint constraint option with an option
to use any number of endpoint conditions n at the end of the prediction
horizon. The value for N_{p}
should be increased by n, as well. Q_{w}
is chosen orders of magnitude larger than the rest of the cost function
penalties to guarantee stability. Q_{w}
also effects performance and may be made smaller to provide less aggressive
control.
For additional field descriptions, see the Generic
MoDMPC Simulink Block Fields section.
5.4 MoDMPC Block 3
MoDMPC Block 3 uses weighted endpoint conditions in the same manner
as MoDMPC Block 2. MoDMPC Block 3 allows the user to specify preprogrammed
setpoint trajectories to provide anticipative action to setpoint changes.
If the controller is run beyond the last specified setpoint value of r,
the MoDMPC Block 2 will hold the setpoint constant for the rest of the
simulation. r must be consistent with the number of outputs to be
controlled, and must be sampled at a constant rate, consistent with the
estimation database sample time.
For additional field descriptions, see the Generic
MoDMPC Simulink Block Fields section.
5.5 MoDMPC Block 4
Output constraint infeasibilities are handled in a 'soft' manner by
MoDMPC Block 4. The same endpoint condition strategy used in MoDMPC
Block 2 is employed. Note MoDMPC Block 4 requires setpoint trajectories
be fed in a nonanticipative manner from the Simulink window.
For additional field descriptions, see the Generic
MoDMPC Simulink Block Fields section.
5.6 MoDMPC Block 5
In comparison, MoDMPC Block 5 is the most computationally intensive
block. MoDMPC Block 5 can handled measured disturbances, and nonsquare
systems. The soft endpoint and soft output constraint infeasibility
handling is managed in the same way as MoDMPC Block 2, and MoDMPC Block
4, respectively. The initial input, initial disturbance, and initial
output variables must be column vectors.
For additional field descriptions, see the Generic
MoDMPC Simulink Block Fields section.
References

[1]

M. W. Braun, D. E. Rivera, and A. Stenman, 2002, Application of Minimum
Crest Factor Multisinusoidal Signals for "PlantFriendly" Identification
of Nonlinear Process Systems. Control Engineering Practice, in press.

[2]

M. W. Braun, D. E. Rivera, and A. Stenman, 2001, A 'ModelonDemand'
identification methodology for nonlinear process systems. International
Journal of Control, Vol.74, No. 18, 17081717.

[3]

H. Demircioglu and D. W. Clarke, 1993, Generalised predictive control
with endpoint state weighting. IEE ProceedingsD, Vol. 140, No.
4, 275282.

[4]

E. S. Meadows and J.B. Rawlings, 1997, Model Predictive Control.
In: Nonlinear Process Control, Eds. M. A. Henson and D. E. Seborg, Prentice
Hall PTR, Upper Saddle River, NJ.

[5]

C. M. Hurvich and J. S. Simonoff, 1998, Smoothing parameter selection
in nonparametric regression using an improved Akaike information criterion.
Journal
of the Royal Statistical Society B, Vol 60, 271293.

[6]

L. Ljung, 1997, MATLAB System Identification Toolbox: User's Guide.
The Mathworks, Inc., Nantucket, MA.

[7]

L. Ljung, 1999, System Identification: Theory for the User.
Prentice Hall PTR, Upper Saddle River, NJ.

[8]

P. O. M. Scokaert and J. B. Rawlings, 1999, Feasibility Issues in
Linear Model Predictive Control. AIChE Journal, Vol. 45, No. 8,
16491659.

[9]

Stenman, A., Gustafsson, F., and Ljung, L., 1996, Justintime models
for dynamical systems. Proceedings of the IEEE Conference on Decision
and Control, Kobe, Japan, 11151120.

[10]

Stenman, A., 1999, Model on Demand: Algorithms, Analysis and Applications.
Dissertation No. 571, Linköping University, SE581 83 Linköping,
Sweden.
This document was translated from L^{A}T_{E}X
by
H^{E}V^{E}A.