y( k) =

s i v ( k i ) + s n v ( k n 1 )
Step response models can be used for both stable and integrating processes. For an integrating process it is assumed that the slope of the response remains constant after n steps, i.e., sn sn 1 = sn + 1 sn = sn + 2 sn + 1 =. For a multi-input multi-output (MIMO) process with nv inputs and ny outputs, one obtains a series of step response coefficient matrices s 1, 1, i Si = s 2, 1, i s n , 1, i s n , 2, i y y s 1, 2, i s 1, n v, i
where sl,m,i is the ith step response coefficient relating the mth input to the lth output.

sn , n , i

The MPC Toolbox stores step response models in the following format:

S1 S2 Sn

plant = nout(1) nout(2) nout ( n ) y
where delt2 is the sampling time and the vector nout indicates if a particular output is integrating or not:
nout(i) = 0 if output i is integrating. nout(i) = 1 if output i is stable.
The step response can be obtained directly from identification experiments, or generated from a continuous or discrete transfer function or state-space model. For example, if the discrete system description (sampling time T = 0.1) is y(k) = 0.5y(k 1) + v(k 3) then the transfer function is z g ( z ) = ------------------------+ 0.5z

(n ny + ny + 2 ) nv

The following commands (see mpctut.m) generate the step response model for this system and plot it:
num = 1; den = [1 0.5]; delt1 = 0.1; delay = 2; g = poly2tfd(num,den,delt1,delay); % Set up the model in tf format tfinal = 1.6; delt2 = delt1; nout = 1; plant = tfd2step(tfinal,delt2,nout,g); % Calculate the step response plotstep(plant) % Plot the step response

u1 step response : y1 1

0.8 TIME
Alternatively, we could first generate a state-space description applying the command tf2ss and then generate the step response with ss2step. In this case, we need to pad the numerator and denominator polynomials to account for the time delay.
num = [0 num]; den = [den 0 0]; [phi,gam,c,d] = tf2ss(num,den); % Convert to state-space plant = ss2step(phi,gam,c,d,tfinal,delt1,delt2,nout); % Calculate step response
We can get some information on the contents of a matrix in the MPC Toolbox via the command mpcinfo. For our example, mpcinfo(plant) returns:
This is a matrix in MPC Step format. sampling time = 0.1 number of inputs = 1 number of outputs = 1 number of step response coefficients = 16 All outputs are stable.

Model Identification

The identification routines available in the MPC Toolbox are designed for multi-input single-output (MISO) systems. Based on a historical record of the output yl(k) and the inputs v1(k); v2(k),. , v n (k),

5 y-5 -y1

u5 Time 10
Figure 2-5 Responses to a Unit Torque Disturbance for Plant #1 (nominal model = Plant #1)

-2 -4 0

Outputs y2 yTime 10

3 u1 u2

Figure 2-6 Responses to a Unit Torque Disturbance for Plant #2 (nominal model = Plant #1)
As expected, the nominal performance for both Plant #1 and Plant #2 has deteriorated when compared to the simulations shown in Figure 2-1 and Figure 2-2. A similar effect would be observed if we had detuned the controller which uses Plant #2 as the nominal model.
Simulation #4. The parameter values are the same as in Simulation #3. Instead of using cmpc, we use mpccon and mpcsim for simulating the closed loop responses. Figure 2-2 compares the responses for Plant #1 using mpccon and mpcsim, and cmpc. As we can see, for this example and these tuning parameters, the improvement obtained through the on-line optimization in cmpc is small. However, the difference could be large, especially for ill-conditioned systems and other tuning parameters. For example, by reducing the output horizon to P = 80 while keeping the other parameters the same, the responses for Plant # 1 found with mpccon and mpcsim are significantly slower than those obtained with cmpc (Figure 2-8).
-1 Output -2 y1 -3 -4 -5 -Time
Figure 2-7 Comparison of Responses From cmpc, and mpccon and mpcsim for Plant #1 P = inf
Solid: cmpc Dashed: mpccon and mpcsim
-1 Output -2 -3 -4 -5 -6 -Time
Figure 2-8 Comparison of Responses From cmpc, and mpccon and mpcsim for Plant #1 (P = 80)
Application: Control of a Fluid Catalytic Cracking Unit
Fluid Catalytic Cracking Units (FCCUs) are widely used in the petroleum refining industry to convert high boiling oil cuts (of low economic value) to lighter more valuable hydrocarbons including gasoline. Cracking refers to the catalyst enhanced thermal breakdown of high molecular weight hydrocarbons into lower molecular weight materials. A schematic of the FCCU studied4 is given in Figure 2-9. Fresh feed is contacted with hot catalyst at the base of the riser and travels rapidly up the riser where the cracking reactions occur. The desirable products of reaction are gaseous (lighter) hydrocarbons which are passed to a fractionator and subsequently to separation units for recovery and purification. The undesirable byproduct of cracking is coke which is deposited on the catalyst particles, reducing their activity. Catalyst coated with coke is transported to the regenerator section where the coke is burned off thereby restoring catalytic activity and raising catalyst temperature. The regenerated catalyst is then transported to the riser base where it is contacted with more fresh feed. Regenerated catalyst at the elevated temperature provides the heat required to vaporize the fresh feed as well as the energy required for the endothermic cracking reaction.

In many applications, all outputs are measured. In some cases, however, one has nym measured and nyu unmeasured outputs in y, where nym + nyu = ny. If so, the MPC Toolbox assumes that the y vector and the C and D matrices are arranged such that the measured outputs come first, followed by the unmeasured outputs.

Mod Format

The MPC Toolbox works with state-space models in a special format, called the mod format. The mod format is a single matrix that contains the state-space , , C, and D matrices plus some additional information (see mod format in the "Command Reference" chapter for details). The MPC Toolbox includes a number of commands that make it easy to generate models in the mod format. The following sections illustrate the use of these commands.
SISO Continuous-Time Transfer Function to Mod Format
The MPC Toolbox uses a format called the tf format. Let the continuous-time transfer function be + + bn Td s b0s + b1s G ( s ) = --------------------------------------------------------------- e n n1 a0 s + a 1 s + + an where Td is the time delay. The tf format is a matrix consisting of three rows: row 1: row 2: row 3: The n coefficients of the numerator polynomial, b0 to bn. The n coefficients of the denominator polynomial, a0 to an. column 1: The sampling period. This must be zero for a continuous system. (It must be positive for discrete transfer functions see next section). column 2: The time delay in time units. It must satisfy Td 0. The tf matrix will always have at least two columns, since that is the minimum width of the third row.
You can either define a model in the tf format directly or use the command poly2tfd. The general form of this command is
g = poly2tfd(num,den,delt,delay)
For example, consider a SISO system modeled by the transfer function
5.3s 13.6s + 1 G ( s ) = ---------------------------------------------------e 2 54.3 s + 113.5s + 1
To create the tf format directly you could use the command
G = [0 -13.6 1; 54.3 113.5 1; 0 5.3 0];
which defines a matrix consisting of three rows and three columns. Note that all rows must have the same number of columns so you must be careful to insert zeros where appropriate. The poly2tfd command is more convenient since it does that for you automatically:
G = poly2tfd([-13.6 1],[54.3 113.5 1],0,5.3);
Either command would define a variable G in your workspace, containing the matrix

0 54.3000 0

-13.6000 113.5000 5.3000

1.0000 1.0000 0

To convert this to the mod format, use the command tfd2mod, which has the form
model = tfd2mod(delt,ny,g1,g2,g3,.,gN)

where:

The sampling period. tfd2mod will convert your continuous time transfer function(s) g1,., gN to discrete-time using this sampling period. is the number of output variables in the plant you are modeling. A sequence of N transfer functions in the tf format, where N 1. tfd2mod assumes that these are the individual elements of a transfer-function matrix:

