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:
 

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,
N
å
i=1
(Yi-m(Xi,q1,q2))     (3)
Such dynamical models as the AutoRegressive with Exogenous input (ARX), which has a linear structure
A(q-1)y(k)=q
-nk
 
B(q-1)u(k),     (4)
where
A(q-1)=1+a1q-1+...+a
 
na
q
-na
 
B(q-1)=b0+b1q-1+...+b
 
nb
q
-nb
 
,
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 disadvantages of this approach include:

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 
 
min
b
 
å
iÎWk(x)
(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 æ
ç
ç
è
d(Xi,x)
h
ö
÷
÷
ø
.     (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
yk
 
=
 
(Y
 
i1
  ...  Y
 
ik
)T,     (12)
and define the local design matrix of basis functions as
Xk
 
=
 
æ
ç
ç
ç
ç
ç
ç
è
B T(X
 
i1
-x)
·
·
·
B T(X
 
ik
-x)
ö
÷
÷
÷
÷
÷
÷
ø
.     (13)
The local weight matrix can be written
Wk
 
=
 
diag(w
 
i1
  ...  w
 
ik
).     (14)
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
m(Xi,k)
 
=
 
BT(Xi-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
 
=
 
 
å
iÎWk(x)
wi(x)(Yi-m(Xi,k))2
tr(Wk)
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
FPE(x,k)
 
=
 
 
å
iÎWk(x)
wi(x)(Yi-m(Xi,k))2
tr(Wk)
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.

  1. Fit at a very small bandwidth h0, for which the estimation problem is well posed.
  2. Increase bandwidth exponentially according to

  3.  
    hi=Ch·hi-1, Ch=1+
    0.3
    d
        (19)
    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.
  4. 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=
Np-1
 
å
k=0
 
Qe(k)(r(t+k+1)-yest(t+k+1))2 +Qu(k)u2(t+k)+Q
 
Du
(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
y
 
=
 
(y(t+1)  ...  y(t+N))T,
Du
 
=
 
Du(t)  ... Du(t+N-1))T,
y
 
=
 
(y(t+1)  ...  y(t+N))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
S
 
=
 
SL,   L
 
=
 
æ
è
I
0
ö
ø
.     (30)
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
J=r-y -SDu
2
 
Qe
+
TDu+u
2
 
Qu
+Du
2
 
Q
 
Du
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

 
min
Du,e
||r-y -SDu ||
2
 
Qe
+ ||TDu+u ||
2
 
Qu
+||Du ||
2
 
Q
 
Du
+
(r(t+Np)-yest(t+Np))TQw(r(t+Np)-yest(t+Np)) +||e ||
2
 
Qy
+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: 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 ···
y
 
ny
u1 ···
u
 
nu
u1 ···
u
 
nu
y1
n
 
a1
···
n
 
a
 
ny
n
 
b1
···
n
 
b
 
nu
n
 
k1
···
n
 
k
 
nu
·
·
·
·
·
·
  ·
·
·
·
·
·
  ·
·
·
·
·
·
  ·
·
·
y
 
ny
n
 
a1
···
n
 
a
 
ny
n
 
b1
···
n
 
b
 
nu
n
 
k1
···
n
 
k
 
nu
    (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
é
ë
1 0 2 1 2 1
1 2 1 1 1 1
ù
û
.     (43)
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.

 

F(z)=
(1-t)z-1
1-tz-1
    (44)


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 rNp 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.