Model-on-Demand 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, SE-581 83 Linköping, Sweden
**Department of Chemical and Materials
Engineering
Control Systems Engineering Laboratory,
Arizona State University, Tempe, Arizona 85287-6006
This file is best viewed with Internet Explorer or Netscape Communicator 4.7.
1 Introduction
The Model-on-Demand Model Predictive Control (MoDMPC) Toolbox
for Simulink provides researchers and industrial practitioners an easy-to-use
interface to try out MoD and MoDMPC on both SISO and MIMO nonlinear process
system simulations. This help section is organized as follows:
-
Introduction
-
Model-on-Demand Estimation and Control: A Brief Review
-
The MoD Visualization and Validation Toolbox
-
The MoDMPC Simulink Toolbox
2 Model-on-Demand Estimation and Control: A Brief Review
The Model-on-Demand 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 Model-on-Demand technique is data-driven, 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 black-box 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 user-selected
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 "Model-on-Demand" estimation, one must first
understand the concepts of local and global modeling. Suppose we are presented
with a noisy dataset (Yi,Xi)i=1N
and we are looking for some function m to construct a relationship
Yi=m(Xi)+ei,
i=1,...,N,
(1)
where Xi represents an independent
variable, Yi represents a response
variable, and ei
depicts a set of independent, identically distributed random variables
with zero means and known variances si2.
m
is some nonlinear function that maps Xi
to Yi. 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 first-principles knowledge of the system. For example,
linear regression assumes that the function m is of the form
m(Xi,q1,q2)=q1+XiTq2,
(2)
and that the parameters q1
and q2 have
been fit by minimizing the sum of square errors over the entire dataset,
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
|
|
|
| B(q-1)=b0+b1q-1+...+b |
|
q |
|
, |
|
are considered global models. Note that q-1
is the backward time-shift 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 Model-on-Demand 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 on-line 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 Model-on-Demand 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 data-driven 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:
-
brute-force database search for local subsets of data
-
the estimator suffers from the "curse of dimensionality"
-
the estimate is dependent on on-line computational resources
-
suboptimal estimates of the system can arise when the current operating
condition falls outside of data support
2.3 Model-on-Demand 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 user-decisions 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 least-squares solution.
Note that we are now assuming that x and Xi
can take on a vector form and they have the same length. Using an l2
norm, we can write
| b=arg |
|
|
|
(Yi-BT(Xi-x)b)2wi(x),
(5) |
where
b=(b0b1
...
bp)T
(6)
b represents the parameter vector, where p
is dependent on d the length of the Xi
vector. For the linear case, p=d. For the quadratic case
p=d+d(d+1)/2.
B(Xi-x)
represents the vector of regressors for the polynomial model. wi
represents the corresponding smoothing weight. Wk(x)
denotes the selected regressor neighborhood. Hence, the linear and quadratic
local models, respectively, can be represented as
|
B(Xi-x) |
= |
[1 (Xi-x)T]T |
(7) |
|
and |
|
B(Xi-x) |
= |
[1 (Xi-x)TvechT((Xi-x)(Xi-x)T)]T |
|
where vech(A) is the operation of obtaining a vector by eliminating
the above-diagonal entries of matrix A and stacking the remaining
entries column-wise.
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 wi takes
place locally according to a multivariable kernel function of the form
| wi(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(Xi,x)=(Xi-x)TM(Xi-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
Wk,
such that the least squares problem is not ill-conditioned. The neighborhood
size is then expanded until a reasonable tradeoff between bias and variance
errors is reached. Ideally, one seeks to minimize the point-wise 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 goodness-of-fit. 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
and define the local design matrix of basis functions as
| Xk |
|
æ
ç
ç
ç
ç
ç
ç
è |
|
ö
÷
÷
÷
÷
÷
÷
ø |
|
. (13) |
The local weight matrix can be written
The solution to the local regression problem and the corresponding estimate
can now be written as
|
b=(XkTWkXk)-1XkTWkyk |
|
m(x,k)=[1 0 ...
0]Tb. |
For notational rigor, also note
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.
|
|
|
| x exp( |
| a tr((XkTWkXk)-1(XkTWk2Xk)) |
|
|
| tr(Wk) |
|
) |
|
Akaike's FPE A localized version of Akaike's Final Prediction Error
(FPE) is also available as
|
|
|
| x |
| 2tr(Wk)+a((XxTWkXk)-1(XkTWk2Xk)) |
|
|
| 2tr(Wk)-a((XxTWkXk)-1(XkTWk2Xk)) |
|
. |
|
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 brute-force through the entire dataset,
computing a model and evaluating the goodness-of-fit 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 h0,
for which the estimation problem is well posed.
-
Increase bandwidth exponentially according to
where Ch>1, and d is the
dimension of the regressor space. Compute the corresponding fits. Proceed
until the goodness-of-fit falls below a predefined significance level.
-
Use the model of the bandwidth with lowest goodness-of-fit cost in 2.
In actuality, the user can set a range for the number of datapoints used
by the algorithm, kmin to
kmax,
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 kopt
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 kmin.
In practice, there is no real limit on kmaxother
than computation time and database length. kmax
can be set beyond the highest value for kopt
observed for the validation dataset. Choosing kmin
and kmax 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 MATLABTM.
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 MATLABTM 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 discrete-time 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= |
|
Qe(k)(r(t+k+1)-yest(t+k+1))2 |
+Qu(k)u2(t+k)+Q |
|
(k)Du2(t+k)
(20a) |
|
where Qe(k), Qu(k)
and QDu(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 yest the
estimated output.
Of the Nu 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 Np. This is advantageous
since optimization of Equation 20a can be computationally
intensive for large prediction horizons Np.
The practice of using smaller control horizons Nu
has the effect of producing less aggressive controllers and providing stable
control for non-minimum phase systems [4].
A special feature of the formulation is the presence of the move size,
Du(t+k)=u(t+k)-u(t+k-1),
(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.,
|
umin£u(t+k)£umin, |
|
Dumin£Du(t+k)£Dumin, |
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 yest(t)
based on information at time t-1.
d(t)=y(t)-yest(t)
the output prediction is then updated with d(t) such that
the estimate over the prediction horizon is
for k=0...Np,
yest(t+k)=d(t)+yest(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 Model-on-Demand 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 = 1-q-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)Fk(q-1)D+q-kGk(q-1) |
i.e.,
|
y(t+k)=B(q-1)Fk(q-1)Du(t+k-1) |
|
+Gk(q-1)y(t). |
By partitioning B(q-1)Fk(q-1)
as
B(q-1)Fk(q-1)=Sk(q-1)+q-kSk(q-1),
(25)
where deg Sk(q-1)=k-1
and deg
Sk(q-1)=nb-2,
the output prediction above can be rewritten as
y(t+k)=Sk(q-1)Du(t+k-1)+y(t+k).
(26)
The first term depends on future control actions whereas the remaining
terms depend on measured quantities only. By introducing
|
|
|
| Du |
|
Du(t) ... Du(t+N-1))T, |
|
|
|
|
| S |
|
|
æ
ç
ç
ç
ç
ç
è |
| s0 |
0 |
... |
0 |
| s1 |
s0 |
... |
0 |
·
·
· |
|
·
·
· |
·
·
· |
| sN-1 |
sN-2 |
... |
s0 |
|
ö
÷
÷
÷
÷
÷
ø |
|
, |
|
where si are the coefficients of
Sk(q-1)
and
y=y+SDu
(27)
Therefore
y=y+SDu
(28)
where
| Du |
|
(Du(t) ... Du(t+Nu-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(t-1)'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 re-formulated 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
|
|
|
||r-y -SDu |
|| |
|
+ ||TDu+u |
|| |
|
+||Du |
|| |
|
+ |
|
|
| (r(t+Np)-yest(t+Np))TQw(r(t+Np)-yest(t+Np))
+||e |
|| |
|
+Qy'e |
|
s.t.
|
I(TDu+u) |
£ |
umax |
(32) |
|
-I(TDu+u) |
£ |
umin |
(33) |
|
D(TDu+u) |
£ |
umax+u |
(34) |
|
-D(TDu+u) |
£ |
umin-u |
(35) |
|
y+SDu |
£ |
ymax+emax |
(36) |
|
-y-SDu |
£ |
-ymin+emin |
(37) |
where
e=[emaxTeminT]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+Np
) 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 Qw, 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 closed-loop 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 Qy on the 2-norm
of the output constraint slack variables, is generally chosen to be an
order of magnitude or more larger than the 1-norm penalty, Qy'.
The 1-norm 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 user-friendly 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 step-by-step
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 MATLABTM
Z
format [6]. Therefore, the user must place
both the estimation and validation data in the following format
| Z= |
é
ê
ê
ê
ê
ë |
|
| y1(k) |
... |
yny(k) |
u1(k) |
... |
unu(k) |
·
·
· |
|
·
·
· |
·
·
· |
|
·
·
· |
| y1(k+N-1) |
... |
yny(k+N-1) |
u1(k+N-1) |
... |
unu(k+N-1) |
|
|
ù
ú
ú
ú
ú
û |
, (40) |
where ny is the number of outputs
from the system and nu 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 ill-conditioned 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 2nd
order properties of the data (note by 2nd
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 non-periodic 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 element-by-element
multiplication of the fft of the output multiplied by the element-by-element
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 Model-on-Demand Crossvalidation
In the MoD Visualization and Validation GUI, the
user may crossvalidate the Model-on-Demand 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 na, lagged inputs
nb,
and how long to delay the input
nk
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 ni below (note the
yi and ui
are
not part of the matrix passed to the toolbox).
|
| |
y1 |
··· |
|
u1 |
··· |
|
u1 |
··· |
|
| y1 |
|
··· |
|
|
··· |
|
|
··· |
|
·
·
· |
·
·
· |
|
·
·
· |
·
·
· |
|
·
·
· |
·
·
· |
|
·
·
· |
|
|
|
··· |
|
|
··· |
|
|
··· |
|
|
(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.
|
y1(k) |
= |
f(y1(k-1),u1(k-2),u1(k-3),u2(k-1)) |
(42) |
|
y2(k) |
= |
f(y1(k-1),y2(k-1),y2(k-2),u1(k-1),u2(k-1)), |
|
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
steady-state 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 kmin
and the maximum number of data kmax
to be considered by the MoD algorithm. The user should pick kmin,
such that the performance of MoD does not ``burst'' due to an ill-conditioned
local model. The user will be able to pick an appropriate kmin
value by inspecting the lower plot in the validation figure.
Initially, kmax can be made
the same value as the length of the estimation data. kmax
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 Goodness-of-fit Measure
The goodness-of-fit 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].
|
Goodness-of-fit 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 goodness-of-fit 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 Goodness-of-fit Optimization
This option determines how the optimal bandwidth at each
time step will be determined. The global option finds the best goodness-of-fit
over all bandwidths considered. The polynomial option fits a polynomial
to the values of goodness-of-fit 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 kmin
and kmax 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 two-input/two-output
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 Np 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 Np.
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 Nu
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 non-minimum 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 £ y1
£
1, and 1 £y2£
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 £ u1 £
2, 8 £ u2
£
12, and 100 £
u3
£
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 degree-of-freedom 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 [[y1
; y2 ;... ] [u1
;
u2
;...
]].
The neighborhood size entry [kmin kmax
] 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 yest,
to the setpoint r. Np
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 Np
should be increased by n, as well. Qw
is chosen orders of magnitude larger than the rest of the cost function
penalties to guarantee stability. Qw
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 non-anticipative 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 non-square
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 "Plant-Friendly" Identification
of Nonlinear Process Systems. Control Engineering Practice, in press.
-
[2]
-
M. W. Braun, D. E. Rivera, and A. Stenman, 2001, A 'Model-on-Demand'
identification methodology for non-linear process systems. International
Journal of Control, Vol.74, No. 18, 1708--1717.
-
[3]
-
H. Demircioglu and D. W. Clarke, 1993, Generalised predictive control
with end-point state weighting. IEE Proceedings-D, Vol. 140, No.
4, 275--282.
-
[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, 271-293.
-
[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,
1649--1659.
-
[9]
-
Stenman, A., Gustafsson, F., and Ljung, L., 1996, Just-in-time models
for dynamical systems. Proceedings of the IEEE Conference on Decision
and Control, Kobe, Japan, 1115--1120.
-
[10]
-
Stenman, A., 1999, Model on Demand: Algorithms, Analysis and Applications.
Dissertation No. 571, Linköping University, SE-581 83 Linköping,
Sweden.
This document was translated from LATEX
by
HEVEA.