3. If a measurement were known to be inaccurate, its R(i,i) value should be relatively large. 4. You can see how serious the degradation would be by setting R(2,2) to a large value, e.g.,10000.
of the estimator simplifies the design procedure. It requires only a rough idea of the characteristic times for the disturbances, and the signal-to-noise ratio for each output. For example, you can verify that the following design rejects the d disturbance almost as well as the optimal Kalman design:
% Alternative estimator design -- output disturbances taus = [5]; signoise = [10]; [Kest, newmod] = smpcest(imod,taus,signoise); % Simulation using scmpc -- no model error [y,u,ym] = scmpc(pmod,newmod,ywt,uwt,M,P,tend,. setpts,ulim,ylim,Kest,z,v,d); plotall(y,u,dt)

MPC of Nonlinear Plant

We are now ready to test the controller design on the real (nonlinear) plant. A special version of the scmpc function (called scmpcnl) is available for this purpose. It uses a nonlinear plant model in the S-function format required by Simulink. (See the Simulink documentation for more information on how to write such models.) The model of the paper machine is in the file pap_mach.m. Simulations with Simulink involving nonlinear models usually take much longer (by an order of magnitude) than linear simulations of a plant of comparable complexity. This is especially likely if the plant model is in the form of an.M file, as is the case here. If such models are to be used extensively, it may be worthwhile to code them as a.mex file (see MATLAB documentation). To see how well the MPC design rejects the d disturbance of Figure 3-8, we could use the commands found in the file pm_nonl.m in the directory mpcdemos. The only differences between these commands and those for the original linear simulation are: We have defined the initial values of the plant state and manipulated variables (x0 and u0, respectively). A step size for numerical integration has been specified. The value of 0.05 minutes provides reasonable accuracy in this application. In general, one must choose the step size to fit the problem (or use a variable step-size integration method, as provided by Simulink). You can verify that the results are nearly identical to those shown in Figure 3-8. In other words, the nonlinearities in the plant have caused negligible

addumd converts pmod and dmod into their state-space form, then uses the mpcparal function to build the composite model.
pmod and dmod must have been created with equal sampling periods and number of output variables. pmod must not include unmeasured disturbances, i.e., its mod format must specify nw = 0. All inputs to dmod must be classified as manipulated variables. (They will be reclassified automatically as unmeasured disturbances in model.) So the mod format of dmod must specify nd = nw = 0 (which is the default for all model creation functions).
addmod, addmd, appmod , paramod, sermod, smpcest

appmod

4appmod
Appends two models to form a composite model that retains the inputs and outputs of the original models. In other words, for models in the MPC mod format appmod replaces the append function of the Control Toolbox.

pmod = appmod(mod1,mod2)

The two input models combine as shown in the following block diagram: u1 d1 w1 mod1 y1 u1 u2 d1 u2 d2 w2 mod2 y2 d2 w1 w2 pmod y1 y2
mod1, mod2 and pmod are in the MPC mod format (see mod in the online
MATLAB Function Reference for a detailed description). You would normally create mod1 and mod2 using either the tfd2mod, ss2mod, or th2mod function.

Restriction See Also

mod1 and mod2 must have been created with equal sampling periods. addmod, addmd, addumd, paramod, sermod

autosc, scal, rescal

Purpose Syntax

4autosc, scal, rescal

Scales a matrix automatically or by specified mean and standard deviation.
[ax,mx,stdx] = autosc(x) sx = scal(x,mx) sx = scal(x,mx,stdx) rx = rescal(x,mx) rx = rescal(x,mx,stdx) autosc scales an input matrix or vector x by its column means (mx) and standard deviations (stdx) automatically and outputs mx and stdx as options. By using scal , the input can also be scaled by some specified means and/or

Description

standard deviations. rescal converts scaled data back to original data. Output mx is a row vector containing the mean value for each column of x while stdx is a row vector containing the standard deviation for each column. Outputs ax and sx are obtained by dividing the difference of each column of x and the mean for the column by the standard deviation for the column, i.e., ax(:,i) = (x :,i) mx(i)/stdx(i). Output rx is determined by multiplying each column of x by the corresponding standard deviation and adding the corresponding mean to that product. If only two arguments are specified in scal or rescal , x is scaled by specified means (mx) only.

Example See Also

See mlr for an example of the use of these functions.

mlr, plsr, wrtreg

Simulates closed-loop systems with hard bounds on manipulated variables and/or outputs using models in the MPC step format. Solves the MPC optimization problem by quadratic programming.

Algorithm

The controller gain is a component of the solution to the optimization problem:

p ny 2

Minimize J ( k ) =
( ywti ( j ) [ ri ( k + j ) y i ( k + j ) ] )

j = 1i = 1 nb nu

( uwt i ( j )ui ( j ) )

j = 1i = 1

with respect to u i ( j ) (a series of current and future moves in the manipulated variables), where y i (k + j) is a prediction of output i at a time j sampling periods into the future (relative to the current time, k), which is a function of
u i ( j ) , ri(k + j) is the corresponding future setpoint, and nb is the number of blocks or moves of the manipulated variables.
Consider the linear system: 18.9e 12.8e ----------------------- ------------------------- u ( s ) 16.7s + 1 21.0s + = 7 s 3s u ( s ) y2 ( s ) 2 19.4e 6.6e ----------------------- ------------------------10.9s + 1 14.4s + 1 y1 ( s ) See the mpccl example for the commands that build the model and a simple controller for this process. Here is a slightly more complex design with blocking and time-varying weights on the manipulated and output variables:
P=6; M=[2 4]; uwt=[1 0; 0 1]; ywt=[1 0.1; 0.8 0.1; 0.1 0.1]; Kmpc=mpccon(model,ywt,uwt,M,P); tend=30; r=[1 0]; [y,u]=mpcsim(plant,model,Kmpc,tend,r);
There is no particular rationale for using time varying weights in this case it is only for illustration. The manipulated variables will make 2 moves during the prediction horizon (see value of M, above). The uwt selection gives u1 a unity weight and u2 a zero weight for the first move, then switches the weights for the second move. If there had been any additional moves they would have had the same weighting as the second move.
The ywt value assigns a constant weight of 0.1 to y2, and a weight that decreases over the first 3 periods to y1. The weights for periods 4 to 6 are the same as for period 3. The resulting closed-loop (servo) response is:

0.5 0.4 0.3 0.2 0.1 0

cmpc, mpccl, mpcsim

mpcinfo

Purpose Syntax: Description

4mpcinfo

Determines the type of a matrix and returns information about the matrix.
mpcinfo(mat) mpcinfo returns information about the type and size of the matrix, mat. The information is determined from the matrix structure. The matrix types include MPC step format, MPC mod format, varying format and constant. mpcinfo returns text output to the screen.
If the matrix is in MPC step format, the output includes the sampling time used to create the model, number of inputs, number of outputs and number of step response coefficients; it also indicates which outputs are stable and which outputs are integrating. If the matrix is in MPC mod format, the output includes the sampling time used to create the model, number of states, number of manipulated variable inputs, number of measured disturbances, number of unmeasured disturbances, number of measured outputs and number of unmeasured outputs. For a matrix in varying format, as formed in mod2frsp, the number of independent variable values, and the number of rows and number of columns of each submatrix are output. For a constant matrix, the text output consists of the number of rows and number of columns.
1 MPC step format: After running the mod2step example mpcinfo(plant)

returns:

This is a matrix in MPC Step format. sampling time = 1.5 number of inputs = 3 number of outputs = 4 number of step response coefficients = 3 All outputs are stable.
2 MPC mod format: After running the ss2mod example mpcinfo(pmod)
This is a matrix in MPC Mod format. minfo = [0 ] sampling time = 2 number of states = 3 number of manipulated variable inputs = 1 number of measured disturbances = 1 number of unmeasured disturbances = 1 number of measured outputs = 1 number of unmeasured outputs = 0
3 varying format: After running the mod2frsp example mpcinfo(eyefrsp)
varying: 30 pts 2 rows 2 cols
mod, step, mod2frsp , varying format

4mpcsim

Simulates closed-loop systems with saturation constraints on the manipulated variables using models in the MPC step format. Can also be used for open-loop simulations.
yp = mpcsim(plant,model,Kmpc,tend,r) [yp,u,ym] = mpcsim(plant,model,Kmpc,tend,r,usat,. tfilter,dplant,dmodel,dstep)
mpcsim provides a convenient way to simulate the performance of the type of system shown in the above diagram. Measurement noise can be simulated by treating it as an unmeasured disturbance. The required input variables are as follows: plant
Is a model in the MPC step format that is to represent the plant.
Is a model in the MPC step format that is to be used for state estimation in the controller. In general, it can be different from plant if you want to simulate the effect of plant/controller model mismatch. Note, however, that model should be the same as that used to calculate Kmpc.

You may use a different number of rows in the matrices r and usat, should that be appropriate for your simulation.
Let us now demonstrate the use of the controller nlmpcsim. Since the plant used in Example 1 is linear, using mpcsim would be much faster. The point, however, is to show how masked variables are specified for nlmpcsim.
1 The plant is linear with two inputs and two outputs. It is represented by
1.dx 0.2 u + 50 ------ = 1 x+ dt 0 ------1.5 y = x The Simulink S-function for this plant is in mpcplant.m. The nominal steady-state operating condition is y0 = [58.3 1.5] and u0 = [100 1]. The Simulink block to simulate this plant using nlmpcsim is in nlmpcdm1.m and shown in Figure 4-1.
nlmpcsim Demo #1 (Double click on the "?" for more info)

? Load Data Plot Results

t Clock u Save Inputs mpcplant Sfunction y1, y2 Save Time

y Save Outputs

nlmpcsim nlmpcsim
Figure 4-1 Simulink Block for Example 1
The following statements build the step response model and specify the parameter values. Note that model does not equal the plant model stored in
mpcplant.m. The important thing to notice is that both r and usat are actual variables. They are not deviation variables. g11=poly(0.4,[1 2]); g21=poly2tfd(0,1); g12=poly2tfd(0,1); g22=poly2tfd(1,[1 1]); tfinal=8; delt=0.2; nout=2; model=tfd2step(tfinal,delt,nout,g11,g21,g12,g22); ywt=[1 1]; uwt=[0 0]; M=4; P=10; r=[68.3 2]; usat=[200 200]; tfilter=[ ]; Kmpc = mpccon(model,ywt,uwt,M,P); dmodel = [ ];
There are two ways to simulate the closed loop system. We can set the simulation parameters and click on Start under Simulation or via the following statements.
plant= nlmpcdm1; y0=[58.3 1.5]; u0=[100 1]; tfsim = 2; tol=[1e-3]; minstep=[ ]; maxstep=[ ]; [t,yu]=gear(plant,tfsim,[y0 u0],[tol,minstep,maxstep]);
Figure 4-2 shows the response for the setpoint change.
Figure 4-2 Output responses for a setpoint change for Example 1

estimator), nu integrators that operate on the u signal produced by the standard MPC formulation to yield a u signal that can be used as input to the plant and as a closed-loop output, and nd differencing elements that operate on the d signal to produce the d signal required in the standard MPC formulation. If there are no measured disturbances, these states are omitted.
The closed-loop input and output variables are: r( k) z(k ) u cl ( k ) = w u ( k ) d(k ) w( k)

yp( k ) u(k) (k k) y

where y = ( k k ) is the estimate of the noise-free plant output at sampling period k based on information available at period k. This estimate is generated by the controller element. Note that ucl will include d and/or w automatically whenever pmod includes measured disturbances and/or unmeasured disturbances. Thus the length of the ucl vector will depend on the inputs you have defined in pmod and imod. Similarly, ycl depends on the number of outputs and manipulated variables. Let m and p be the lengths of ucl and ycl, respectively. Then m = 2ny + nu + nd + nw p = 2ny + nu The state-space form of the controller model, cmod, can be written as: xc(k + 1) = cxc(k) + cluc(k) yc(k) = Ccxc(k) + D cuc(k) where r( k) y (k ) d(k )

uc (k ) =

yc( k ) = u( k )
Consider the linear system: 18.9e 12.8e ----------------------- ------------------------- u ( s ) y1(s ) 16.7s + 1 21.0s + = 7s 3s u ( s ) y2(s ) 2 19.4e 6.6e ----------------------- ------------------------10.9s + 1 14.4s + 1 We build this model using the MPC Toolbox functions poly2tfd and tfd2mod.
g11=poly2tfd(12.8,[16.7 1],0,1); g21=poly2tfd(6.6,[10.9 1],0,7); g12=poly2tfd(-18.9,[21.0 1],0,3); g22=poly2tfd(-19.4,[14.4 1],0,3); delt=3; ny=2; imod=tfd2mod(delt,ny,g11,g21,g12,g22); pmod=imod; % No plant/model mismatch
P=6; % Prediction horizon. M=2; % Number of moves (input horizon). ywt=[ ]; % Output weights (default - unity on % all outputs). uwt=[ ]; % Man. Var weights (default - zero on % all man. vars). Ks=smpccon(imod,ywt,uwt,M,P);
Now we can calculate the model of the closed-loop system and check its poles for stability:
clmod=smpccl(pmod,imod,Ks); maxpole=max(abs(smpcpole(clmod)))

maxpole = 0.8869

Since this is less than 1, the plant and controller combination will be closed-loop stable. (The closed-loop system has 20 states in this example).

The following input variables are optional. In general, setting one of them equal to an empty matrix causes smpcsim to use the default value, which is given in the description.
Is a matrix giving the saturation limits on the manipulated variables. Its format is as follows: u min, 1 ( 1 ) ,
u min, n u ( 1 ) u min, n ( 2 )

u min, 1 ( 2 ) ,

u min, 1 ( N ) ( ) u max, 1 ( 1 ) u max, 1 ( 2 )

, ( )

u max, 1 ( N ) u m ax, 1 ( 1 ) u m ax, 1 ( 2 ) u max, 1 ( N )
Note that it contains three matrices of N rows. N may be different than that for the setpoint matrix, r, but the idea is the same: the saturation limits will vary for the first N sampling periods of the simulation, then be held constant at the values given in the last row of usat for the remaining periods (if any). The first matrix specifies the lower bounds on the nu manipulated variables. For example, umin,j(k) is the lower bound for manipulated variable j at time t = kT in the simulation. If umin,j(k) = inf, manipulated variable j will have no lower bound at t = kT. The second matrix gives the upper bounds on the manipulated variables. If umax,j(k) = inf, manipulated variable j will have no upper bound at t = kT.

u m ax, n ( N )

u ) (N) , ( m ax, nu, ( )
The lower and upper bounds may be either positive or negative (or zero) as long as umin,j(k) umax,j(k). The third matrix gives the limits on the rate of change of the manipulated variables. In other words, smpcsim will force|uj(k) uj(k 1)| umax,j(k). The limits on the rate of change must be nonnegative. The default is no saturation constraints, i.e., all the umin values will be set to inf, and all the umax and umax values will be set to inf. Saturation constraints are enforced by simply clipping the manipulated variable moves so that they satisfy all constraints. This is a nonoptimal solution that, in general, will differ from the results you would get using the ulim variable in scmpc.
Is a matrix of measured disturbances (see above diagram). The format is the same as for r, except that the number of columns is nd rather than ny. The default is a row of nd zeros.

a matrix function F(). If the smaller dimension of F(i) is m, and if 1(i),. , m(i) are the singular values of F(i), in decreasing magnitude, then the output sigma is a matrix of singular values arranged as follows: 1 ( 1 ) 2 ( 1 ) 1 ( 2 ) 2 ( 2 ) m ( 1 ) m ( 2 )

sigma =

1 ( N ) 2( N )
The output omega is a column vector containing the frequencies 1 ,. , N.

m ( N )

See mod2frsp, varying format for an example of the use of this function.

mod2frsp

tfd2mod, tf format

4tfd2mod, tf format

tfd2mod converts a transfer function (continuous or discrete) from the MPC tf
format into the MPC mod format, converting to discrete time if necessary.
model = tfd2mod(delt2,ny,g1,g2,g3,.,g25)
Consider a transfer function such as b 0 s + b1 s + + bn G ( s ) = --------------------------------------------------------------n n1 a0 s + a 1 s + + an or b0 + b 1 z + + bn z G ( z ) = ------------------------------------------------------------1 n a0 + a1 z + + an z The MPC tf format is a matrix consisting of three rows: row 1 row 2 row 3 The n coefficients of the numerator polynomial, b0 to bn. The n coefficients of the denominator polynomial, a0 to an. column 1: The sampling period. This must be zero if the coefficients in the above rows are for a continuous system. It must be positive otherwise. column 2: The time delay. For a continuous-time transfer function, it is in time units. For a discrete-time transfer function, it is the integer number of sampling periods of time delay. The tf matrix will always have at least two columns, since that is the minimum width of the third row. The input arguments for tfd2mod are:
The sampling period for the system. If any of the transfer functions g1,., gN are continuous-time or discrete-time with sampling period not equal to delt2, tfd2mod will convert them to discrete-time with this sampling period.
The number of output variables in the plant you are modeling.

g1, g2,.gN

A sequence of N transfer functions in the tf format described above, where N 1. These are assumed to be the individual elements of a transfer-function matrix: g 1, 1 g 2, 1 g 1, 2 g 2, 2

g 1, n g 2, n

g n y, 1 g n y , 2
Thus it should be clear that N must be an integer multiple (nu) of the number of outputs, ny. Also, tfd2mod assumes that you are supplying the transfer functions in a column-wise order. In other words, you should first give the ny transfer functions for input 1 (g1,1 to gny, 1), then the ny transfer functions for input 2 (g1,2 to gny, 2), etc.

tfd2mod converts the transfer functions to discrete-time, if necessary, and
combines them to form the output variable, model, which is a composite system in the MPC mod form.

g ny, nu

Consider the linear system: 18.9e 12.8e 3.8e ----------------------- ------------------------- u ( s ) ----------------------16.7s + 1 21.0s + 1 14.9s + 1 w ( s ) 1 = + 7s 3s u ( s ) 3s y2( s ) 2 19.4e 6.6e 4.9e ----------------------- ----------------------------------------------10.9s + 1 14.4s + 1 13.2s + 1

y1( s )

The following commands build separate models of the response to the manipulated variables, u, and the unmeasured disturbance, w, all for a sampling period T = 3 then combines them using addumd to get a model of the entire system (the pmod variable):
g11=poly2tfd(12.8,[16.7 1],0,1); g21=poly2tfd(6.6,[10.9 1],0,7); g12=poly2tfd(-18.9,[21.0 1],0,3); g22=poly2tfd(-19.4,[14.4 1],0,3); delt=3; ny=2; umod=tfd2mod(delt,ny,g11,g21,g12,g22); gw1=poly2tfd(3.8,[14.9 1],0,8); gw2=poly2tfd(4.9,[13.2 1],0,3); wmod=tfd2mod(delt,ny,gw1,gw2); pmod=addumd(umod,wmod);
The current limit on the number of input transfer functions is N = 25.

mod, poly2tfd, tfd2step

tfd2step

4tfd2step

Calculates the MIMO step response of a model in the MPC tf format. The resulting step response is in the MPC step format.
plant = tfd2step(tfinal,delt2,nout,g1) plant = tfd2step(tfinal,delt2,nout,g1,.,g25)
The input variables are as follows:

tfinal

Truncation time for step response.
Desired sampling period for step response.
Output stability indicator. For stable systems, this argument is set equal to the number of outputs, ny. For systems with one or more integrating outputs, this argument is a column vector of length ny with nout(i)=0 indicating an integrating output and nout(i)=1 indicating a stable output.
A sequence of N transfer functions in the tf format (see tf format section), where N 1. These are assumed to be the individual elements of a transfer-function matrix: g 1, 1 g 2, 1 g 1, 2 g 2, 2

gn , 1 gn , 2 y y

Thus it should be clear that N must be an integer multiple (nu) of the number of outputs, ny.
tfd2step assumes that you are supplying the transfer functions in a
column-wise order. In other words, you should first give the ny transfer functions for input 1 (g1,1 to gny, 1), then the ny transfer functions for input 2 (g1,2 to gny, 2), etc.

gn , n

The output variable plant is the calculated step response of the ny outputs with respect to all inputs. The format is as described in the step section.

Consider the linear system: 18.9e 12.8e ----------------------- ------------------------- u ( s ) y1 ( s ) 16.7s + 1 21.0s + = 7s 3s u ( s ) y2 ( s ) 2 19.4e 6.6e ----------------------- ------------------------10.9s + 1 14.4s + 1 which is the same as that considered in the mpcsim example. We build the individual tf format models, then calculate and plot the MIMO step response.
g11=poly2tfd(12.8,[16.7 1],0,1); g21=poly2tfd(6.6,[10.9 1],0,7); g12=poly2tfd(-18.9,[21.0 1],0,3); g22=poly2tfd(-19.4,[14.4 1],0,3); delt=3; ny=2; tfinal=90; plant=tfd2step(tfinal,delt,ny,g11,g21,g12,g22,gw1,gw2); plotstep(plant)
The plots should match the example output in the plotstep description.
mod2step, plotstep, ss2step

th2mod, theta format

Converts a SISO or MISO model from the theta format (as used in the System Identification Toolbox) to one in the MPC mod format. Can also combine such models to form a MIMO system.
umod = th2mod(th) [umod,emod] = th2mod(th1,th2,.,thN)

4th2mod, theta format

The System Identification Toolbox allows you to identify single-input, single-output (SISO) and multi-input, single-output (MISO) transfer functions from data. The MISO form relating an output, y, to m inputs, u1 to um, and a noise input, e, is: B2 ( z ) Bm( z) B1 ( z ) C( z) -e A ( z )y ( k ) = -------------- u 1 ( k ) + -------------- u 2 ( k ) + + ----------------u m ( k ) + ----------- ( k ) D(z) F1 ( z ) F2 ( z ) Fm ( z ) where A, Bi, C, D, and Fi are polynomials in the forward-shift operator, z. The System Identification Toolbox automatically stores such models in a special format, the theta format. See the System Identification Toolbox Users Guide for details.
th2mod converts one or more MISO theta models into the MPC mod format,
which you can then use with the MPC Toolbox functions. If you supply a single input argument, th, and a single output argument, umod, then umod will model the response of a single output, y, to m inputs, u1 to um, where m 1. The value of m depends on the number of inputs included in the input model, th. Note that umod will reflect the values of the A(z), B(z), and F(z) polynomials in eq. 1. If you supply a second output argument, emod, it will model the response of y to the noise, e, i.e., the A(z), C(z) and D(z) polynomials in eq. 1. If you supply p input models (1 p 8), tfd2mod assumes that they define a MIMO system in the following form: B 11 ( z ) B 1m ( z ) C1 ( z ) A 1 ( z )y 1 ( k ) = -----------------u 1 ( k ) + + -------------------u m ( k ) + -------------- e 1 ( k ) F 11 ( z ) F 1m ( z ) D1( z )
B p1 ( z ) B pm ( z ) Cp( z ) A p ( z )y p ( k ) = ----------------- u 1 ( k ) + + ------------------- u m ( k ) + -------------- e p ( k ) F p1 ( z ) F pm ( z ) Dp (z )

Model Predictive Control Problem Setup
Prediction Model. Offsets. Optimization Problem. Standard Form. Alternative Cost Function. Terminal Weights and Constraints. Custom Constraints on Inputs and Outputs. State Estimation. Measurement Noise Model. Output Disturbance Model. State Observer. QP Matrices. Prediction. Optimization Variables. Cost Function. 2-2 2-4 2-5 2-5 2-8 2-8 2-11 2-12 2-12 2-13 2-13 2-16 2-16 2-17 2-19
Constraints. Model Predictive Control Computation. Unconstrained MPC. Constrained Model Predictive Control. Using Identified Models.

2-20 2-22 2-22 2-22 2-23

Model Predictive Control Simulink Library
MPC Library. MPC Controller Block. MPC Controller Block Mask. Input Signals. Output Signals. Look Ahead and Signals from the Workspace. Initialization. Using Model Predictive Control Toolbox Software with Simulink Coder Software. Multiple MPC Controllers. 3-2 3-3 3-3 3-4 3-5 3-6 3-7 3-7 3-9

Case-Study Examples

Introduction. Servomechanism Controller. System Model. Control Objectives and Constraints. Defining the Plant Model. Controller Design Using MPCTOOL. 4-2 4-3 4-3 4-5 4-5 4-6
Using Model Predictive Control Toolbox Commands. Using MPC Tools in Simulink. Paper Machine Process Control. Linearizing the Nonlinear Model. MPC Design. Controlling the Nonlinear Plant in Simulink. Bumpless Transfer in MPC. Nonlinear Control Using Multiple Models. Custom Constraints in a Blending Process. About the Blending Process. MPC Controller with Custom Input/Output Constraints. Using the Tuning Advisor. Using Terminal Penalty to Provide LQR Performance. Bibliography.
4-20 4-24 4-27 4-28 4-29 4-36 4-40 4-47 4-55 4-55 4-56 4-66

4-71 4-77

Reference for the Design Tool GUI
Opening the MPC Design Tool. Creating a New MPC Design Task. Menu Bar. File Menu. MPC Menu. Toolbar. 5-3 5-5 5-6 5-6 5-7 5-9
Tree View. Node Types. Renaming a Node. Importing a Plant Model. Import from. Import to. Buttons. Importing a Linearized Plant Model. Importing a Controller. Import from. Import to. Buttons. Exporting a Controller. Dialog Box Options. Buttons. Signal Definition View. MPC Structure Overview. Buttons. Signal Properties Tables. Right-Click Menu Options. Plant Models View. Plant Models List. Model Details. Additional Notes. Buttons. Right-Click Options. Controllers View. Controllers List. Controller Details. Additional Notes. Buttons. Right-Click Options. Simulation Scenarios List. Scenarios List.
5-10 5-10 5-10 5-12 5-13 5-14 5-14 5-15 5-18 5-18 5-20 5-20 5-21 5-21 5-22 5-23 5-24 5-24 5-24 5-26 5-28 5-29 5-29 5-30 5-30 5-30 5-31 5-32 5-32 5-33 5-33 5-34 5-35 5-36
Scenario Details. Additional Notes. Buttons. Right-Click Options. Controller Specifications View. Model and Horizons Tab. Constraints Tab. Constraint Softening. Weight Tuning Tab. Entering Vectors in Table Cells. Estimation Tab. Right-Click Menus. Simulation Scenario View. Simulation Settings. Setpoints. Measured Disturbances. Unmeasured Disturbances. Signal Type Settings. Simulation Button. Tuning Advisor Button. Right-Click Menus. Tuning Advisor. Defining the Performance Metric. Baseline Performance. Sensitivities and Tuning Advice. Updating the Controller. Restoring Baseline Tuning. Modal Dialog Behavior. Scenarios for Performance Measurement. Response Plots. Data Markers. Displaying Multiple Scenarios. Viewing Selected Variables. Grouping Variables in a Single Plot. Normalizing Response Amplitudes.

compensate for vs impact on y immediately rather than waiting until the effect appears in the y measurement. This approach is called feedforward control. Model Predictive Control Toolbox design always provides feedback compensation for unmeasured disturbances and feedforward compensation for any measured disturbance. The design requires a model of the impact that v and u have on y (symbolically, v y and u y ). It uses this plant model to calculate the u adjustments needed to keep y at its setpoint. This calculation considers the effect of any known constraints on the adjustments (typically an actuator upper or lower bound or a constraint on how rapidly u can vary). You can also specify bounds on y. The ability to create such constraint specifications is a distinguishing feature of Model Predictive Control Toolbox design and can be particularly valuable when you have multiple control objectives to be achieved via multiple adjustments (a MIMO plant). In the context of a SISO system, such constraint handling is often termed as anti-windup. If the plant model is accurate, the plant responds quickly to adjustments in u, and no constraints are encountered, feedforward compensation can counteract the impact of v perfectly. In reality, model imperfections, physical limitations, and unmeasured disturbances cause the y to deviate from its setpoint. Therefore, Model Predictive Control Toolbox design includes a disturbance model ( d y ) to estimate d and predict its impact on y. It then uses its u y model to calculate appropriate adjustments (feedback). This calculation also considers the known constraints. Various noise effects can corrupt the measurement. The signal z in the Block Diagram of a SISO Model Predictive Control Toolbox Application on page 1-2 figure represents such effects. Theses effects can vary randomly with a zero mean, or can exhibit a nonzero, drifting bias. Model Predictive Control Toolbox design uses a z y model in combination with its d y model to remove the estimated noise component (filtering).
The preceding feedforward/feedback actions comprise the controllers regulator mode. The design also provides a servo mode, i.e., it adjusts u such that y tracks a time-varying setpoint. The tracking accuracy depends on the plant characteristics (including constraints), the accuracy of the u y model, and whether or not future setpoint variations can be anticipated, i.e., known in advance. If so, it provides feedforward compensation for these.

Typical Sampling Instant

Model Predictive Control Toolbox design generates a discrete-time controllerone that takes action at regularly spaced, discrete time instants. The sampling instants are the times at which the controller acts. The interval separating successive sampling instants is the sampling period, t (also called the control interval). This section provides more details on the events occuring at each sampling instant.

Alternative Cost Function
You have the option to use the following quadratic objective instead of the standard one (Equation 2-3):

J (u, ) =

p1 i =0
y ( k + i + 1 k ) r ( k + i + 1)
Q y ( k + i + 1 k ) r ( k + i + 1) + u ( k + i k ) Ru ( k + i k )
+ u ( k + i k ) utarget ( k + i )
Ru u ( k + i k ) utarget ( k + i ) + e 2
where Q is an ny by ny matrix, and Ru and Ru are nu by nu matrices, all positive semidefinite. Equation 2-5 allows nonzero off-diagonal weights but uses the same weights at each step in the prediction horizon. Equation 2-3 and Equation 2-5 are equivalent when:
u Weights wiy j , wi, u , and wi, j are constant for all i = 1,., p. , j
Matrices Q, Ru and Ru are diagonal with the squares of the weights wiy j , ,
u wi, u , and wi, j respectively as their diagonal elements. j
Note When using the alternative cost function, you must define the controller using MATLAB commands. The Model Predictive Control Toolbox design tool does not provide this option.
Terminal Weights and Constraints
Terminal weights are the quadratic weights Wy on y(t+p) and Wu on u(t + p 1), where p is the prediction horizon, and apply at time k +p only, i.e., the prediction horizons final step. Using terminal weights, you can achieve infinite horizon control that guarantees closed-loop stability. However, before using terminal weights, you must distinguish between problems with and without constraints. Terminal constraints are the constraints on y(t + p) and u(t + p 1), where p is the prediction horizon. You can use terminal constraints as an alternative way to achieve closed-loop stability by defining a terminal region.
Note Terminal weights and constraints can be used only at the command-line. See setterminal. For the relatively simple unconstrained case, a terminal weight can make the finite-horizon Model Predictive Controller behave as if its prediction horizon were infinite. For example, the MPC controller behavior is identical to a linear-quadratic regulator (LQR). The standard LQR derives from the cost function:
J (u) x(k i)T Qx(k i) u(k i 1)T Ru(k i 1)
where x is the vector of plant states in the standard state-space form: x(k + 1) = Ax + Bu(k) (2-7)
The LQR provides guaranteed nominal stability provided matrices Q and R meet certain conditions. The LQR can be converted to a finite-horizon form as follows:

J (u)

[ x(k i)T Qx(k i) u(k i 1)T Ru(k i 1)] x(k p)T Qp x(k p)
where Qp , the terminal penalty matrix, is the solution of the Riccati equation:
Q p AT Q p A AT Q p B( BT Q p B R)1 BT Q p A Q
which you can obtain using the lqr command in Control System Toolbox software. In general, Qp is a full (symmetric) matrix. You cannot use the standard Model Predictive Control Toolbox cost function to implement the LQR cost function, except only for the first p 1 steps if Q and R are diagonal matrices. Also, you cannot use the alternative cost function because it employs identical weights at each step in the horizon, and by definition, the terminal weight differs from those in steps 1 to p 1. Instead, use the following steps:

y y y y y y Wy = diag w1,1 , w1,2 ,., w1,n ,., w p,1 , w p,2 ,., w p,n y y

(2-12)

Finally, after substituting u(k), u(k), y(k), J(z) can be rewritten as
T T r(1) T utarget (0) v(0) K + x(0)T K z K + K + u(1)T K + J ( z, ) = + z K u z + 2 r x v u ut r( p) utarget ( p 1) v( p) (2-13) + constant 2 T
Note In order to keep the QP problem strictly convex, if the condition number of the Hessian matrix KU is larger than 1012, the quantity 10*sqrt(eps) is added on each diagonal term. This can only occur when all input rates are unpenalized (Wu=0) (see Weights in the Model Predictive Control Toolbox Reference).
If the alternative cost function shown in Equation 2-5 is being used, Equation 2-12 is replaced by the following:
Wu = blkdiag ( Ru ,., Ru ) Wy = blkdiag ( Q,., Q )
Wu = blkdiag ( Ru ,., Ru )

(2-14)

where the block-diagonal matrices repeat p times, i.e., once for each step in the prediction horizon. You also have the option to use a combination of the standard and alternative forms. See Weights in the Model Predictive Control Toolbox Reference for more details.
Now consider the limits on inputs, input increments, and outputs along with the constraint 0.
y y ymin (1) Vmin (1) ymax (1) + Vmax (1) y(1) y y ymin ( p) Vmin ( p) ymax ( p) + Vmax ( p) y( p) u u umin (0) Vmin (0) umax (0) + Vmax (0) u(0) u umin ( p 1) V u ( p 1) u( p 1) umax ( p 1) + Vmax ( p 1) min u(0) u u umin (0) Vmin (0) umax (0) + Vmax (0) u( p 1) u u u umax ( p 1) + Vmax ( p 1) min ( p 1) Vmin ( p 1)
Note Upper and lower bounds that are not finite are removed, as well as the input and input-increment bounds over blocked moves.
Similarly to what was done for the cost function, we can substitute u(k), u(k), y(k), and obtain

M z z + M Mlim

v(0) + Mv + Mu u(1) + M x x(0) v( p)

(2-15)

where matrices Mz,M,Mlim,Mv,Mu,Mx are obtained from the upper and lower bounds and ECR values.
Model Predictive Control Computation
This section describes how the model predictive control optimization problem is solved at each time step k (in mpcmove, mpc_sfun.mex, and mpcloop_engine.mex) by using the matrices built at initialization described in QP Matrices on page 2-16.

Unconstrained MPC

The optimal solution is computed analytically
T r(1) T utarget (0) v(0) K + x(0)T K 1 z* = K u K r + K v + u(1)T K u + ut x r( p) utarget ( p 1) v( p) T
and the model predictive controller sets u(k)=z*0, u(k)=u(k1)+u(k).
Constrained Model Predictive Control
The optimal solution z*, * is computed by solving the quadratic program described in Equation 2-13 and Equation 2-15, using the QP solver coded in the qpsolver.mex function (see qpdantz for more details).

Using Identified Models

Model Predictive Control Toolbox software is able to employ plant models generated by System Identification Toolbox software from input/output measurements. Model Predictive Control Toolbox software labels control input signals as Manipulated, measured input disturbances as Measured, and unmeasured input disturbances as Unmeasured. On the other hand, System Identification Toolbox software has a different naming rule, as it calls Measured the inputs that are measurable quantities, and Noise those that are not. When you specify an identified model in the Model Predictive Control constructor as the plant model, Model Predictive Control Toolbox software treats Noise signals as Unmeasured input signals, and Measured signals as Manipulated signals, assuming that all measured inputs are also manipulated variables. You can later change signal types, for instance, to specify that some measured inputs are measured disturbances, rather than manipulated variables (see setmpcsignals). You can also change their names, which affects displays in graphics and tables (see setname). Model Predictive Control Toolbox software internally converts the identified model you have provided as a plant model into the classical (A,B,C,D) state-space format. The columns of the B matrix originally related to Noise channels are treated as the effect of unmeasured input disturbances on the state of the plant. On the other hand, the columns of the D matrix related to Noise channels as treated as the effect of measurement noise superimposed on the output signal. Accordingly, Model Predictive Control Toolbox software treats as the plant model the state-space model obtained from (A,B,C,D) by zeroing the columns of D related to Noise channels. Those columns are instead used as a static noise model, or cascaded to an existing noise model if you have specified one. A unit static gain is assumed as the disturbance model, unless you have specified another one.
MPC Library on page 3-2 MPC Controller Block on page 3-3 Multiple MPC Controllers on page 3-9

MPC Library

The MPC Simulink Library provides two blocks you can use to model MPC control in Simulink. Access the library using the Simulink Library Browser or by typing mpclib at the command prompt. The latter reveals the librarys contents as shown in the following figure.

being used, the controller sampling period, and the prediction and control horizons. All are editable. For now, leave them at their default values. The buttons on the Controllers view allow you to: Import a controller designed previously and stored either in your workspace or in a MAT-file. Export the selected controller to your workspace. Create a New controller, which will be initialized to the Model Predictive Control Toolbox defaults. Copy the selected controller to create a duplicate that you can modify. Delete the selected controller.
Specifying Controller Properties
Select the MPC1 subnode. The main pane should change to the controller design. If the selected Prediction model is continuous-time, as in this example, the Control interval (sampling period) defaults to 1. You need to change this to an application-appropriate value. Set it to 0.1 seconds (as shown in Controller Design View, Models and Horizons Pane on page 4-11). Leave the other values at their defaults for now.
Controller Design View, Models and Horizons Pane

Specifying Constraints

Next, click the Constraints tab. The view shown in Controller Design View, Constraints Pane on page 4-11 appears. Enter the appropriate constraint values. Leaving a field blank implies that there is no constraint.
Controller Design View, Constraints Pane
In general, its good practice to include all known manipulated variable constraints, but its unwise to enter constraints on outputs unless they are an essential aspect of your application. The limit on applied torque is such a constraint, as are the limits on applied voltage. The angular position has physical limits but the controller shouldnt attempt to enforce them, so you should leave the corresponding fields blank (see Controller Design View, Constraints Pane on page 4-11). The Max down rate should be nonpositive (or blank). It limits the amount a manipulated variable can decrease in a single control interval. Similarly, the Max up rate should be nonnegative. It limits the increasing rate. Leave both unconstrained (i.e., blank). The shaded columns cant be edited. If you want to change this descriptive information, select the root node view and edit its tables. Such changes apply to all controllers in the design.

u = G p

where Gp is the flow rate of stock entering the feed tank, and Gw is the recycled white water flow rate. The consistency of stock entering the feed tank, Np, is a measured disturbance. v = Np
The white water consistency is an unmeasured disturbance. d = Nw Variables are normalized. All are zero at the nominal steady state and have comparable numerical ranges. Time units are minutes. The process is open-loop stable. The mpcdemos folder contains the file mpc_pmmodel.m, which implements the nonlinear model equations as a Simulink S-function. The input sequence is Gp, Gw, Np, Nw, and the output sequence is 2, N1, N2.
Linearizing the Nonlinear Model
The paper machine headbox model is easy to linearize analytically, yielding the following state space matrices:
A = [-1.0.3940 -0.0 -0.0.8200 -0.7840 0.4130 -0.4260]; B = [1.2740 1.0 1.3400 -0.6500 0.2030 0.0]; C = [0 1.0 1.0 1.0000]; D = zeros(3,4);
Use these to create a continuous-time LTI state-space model, as follows:
PaperMach = ss(A, B, C, D); PaperMach.InputName = {'G_p', 'G_w', 'N_p', 'N_w'}; PaperMach.OutputName = {'H_2', 'N_1', 'N_2'};
(The last two commands are optional; they improve plot labeling.) As a quick check of model validity, plot its step responses as follows:

step(PaperMach);

The results appear in the following figure. Note the following: The two manipulated variables affect all three outputs. They have nearly identical effects on H2. The GwN2 pairing exhibits an inverse response. These features make it difficult to achieve accurate, independent control of H2 and N2.
Linearized Paper Machine Models Step Responses

MPC Design

to open the MPC design tool. Import your LTI Paper Mach model as described in Opening MPCTOOL and Importing a Model on page 4-7.
Next, define signal properties, being sure to designate Np and Nw as measured and unmeasured disturbances, respectively. Your specifications should resemble Signal Properties for the Paper Machine Application on page 4-30.

Unmeasured Input Disturbance
To test the response to unmeasured disturbances, define another new scenario called Feedback. Configure it as for Feedforward, but set the measured disturbance, Np, to zero (constant), and the unmeasured disturbance, Nw, to 1.0 (constant). This simulates a sudden, sustained, unmeasured disturbance occurring at time zero. Running the simulation should yield results like those in Feedback Scenario: Unmeasured Disturbance Rejection on page 4-34. The two controlled outputs
(H2 and N2 ) exhibit relatively small deviations from their setpoints (which are zero). The settling time is longer than for the servo response (compare to Servo Response for Unit Step in Headbox Level Setpoint on page 4-31) which is typical.
Feedback Scenario: Unmeasured Disturbance Rejection
One factor limiting performance is the chosen control interval of 2 time units. The controller cant respond to the disturbance until it first appears in the outputs, i.e., at t = 2. If you wish, experiment with larger and smaller intervals (modify the specification on the controllers Model and Horizons tab).
Effect of Estimator Assumptions
Another factor influencing the response to unmeasured disturbances (and model prediction error) is the estimator configuration. The results shown in Feedback Scenario: Unmeasured Disturbance Rejection on page 4-34 are for the default configuration. To view the default assumptions, select the controller node (MPC1), and click its Estimation tab. The resulting view should be as shown in Default Estimator Assumptions: Output Disturbances on page 4-35. The status
message (bottom of figure) indicates that Model Predictive Control Toolbox default assumptions are being used.
Default Estimator Assumptions: Output Disturbances
Now consider the upper part of the figure. The Output Disturbances tab is active, and its Signal-by-signal option is selected. According to the tabular data, the controller is assuming independent, step-like disturbances (i.e., integrated white noise) in the first two outputs. Click the Input Disturbances tab. Verify that the controller is also assuming independent step-like disturbances in the unmeasured disturbance input.
Thus, there are a total of three independent, sustained (step-like) disturbances. This allows the controller to eliminate offset in all three measured outputs. The disturbance magnitudes are unity by default. Making one larger than the rest would signify a more important disturbance at that location. Click the Measurement Noise tab. Verify that white noise (unit magnitude) is being added to each output. The noise magnitude governs how much influence each measurement has on the controllers decisions. For example, if a particular measurement is relatively noisy, the controller will give it less weight, relying instead upon the model predictions of that output. This provides a noise filtering capability. In the paper machine application, the default disturbance assumptions are reasonable. It is difficult to improve disturbance rejection significantly by modifying them.

y(k i)T Qy(k i) u(k i 1)T Ru(k i 1) x(k p)T Qp x(k p)
where Qp is a penalty applied at the last (i.e., terminal) prediction horizon step, and the prediction and control horizons are equal, i.e., p = m. The required Qp is the Ricatti matrix that you can calculate using the Control System Toolbox lqr and lqry commands. The value is a positive definite symmetric matrix. The following procedure demonstrates how to design an unconstrained MPC controller that provides performance equivalent to a LQR controller:
1 Define a plant with one input and two outputs.
The plant is a double-integrator, represented as a state-space model in discrete-time with sampling interval 0.1 seconds.
A = [1 0;0.1 1]; B = [0.1; 0.005]; C = eye(2); D = zeros(2,1); Ts = 0.1; Plant = ss(A, B, C, D, Ts); Plant.InputName = {'u'}; Plant.OutputName = {'x_1', 'x_2'};
2 Design a LQR controller with output feedback for the plant.
Q = eye(2); R = 1; [K, Qp] = lqry(Plant, Q, R); Q and R are output and input weight matrices, respectively. Qp is the Ricatti matrix.
3 Design an MPC controller equivalent to the LQR controller.
To implement Equation 4-2, compute L, the Cholesky decompostion of Qp, such that LTL= Qp. Then, define auxiliary unmeasured output variables ya(k) = Lx(k) such that yaTya = xTQpx. For the first p-1 prediction horizon steps, the standard Q and R weights apply to the original u and y, and ya has a zero penalty. On step p, the original u and y have zero penalties, and ya has a unity penalty.
a Augment the plant model, and specify the augmented outputs as

unmeasured.

NewPlant = Plant; cholP = chol(Qp); set(NewPlant,'C',[C;cholP],'D',[D;zeros(2,1)],'OutputName', {'x_1', 'x_2', 'Cx_1', 'Cx_2'}); NewPlant.InputGroup.MV = 1; NewPlant.OutputGroup.MO = [1 2]; NewPlant.OutputGroup.UO = [3 4];
b Create an MPC controller with prediction horizon equals control horizon.
P = 3; M = 3; MPCobj = mpc(NewPlant, Ts, P, M);
When there are no constraints, you can use a rather short horizon (in this case, p 1 gives identical results).
c Specify weights for manipulated variable (MV) and output variable (OV).
ywt = sqrt(diag(Q))'; uwt = sqrt(diag(R))'; MPCobj.Weights.OV = [ywt 0 0]; MPCobj.Weights.MV = uwt; MPCobj.Weights.MVrate = 1e-6;
The two augmented outputs have zero weights during the prediction horizon.
d Specify terminal weights.
To obtain the desired effect, define unity weights for these at the final point in the horizon.
U = struct('Weight', uwt); Y = struct('Weight', [1 1]); setterminal(MPCobj, Y, U); setoutdist(MPCobj,'remove'); % Remove added state estimator setestim(MPCobj, C); % State estimates = measured values

1 Configures the Control and Estimation Tools Manager. 2 Temporarily inserts linearization input and output points in the Simulink
model at the inputs and outputs of the MPC block.
3 When the Create a new operating condition from MPC I/O values
option is selected, the Model Predictive Control Toolbox software temporarily inserts output constraints at the inputs/outputs of the MPC block.
4 Finds a steady state operating condition based on the constraints or uses
the specified operating condition.
5 Linearizes the plant model about the operating point.
The linearized plant model appears as a new node under Plant Models. For details of the linearization process, refer to the Simulink Control Design documentation.

Linearization Options

You can customize the linearization process in several ways: To specify a name for the linearized plant model, enter the name in the Linearization model name edit box. To use an alternative operating condition, you can:
Select one from the menu next to Use the previously computed operating condition. This list contains all operating conditions that exist within the current project. Select Create a new operating condition from MPC I/O values to compute an operating condition by optimization, using the nominal plant values as constraints. See the Getting Started for an example involving a nonlinear chemical reactor (CSTR).
To replace the nominal plant values with the operating point used in the linearization, select the check box next to Replace the MPC nominal I/O values with those derived from the operating condition. When there are multiple MPC blocks in the Simulink diagram, use the Import to menu to select the one that will receive the plant model. Note The above linearization process automatically identifies the plants input and output variables according your signal connections in the Simulink model. The controller block does not allow signals corresponding to unmeasured disturbance or unmeasured output variables. Consequently, such variables cannot be included in a model created via the above linearization procedure. If you must include such variables in your controller design, use the Simulink Control Design tool to designate the signals to be used, linearize the plant, and then import this linearized model into the MPC design tool. See the Getting Started for an example of this procedure.

Importing a Controller

To import a controller, do one of the following: Select the MPC/Import/Controller menu option. Select the MPC project/task node in the tree (see Tree View on page 5-10), and then click the Import Controller button. Right-click the MPC project/task node and select the Import Controller menu option. If youve already designed a controller, select the Controllers node and then click the Import button, or right-click the Controllers node and select the Import Controller menu option. The MPC Controller Importer dialog box opens. The following sections describe its options: Import from on page 5-18 Import to on page 5-20 Buttons on page 5-20

The Name field is editable. Each model must have a unique name. The name you assign here will be used within the design tool only. The Type field is noneditable and indicates the models LTI object type (see the Control System Toolbox documentation for a detailed discussion of LTI models). The Sampling Period field is zero for continuous-time models, and a positive real value for discrete-time models. The Imported on field gives the date and time the model was imported.

Model Details

This scrollable viewport shows details of the model currently selected in the plant models list (see Plant Models List on page 5-29). An example appears below.

Additional Notes

You can use this editable text area to enter comments, distinguishing model features, etc.
Opens the Plant Model Importer dialog box (see Importing a Plant Model on page 5-12).

Delete

Deletes the selected model. If the model is being used elsewhere (i.e., in a controller or scenario), the first model in the list replaces it and a warning message appears.

Right-Click Options

Right-clicking the Plant models node causes the following menu option to appear.

Import Model

Controllers View
Selecting Controllers in the tree displays this view (see example below).
The following sections describe the views main features: Controllers List on page 5-32 Controller Details on page 5-32 Additional Notes on page 5-33 Buttons on page 5-33 Right-Click Options on page 5-34

Controllers List

This table lists all the controllers in your project. The example below lists two controllers, MPC1 and MPC2.
The Name field is editable. The name you assign here must be unique. You will refer to it elsewhere in the design tool, e.g., when you use the controller in a simulation scenario. Each listed controller corresponds to a subnode of Controllers (see Tree View on page 5-10). Editing the name in the table will rename the corresponding subnode. The Plant Model field is editable. To change the selection, click the cell and choose one of your models from the list. (All models appearing in the Plant Models view are valid choices. See Plant Models View on page 5-28.) The Control Interval field is editable and must be a positive real number. You can also set it in the Controller Specifications view (see Model and Horizons Tab on page 5-39 for more details). The Prediction Horizon field is editable and must be a positive, finite integer. You can also set in the Controller Specifications view (see Model and Horizons Tab on page 5-39 for more details). The noneditable Last Update field gives the date and time the controller was most recently modified.

Prior to Time, the signal = Initial Value. At Time, a square pulse of duration Period and magnitude Size occurs. y = y0 for 0 t < t0, where y0 = Initial Value, t0 = Time y = y0+M for t0 t t0, where M = Size, T = Period y = y0 for t t0 + T

Gaussian

Prior to Time, the signal = Initial Value. At Time, the signal begins to vary randomly about Initial Value with standard deviation Size.
y = y0 for 0 t < t0, where y0 = Initial Value, t0 = Time y = y0+Mrandn for t t0, where M = Size.
randn is the MATLAB random-normal function, which generates random
numbers having zero mean and unit variance.

Simulation Button

Click the Simulate button to simulate the scenario. You can also press Ctrl+R, use the toolbar icon (see Toolbar on page 5-9), or use the MPC/Simulate menu option (see Menu Bar on page 5-6).

Tuning Advisor Button

Click Tuning Advisor to open a window that helps you improve your controllers performance. See Tuning Advisor on page 5-69.

Copy Scenario

Creates a new simulation scenario having the same settings and a default name.

Delete Scenario

Deletes the scenario.

Rename Scenario

Opens a dialog box allowing you to rename the scenario. Note Each scenario in a design project/task must have a unique name.

Tuning Advisor

When you design MPC controllers, you can use the Tuning Advisor to help you determine which weight has the most influence on the closed-loop performance. The Tuning Advisor also helps you determine in which direction to change the weight to improve performance. Using the Advisor, you can know numerically how each weight impacts the closed-loop performance, which makes designing MPC controllers easier when the closed-loop responses does not depend intuitively on the weights. To start the Tuning Advisor, click Tuning Advisor in a simulation scenario view (see Tuning Advisor Button on page 5-68). The next figure shows the default Tuning Advisor window for a distillation process in which there are two controlled outputs, two manipulated variables, and one measured disturbance (which the Tuning Advisor ignores). In this case, the originating scenario is Scenario1.
The Tuning Advisor populates the Current Tuning column with the most recent tuning weights of the controller displayed in the Controller in Design. In this case, Obj is the controller. The Advisor also initializes the Performance Weight column to the same values. The Scenario in Design

controller 3-7 observer 2-14 observer form 2-12 state estimation 1-14 states plant 2-3
tasks importing model 5-14 tolerance bands 1-13 toolbar overview 5-9 tree view node types 5-10 overview 5-10 Tuning Advisor 5-69
units specification 5-24 unmeasured disturbances example 1-2 using identified models, system identification model conversion 2-23 Using the Tuning Advisor 4-66
views signal definition 5-23 tree 5-10
weights inputs 5-49 optimization problem 2-5 outputs 5-50

Index-4

specification dialog 5-48 tuning 5-48 workspace
controller export 5-21 importing model from 5-13

Index-5