The MPC Simulink Library User’s Guide Version 1 A. Bemporad M. Morari N. L. Ricker MPC Simulink Library How to Contact A. Bemporad: Automatic Control Laboratory Physikstrasse 3 ETH-Z, ETL I26 8092 Zurich, Switzerland http://control.ethz.ch/~bemporad/ bemporad@aut.ee.ethz.ch The MPC Simulink Library ? COPYRIGHT 2000 by A. Bemporad, M. Morari, and N. L. Ricker The software described in this document is covered by the MPC Tools license agreement. The software may be used or copied only under the terms of the license agreement. This manual may be photocopied and reproduced, but no part may be included in any other document without prior written consent from the authors. U.S. GOVERNMENT: If Licensee is acquiring the Programs on behalf of any unit or agency of the U.S. Government, the following shall apply: (a) For units of the Department of Defense: the Government shall have only the rights specified in the license under which the commercial computer software or commercial software documentation was obtained, as set forth in subparagraph (a) of the Rights in Commercial Computer Software or Commercial Software Documentation Clause at DFARS 227.7202-3, therefore the rights set forth herein shall apply; and (b) For any other unit or agency: NOTICE: Notwithstanding any other lease or license agreement that may pertain to, or accompany the delivery of, the computer software and accompanying documentation, the rights of the Government regarding its use, reproduction, and disclosure are as set forth in Clause 52.227-19 (c)(2) of the FAR. MATLAB, Simulink, Handle Graphics, and Real-Time Workshop are registered trademarks and Stateflow and Target Language Compiler are trademarks of The MathWorks, Inc. Other product or brand names are trademarks or registered trademarks of their respective holders. MPC Graphical User Interface -3 Printing History: Jan 2000 Draft MPC Simulink Library vi Preface Disclaimer State of Development The MPC Simulink Library is in a developmental (beta) stage. It is not an official MathWorks product. It has been tested extensively, but it is likely that some problems remain. Please review all calculations with a critical eye. We will also be improving features in response to your comments. In Case of Difficulty This beta version of the MPC Simulink Library is not being supported by The MathWorks. If you encounter a problem, please check the FAQ (see "Updates" below). If that doesn’t help, please send a detailed description by e-mail to A. Bemporad at the following address: bemporad@aut.ee.ethz.ch He would also be happy to have your suggestions for improvements in the usability of the MPC Simulink Library, its documentation, and the MPC Toolbox in general. Updates The FAQ and latest beta version of the MPC Simulink Library will be maintained for download at http://control.ethz./~bemporad/toolbox/mpclib.html Acknowledgements We appreciate the help of Federica Rusconi, Domenico Mignone, Carles Pedret Ferre, Adrian Toller, Konrad Stadler, Tobias Raithel, Kazuro Tsuda, Jay Lee, Pascal Gahinet, and Greg Wolodkin. 1 Introduction MPC Simulink Library Overview . . . . . . . . . . . . . . . 1-2 Options . . . . . . . . . . . . . . . . . . . . . . . . . . 1-2 Starting the MPC Simulink Library . . . . . . . . . . . . . . 1-3 System Requirements . . . . . . . . . . . . . . . . . . . . 1-3 Installation . . . . . . . . . . . . . . . . . . . . . . . . . 1-3 Quick Start . . . . . . . . . . . . . . . . . . . . . . . . 1-4 Leisurely Start . . . . . . . . . . . . . . . . . . . . . . . 1-4 MPC Fundamentals . . . . . . . . . . . . . . . . . . . . 1-5 1 Introduction 1-2 MPC Simulink Library Overview The MPC Simulink Library is designed to help you analyze and simulate Model Predictive Control (MPC) modules within any Simulink description of the environment. Options The MPC Simulink Library supports four controller blocks, to be connected in feedback with the system to regulate. Each block receives output measurements and returns the control input action to the system. The four blocks are the following: ? Regulator. The MPC block regulates the output of the system to zero. ? Controller. An additional reference signal is received by the MPC controller. The output of the system will track such a signal. ? Controller with Measured Disturbance Rejection. An additional information about measured input disturbance signals entering the system is received by the MPC controller and taken into account in the computation of the control action. ? Controller with Anticipation. Reference and measured disturbance signals are read from file. This allows knowing future samples when computing the control action. Starting the MPC Simulink Library 1-3 Starting the MPC Simulink Library System Requirements You’ll need the following MATLAB software: ? MPC Toolbox for MATLAB Version 5, including the MPC Simulink Library files (which are currently an add-on to the standard MPC Tools release). ? MATLAB Version 5.3 or greater. ? Simulink Version 3.0 or greater. ? Control Toolbox. Installation 1 Make a backup copy of your MPC Tools directory (which should be located in your MATLAB directory, e.g., MATLAB/toolbox/mpc). This will allow you to restore your original state if the GUI installation interferes in some unexpected way. 2 Download the latest version of the MPC Simulink Library from http://control.ethz.ch/~bemporad/toolbox/mpclib.html It’s a ZIP archive. Extract the contents using one of the many available utility programs (e.g., WinZip). This should create a new mpclib directory within the mpc directory. IMPORTANT: The ZIP archive is not a full version of MPC Tools. You must have existing mpc/mpccmds and mpc/mpcdemos directories. The files from the archive only update these directories. 3 Set your MATLAB Path so the new mpc/mpclib directory is included. Verify that the original mpc/mpccmds and mpc/mpcdemos directories are also on the path. 4 Type “mpclib” in the MATLAB Command window. You should see the MPC Simulink Library. 5 If it doesn’t work or you have other problems, see section In Case of Difficulty in the preface to this document. 1 Introduction 1-4 Quick Start 1 Type “mpclibdemo” in the command window. 2 Choose one of the tutorial examples. 3 Examine the Simulink MDL file and the Matlab M file for the demo. 4 Try it with your own applications! Leisurely Start The following sections of this document provide tutorial examples and additional details. If possible, work through the steps in MATLAB/Simulink as you read. MPC Fundamentals 1-5 MPC Fundamentals This section reviews MPC concepts and defines key terms. If you’re already familiar with MPC you may prefer to read Chapter 2, “MPC Worksheet Basics” next. Overview: single-input, single-output (SISO) MPC Figure 1-1 shows a situation in which MPC is trying to hold a single variable, , at a target value, r, by adjusting the “manipulated variable” (or “actuator”) u. See Table 1-1 for a brief description of the signals appearing in Figure 1-1. The box labeled “plant” is the real process or mechanism that produces . It responds to changes in u as well as to two types of disturbance signals: measured, v, and unmeasured, d. (Some applications have unmeasured disturbances only; the v signal is optional.) The unmeasured disturbance represents the usual mysterious events that upset plant operation, causing variations in . The only sign that such an event has occured is a change in the measured output, y. Figure 1-1 Block diagram of a single-input, single-output MPC application. y y y + +PlantMPC v r y d z yyu v Measured Disturbance Measured Output (Controlled Variable) Noise Setpoint Unmeasured Disturbance Actuator Plant Output 1 Introduction 1-6 The measured disturbance also affects . The MPC block includes models of the way in which v and u affect (symbolically, and ). It uses this information to calculate u adjustments that keep at its setpoint in spite of the (known) disturbance. This calculation considers the effect of any known constraints on the adjustments (e.g., an actuator at its upper or lower bound). If the models are accurate and the plant responds quickly to u, this “feedforward” compensation counteracts the effect of v perfectly. Table 1-1 Summary of MPC signals (see also Figure 1-1). Symbol Description d Unmeasured disturbance. A disturbance unknown to MPC that affects the plant output. MPC provides feedback compensation for such disturbances. r Setpoint (or reference signal). The target value for the controlled variable. u Actuator (or manipulated variable). The signal MPC adjusts in order to achieve its control objectives. v Measured disturbance (optional). MPC feedforward compensation adjusts for such disturbances as they occur to minimize their impact on the controlled variable. Controlled variable (or plant output). The signal to be held at the setpoint. This is the “true” value, uncorrupted by measurement noise. y Measurement of the controlled variable. Used to estimate the true value, . z Measurement noise. Represents electrical noise, sampling errors, drifting calibration, and other effects that reduce the accuracy of the measurement. y y y y vy→ uy→ y MPC Fundamentals 1-7 NOTE: One may also specify bounds on . These constraint-handling properties are a distinguishing feature of MPC and can be particularly valuable when one has multiple control objectives to be achieved via multiple adjustments. In reality, however, model imperfections, plant limitations, and unmeasured disturbances cause the measurement, y, to deviate from its expected value. Thus, MPC uses the output measurement and a “disturbance model” () to predict future changes in . It then uses its model to calculate appropriate adjustments (a form of “feedback” compensation). This calculation also considers the known constraints. Various types of “noise” can corrupt the measurement. The signal z in Figure 1-1 represents such effects. They may vary randomly about a zero mean or exhibit a non-zero, drifting bias. MPC uses a model in combination with its model to remove the estimated noise component of the measurement (“filtering”). The above feedforward/feedback actions are MPC’s “regulator” mode. MPC also has a “servo” mode, i.e., it adjusts u such that tracks a time-varying setpoint. The tracking accuracy depends on the plant characteristics (including constraints), the accuracy of the model, and whether or not future setpoint variations can be anticipated, i.e., known in advance. Details: SISO case Sampling period and sampling instants MPC operates at discrete intervals of Dt time units, the “sampling period.” Suppose that MPC starts at time t = 0. The “sampling instants” are the times at which MPC adjusts the manipulated variable, u. They are integer multiples of the sampling period: 0,Dt, 2Dt, 3Dt, ..., kDt, where the integer index, k, represents the current sampling instant. y dy→ y uy→ zy→ dy→ y uy→ 1 Introduction 1-8 Figure 1-2 The MPC problem at the the sampling instant t. An MPC sampling instant Figure 1-2 shows the state of a hypothetical SISO MPC system which has been operating for some time; k is the current sampling instant. The current measured output, yk, and previous measurements, yk-1, yk-2, ..., are known and are the filled circles in Figure 1-2 (a). The current measured disturbance, dk, and its past values are also known (not shown). (b) (a) input horizon output horizon t t+1 futurepast t t+1 t+m t+p Predicted outputs y(t+k|t) t+m t+p Manipulated (t+)uk Inputs Reference r(t+k) umin umax MPC Fundamentals 1-9 Figure 1-2 (b) shows MPC’s previous “moves,” uk-41, ..., uk-1, as filled circles. As is usually the case, a “zero-order hold” receives each move from MPC and sends that value to the plant continuously until the next sampling instant – note the resulting step-wise variations in Figure 1-2 (b). MPC must now calculate the current move, uk. It does so in two phases: 1 Estimation. In order to make an intelligent move, MPC needs to know the current state of the system. This includes the true value of the controlled variable, , and any internal plant variables that influence the future trend, , ..., . To accomplish this MPC uses all past and current measurements and the models , , , and . For details, see Chapter 4, “Disturbance Detection and Estimation”. 2 Optimization. Values of setpoints, measured disturbances, and constraints are specified over a finite “horizon” of future sampling instants, k+1, k+2, ..., k+P, where P (a finite integer ≥ 1) is the “prediction horizon” – see Figure 1-2 (a). MPC then computes the M moves uk, uk+1, ... uk+M-1, where M ( ≥ 1, ≤ P) is the “control horizon” – see Figure 1-2 (b). In our hypothetical example, P = 9 and M = 4. Suppose that the optimal sequence of moves is the series of four open circles in Figure 1-2 (b). MPC’s model predicts that this will result in the output sequence shown as the nine open circles in Figure 1-2 (a). These moves are optimal in the sense that no constraints are violated and the outputs track the setpoint “closely.” See section Optimization for a detailed definition of optimality. MPC sends only the move uk to the plant. The plant operates with this input until MPC’s next sampling instant, Dt time units later. MPC then obtains a new set of measurements and revises completely the plan it had formulated at the previous sampling instant, thus compensating for model error and unknown disturbances. Prediction and control horizons One might wonder why MPC bothers to optimize over P sampling periods into the future and calculate M moves when it discards all but the first move. Indeed, under certain conditions MPC gives identical results for P = M = 1 as it does for P = M = ∞. More often, however, the horizons have an important impact. Some examples are: yk yk 1+ ykP+ uy→ dy→ wy→ zy→ 1 Introduction 1-10 ? Constraints. MPC’s current move, uk, depends on future constraints, i.e., MPC plans ahead to avoid future constraint violations. ? Delays in the plant. Suppose that the plant includes a pure time delay; MPC’s current move has no effect until D sampling periods into the future. Then we need P ≥ D and M ≤ P ? D so that MPC can include the effect of each contemplated move in its optimization calculations. ? Other nonminimum phase plants. Consider a SISO plant with an inverse-response, e.g., its unit-step response is a decrease below the initial condition followed by an increase to a final value above the initial condition. Experience suggests the following rules of thumb: 1 Choose the MPC sampling period such that the plant’s open-loop settling time is approximately 20-30 sampling periods (i.e., the sampling period is approximately one fifth of the dominant time constant). 2 Choose P to be the number of sampling periods used in step 1. 3 Use a relatively small M, e.g., 3-5. The resulting MPC performance will usually be insensitive to minor adjustments in these M and P values, and further “tuning” should not be necessary. If performance is poor, it is probably due to other factors, such as modeling errors. Blocking In Figure 1-2 (b), M = 4 and P= 9, and MPC is optimizing the first M moves of the prediction horizon, after which the manipulated variable remains constant for P – M = 5 sampling instants. Figure 1-3 shows an alternative “blocked” strategy – again with M = 4 – in which the first planned move occurs at sampling instant k, the next at k+2, the next at k+4, and the final at k+6. A “block” is one or more successive sampling periods during which the manipulated variable is constant. The “block durations” are the number of sampling periods in each block. In Figure 1-3 the block durations are 2, 2, 2, and 3. (Their sum must equal P.) MPC Fundamentals 1-11 Figure 1-3 Example of blocking with M = [2, 2, 2, 3]. As for the default (unblocked) mode, only the current move, uk, actually goes to the plant, and MPC repeats the optimization of M moves at its next sampling instant. Thus, the signal going to the plant actually varies every Dt time units. So why use blocking? When P >> M (as is generally recommended), and all M moves are at the beginning of the horizon, the moves tend to be larger (because all but the final move last just one sampling period). Blocking often leads to smoother adjustments, all other things being equal. Optimization For a SISO plant, MPC determines its moves by solving the following optimization problem (formulated for the kth sampling instant): (1-1) such that (1-2) Sampling Instant umin umax Past Moves Planned Moves k +1 +2 +3 +4 +5 +6 +7 +8 +9-1-2-3-4 Min uk … ukP1–+,, w y r ki+ y?ki+–() 2 wu ?u ki1–+() 2+[] i 1= P ∑ y?ki+ fuk … uki1–+,,()= 1 Introduction 1-12 (1-3) (1-4) (1-5) where is the adjustment at sampling instant j, and and are non-negative “weights.” Equation (1-2) represents the model MPC uses to predict the controlled variable at instant k+i. In standard MPC it is a linear function of the adjustments . Note: The prediction is also a function of known and estimated disturbances, but these effects are constant, i.e., independent of the adjustments. They are implicit in Equation (1-2). Equations (1-3) – (1-5) are optional. They specify the bounds on the manipulated variable, the controlled variable, and the change in the manipulated variable. When P > M (the usual case), there are also P – M equality constraints of the form , i.e., to hold the manipulated variables constant for a specified sequence of sampling instants. This depends on the specified block durations. For example, the situation of Figure 1-3 would have , , etc. Choosing the optimization weights The main difficulty in an MPC design (other than modeling) is the selection of the weights and in Equation (1-1). Some trial-and-error tuning is often needed. The following guidelines may be helpful: ? A positive value of penalizes predicted deviations of the controlled variable from its setpoint. The larger this weight, the more closely MPC tries to track the setpoint. If , MPC ignores the setpoint. ? A positive value of penalizes planned changes in the manipulated variable. A moderate value makes the controller “cautious,” which usually increases its “robustness” (making it less sensitive to modeling errors – a positive effect). A large value makes it “sluggish,” however. In the extreme as goes to infinity, MPC ignores deviations of the controlled variable from its setpoint. umin uki+ umax for i = 0:P-1≤≤ ymin y?ki+ ymax for i = 1:P≤≤ ?uki+ ?umax for i = 0:P-1≤ ?uj uj uj 1––= wy wu uk … uki1–+,, y?ki+ ukj+ ukj1–+= uk 1+ uk= uk 3+ uk 2+= wy wu wy wy 0= wu wu MPC Fundamentals 1-13 ? In a SISO problem, only the ratio of to matters. A good strategy is to fix one of them at unity and vary the other, testing the performance in simulations to obtain the desired tradeoff between fast setpoint tracking (large ratio) and cautious adjustments of the manipulated variable (small ratio). (By default, the MPC GUI sets .) ? MPC does whatever it must in order to satisfy any “hard” inequality constraints, equations (1-3)–(1-5). The values of and become less important when constraints are encountered. wy wu wy 1 wu, 0.1== wy wu 1 Introduction 1-14 Model Predictive Control Simulink Library Alberto Bemporad, Manfred Morari, Larry Ricker September 2, 2000 This document describes a new Simulink library for the Model Predictive Control Toolbox. Please report any bug or comment to bemporad@aut.ee.ethz.ch. 1 Description The new MPC Simulink library contains four blocks, see Fig. 1. The MPC blocks are based on the model shown in Fig. 2. This consist of the LTI system 1 braceleftBigg x(k+1) = Ax(k)+B u u(k)+B v v(k)+B d d(k) y(k)=Cx(k)+D v v(k)+D d d(k) (1) where x(k) represents the state of the system, u(k) are manipulated variables (MV) or command inputs, v(k) is a vector of measured disturbances (MD), d(k) are unmeasured disturbances (UD), andy(k) is the output vector, which is composed of measured outputs (MY)y m (k) and unmeasured outputs (UY) y u (k). Note that the matrix D u is assumed to be zero, i.e. no direct feedthrough of MVs on the output vector. 2 The unmeasured disturbance d(k) is modeled as the output of the LTI system braceleftBigg x d (k+1) = Fx d (k)+Gn(k) d(k)=Hx d (k)+Kn(k) (2) System (2) is driven by the random Gaussian noisen(k), which has zero mean and unit covariance matrix. For instance, a step-like unmeasured disturbance is modeled as the output of an integrator. 1 Ifcontinuoustimemodelsaresupplied,internallythesearesampledwiththecontroller’ssamplingtime 2 Thisassumptionisjusti?edbythefollowingargument. Let T s bethesamplingtime. Thetime τ k requiredby the controller to compute u(k) after y m (k)=y m (kT s ) has been acquired, is always ?nite. Then, u(t) ≡ u(k)for kT s +τ k ≤ t<(k+1)T s +τ k+1 . By assuming that zero-order holders are present in the controlsystem, a direct feedthroughwouldgiveacontributionon y m (k)oftheformD u u(kT s )=D u u(k ?1). Therefore,themodelcanbe describedas braceleftbigg x(k+1) = Ax(k)+B u u(k)+B v v(k)+B d d(k) y(k)=Cx(k)+D u u(k ?1)+D v v(k)+D d d(k) Byextendingthestate x(k)toinclude u(k ?1),onereturnstotheform(1). 1 Figure 1: MPC Simulink Library. The MPC controller selects the input u(k) by solving the optimization problem ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? min u(k|k),...,?u(m?1+k|k) p?1 summationdisplay i=0 bardblw u i [u(k+i|k)?u target (k)]bardbl 2 +bardblw ?u i ?u(k+i|k)bardbl 2 + bardblw y i+1 [y(k+i+1|k)?r(k+i+1)]bardbl 2 +ρ epsilon1 epsilon1 2 subj. to ? ? ? ? ? ? ? ? ? ? ? ? ? u min i ≤ u(k+i|k) ≤u max i ?u min i ≤ ?u(k+i|k) ≤ ?u max i ,i=0,...,p?1 ?epsilon1+y min i ≤ y(k+i+1|k) ≤y max i +epsilon1 ?u(k+j|k)=0,j=m,...,p epsilon1≥ 0 (3) with respect to the sequence of input increments {?u(k|k),...,?u(m? 1+k|k)} and the slack variable epsilon1. In (3), ?(k+i|k) denotes the value predicted for time k+i based on the information available at time k; r(k) is the current sample of the output reference. While only the measured outputs are connected to the MPC Simulink block,r(k) is a reference for all the outputs (measured and unmeasured). When the referenceris not known in advance, the current referencer(k)isused over the whole prediction horizon, namely r(k+i+1)≡ r(k) in (3). The exploitation of future references in MPC is referred to as anticipative action, and is only possible when the reference is read from ?le during simulation. A similar anticipative action can be performed with respect to measured disturbances v(k), namely v(k+i)=v(k) if the measured disturbance is not known in advance (e.g. is coming from another Simulink block) or v(k+i) is obtained from a ?le. In the prediction, d(k+i) is instead obtained by setting n(k+i) ≡ 0 in (2). Input and input-variation constraints are treated as hard constraints, output constraints are considered as soft. This is done by introducing the slack variable epsilon1 ≥ 0, and prevents the MPC 2 Model Outputsy(k) (SS, TF, ZPK) Manipulated Variables u(k) Unmeasured Disturbances d(k) Measured Disturbances v(k) x(k) Disturbance Model (SS, TF, ZPK) n(k) x (k) d Figure 2: Model used by the MPC block for prediction/state estimation. controller to get stuck because of infeasibility of the optimization problem. Infeasibility might in fact occur when model and plant mismatch signi?cantly, or just because of numerical round o?. By default, ρ epsilon1 =10 5 max p?1 i=0 {w u i ,w ?u i ,w y i+1 }. u target (k) is a set-point for the input vector; 3 w y i ,w u i ,w ?u i are nonnegative weight coe?cients; and u min i , u max i ,?u min i ,?u max i , y min i , y max i are lower/upper bounds. Only ?u(k|k) is actually used to compute u(k). The remaining samples ?u(k+i|k)aredis- carded, and a new optimization problem based on y m (k+ 1) is solved at the next sampling step k+1. As the true statesx(k), x d (k) are not available to the controller, predictions are obtained from the state estimates ?x(k), ?x d (k). These are computed from the measured outputy m (k) by the linear state observer bracketleftBigg ?x(k|k) ?x d (k|k) bracketrightBigg = bracketleftBigg ?x(k|k?1) ?x d (k|k?1) bracketrightBigg +L[y m (k)?C m ?x(k|k?1)?D vm v(k)?D dm H?x d (k|k?1)] where ? m denotes the rows of C,D relative to measured outputs, and ?x(k+1|k)=A?x(k|k)+B u u(k)+B v v(k)+B d H?x d (k|k) ?x d (k+1|k)=Fx d (k|k) To prevent numerical di?culties in the absence of unmeasured disturbances, the gainLis designed as the Kalman gain for the extended model bracketleftbigg x(k+1) x d (k+1) bracketrightbigg = bracketleftbigg AB d H 0 F bracketrightbiggbracketleftbigg x(k) x d (k) bracketrightbigg + bracketleftbigg B u 0 bracketrightbigg u(k)+ bracketleftbigg B v 0 bracketrightbigg v(k)+ bracketleftbigg B d K G bracketrightbigg n(k)+Iξ(k) (4) y m (k)= bracketleftbig C m D dm H bracketrightbig bracketleftbigg x(k) x d (k) bracketrightbigg +D vm v(k)+D dm Kn(k)+ζ(k) whereξ,ζareadditionalunmeasureddisturbanceshavingcovariancematrixS 1 andS 2 respectively, the covariance matrix of n is assumed to be I,andξ, ζ, n are considered independent. Note that d(k) models both state disturbances and output disturbances (B d =0,D d negationslash=0). 3 One typically uses u target if the number of inputs is greater than the number of outputs, as a lower-priority set-point. 3 Table 1: MPC blocks. Block Description Inputs MPC Regulator Regulate the output to the origin y(t) MPC Controller Track an arbitrary output set-point r(t) y(t), r(t) MPC Controller+MV Track an arbitrary output set-point r(t) y(t), r(t), v(t) and reject a measured input disturbance v(t) MPC Controller+Anticipation Same as MPC Controller+MV, but r(t), v(t) y(t) are read from a .MAT ?le The algorithm implemented in the MPC block uses di?erent procedures depending on the presence of constraints. If all the bounds are in?nite, then the problem (3) is solved analytically, otherwise a Quadratic Programming (QP) solver is used. Since soft output constraints are used, the QP problem should never be infeasible. For robust- ness, in case this should happen for some reason, the previous optimal sequence is applied, i.e. u(k+1)=u(k)+?u(k+1|k), and a warning message is given. 2 Parameters and Dialog Box The parameters of the MPC block are entered through the mask reported in Fig. 3 by double- clicking the MPC block, and are described in Table 2. Table 2: MPC block parameters. Parameter Description Type Default model LTI model used for predictions LTI object mandatory p prediction horizon Real 10 moves (blocking) moves Real or Array p limits upper and lower bounds on u,?u, y Cell Array of Arrays Inf weights weights on u?u target ,?u, y?r Cell Array of Arrays {0,.1,1} filename MAT ?le containing r, v String " MPCadd additional MPC parameters Structure [] ? model is an LTI object which is used for predictions and Kalman ?ltering, where model.b collectsB u , B v ,B d (in general the columns can be scrambled, see mvindex, mdindex below), and similarly model.d collects D u =0,D v , D d . ? p≥ 1is the prediction horizon. The case p=Inf is not implemented yet. ? moves is either a number m ≤p of free moves (i.e. ?u(k+i|k)=0fori = m,...,p), or a vector of blocking moves (see e.g. CMPC.M in the MPC Toolbox Reference Guide). ? limits is a cell array of the form {Umin,Umax,DUmin,DUmax,Ymin,Ymax}, where Umin[i,:] =(u min i ) prime (and similarly for the other matrices, cfr. ulim in CMPC.M). 4 Figure 3: Mask of MPC block. 5 ? weights is a cell array of the form {Uweight,DUweight,Yweight}, where Uweight[i,:] =w u i (and similarly for the other matrices) 4 . ? filename is only present in the MPC block with anticipation, and speci?es the MAT ?le where the reference r(t) and the measured disturbance v(t) signals are stored. The format is the same as in the From File Simulink block. The ?rst variable saved in the MAT ?le is loaded from disk and is used to build the reference and measured disturbance signals. The variable is a matrix whose ?rst row is a vector of timest, the followingn y rows are the referencer(t), and the following n v rows are the measured disturbance v(t). Missing rows are treated as zeros. The signals are resampled with the MPC controller sampling time and stored in the MPC block memory. The ?rst (last) sample is used for simulation time before (after) the speci?ed range of t. Note that the same MAT ?le should be used for generating the actual signals r, v in the simulation, so that the MPC block has a consistent information (see Figure 6). ? MPCadd is a structure containing the ?elds indicated in Table 3. Table 3: Structure of additional MPC parameters. Field Description Type Default utarget inputset-point Array 0 mvindex indicesofmanip. vars(withininputvector) Array mandatory mdindex indicesofmeas. disturb. (withininputvector) Array [] myindex indicesofmeas. outputs(withinoutputvector) Array mandatory noisemodel LTImodelforunmeasureddisturbances Array [] covmat estimatorcovariancematrices Cell Array of Arrays {} rhoeps weightonslackvariable ρ epsilon1 Real 10 5 max p?1 i=0 {w u i ,w ?u i ,w y i+1 } x0 initialstateestimate ?x(0) Array zeros u1 initialpastinput u(?1) Array zeros ts controllersamplingtime Real model.Ts delay simulatedcomputationaldelay Real 0 – MPCadd.utarget is a row vector of set-points for the input vector, as described above. – MPCadd.mvindexandMPCadd.mdindexarerowvectorsindicatingthepositionofMVand MDvariables(u,v)withintheoverallinputvector,forinstanceB u =model.b(:,MPCadd. mvindex). – MPCadd.myindex is a row vector indicating the positions of MY variablesy m within the output vector y. – MPCadd.noisemodel is an LTI object which is used for estimating the unmeasured dis- turbance d(k) (see (2) and Fig. 2). If empty, step-like unmeasured disturbances are assumed by default, i.e. one integrator for each unmeasured disturbance. – MPCadd.covmat is a cell array de?ning the covariance matricesS 1 ,S 2 for the additional disturbancesξ(t),ζ(t) de?ned in (4). By default, if MPCadd.covmat={} then step distur- 4 If Umin, Uweight haveanumberofrowssmallerthan p,thelastrowisautomaticallyrepeated(thesameforthe otherlimitandweightmatrices) 6 bances are added on each output (i.e. integrators driven by white noise of unit covari- ance matrix), S 1 =I, S 2 =I (if there are no unmeasured disturbances), or S 1 =.001I, S 2 = .1I (if there are unmeasured disturbances). If MPCadd.covmat={S1,S2}, then S1 (S2) can be either [] (default value), a scalar (S 1 =S1I), or a matrix. – MPCadd.rhoeps is the weight ρ epsilon1 on the slack variable epsilon1 representing the violation of the constraints. The larger ρ epsilon1 with respect to input and output weights, the more the constraint violation is penalized. Although the same variable epsilon1 is used for softening all output constraints, di?erent penalties can be induced by scaling the output components di?erently. For instance, if the constraint y 1 > 0 is much more important than the constraint ?1 ≤ y 2 ≤ 1, the ?rst row of the C and D matrices of the model can be multiplied by a large scalar, the constraints y min 1 , y max 1 by the same scalar, and the output measurement and reference signal ampli?ed by the same amount before being connected to the MPC controller, for consistency. Note that the output weights must be rescaled as well, accordingly. – MPCadd.x0 is the initial value for the state estimate ?x(0), and u1 the initial value of the input u at time ?1(as u(0) = ?u(0)+u(?1), nonzero values for u1 might be needed if there are tight limits on ?u). The default is ?x(0) = 0. – MPCadd.ts is the sampling time of the controller: the MPC algorithm is executed every MPCadd.ts units of time. If MPCadd.ts is not given, then MPCadd.ts is automatically inherited from the sampling time model.Ts.IfMPCadd.tsnegationslash=model.Ts, then the model is resampled with sampling time MPCadd.ts. – MPCadd.delay,0≤MPCadd.delay≤MPCadd.ts is a time o?set for the execution of the MPCcontroller, i.e.u(t)=u(k)fort∈ [kMPCadd.ts+MPCadd.delay,(k+1)MPCadd.ts+ MPCadd.delay). This allows to take into account in simulation computational delays occurring in real operation. In order to use nonzero MPCadd.delay,aZero-Order Hold Simulink block with sampling time=MPCadd.delay must be inserted between the y, r, v signals and the MPC block (see EX7DELAY.M). 3 Example 1 The script EX1.M prepares the parameters to simulate the Simulink model MPC EX1.MDL,depicted in Fig. 4. % 3x1 system % % MV=[1], MD=[2], UD=[3] % Y=[1] % Prepare parameters for mpc_ex1 % True plant & true initial state sys=ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]})); A=sys.a; 7 Figure 4: Simulink model MPC EX1.MDL. B=sys.b; C=sys.c; D=sys.d; x0=[0 0 0 0 0]’; Ts=.2; %Sample time Tstop=30; %Simulation time % prediction model model=c2d(sys,Ts); p=10; moves=3; 8 umin=0; umax=1; dumin=-Inf; dumax=Inf; ymin=-Inf; ymax=Inf; limits={umin,umax,dumin,dumax,ymin,ymax}; weights={[0],[.1],[1]}; clear MPCadd MPCadd.mvindex=[1]; MPCadd.mdindex=[2]; open_system(’mpc_ex1.mdl’) sim(’mpc_ex1’,Tstop) The LTI system model used for the MPC block is obtained by sampling the LTI continuous time system sys. This is also used as the plant to be controlled in Fig. 4. The reference trajectory r(t) is a unit step. A step measured disturbancev(t) of intensity 1appears at time t= 10, while a step unmeasured disturbanced(t)ofintensity?0.5 appears at timet= 20. The input variableu(t) is constrained within the range [0,1]. The input and output trajectories are reported in Fig. 5. 3.1 Anticipative Action The script EX1A.M simulates the Simulink model MPC EX1A.MDL, depicted in Fig. 6, which is a modi?edversionofMPC EX1.MDL where r and v signals are read from the MAT ?le RV.MAT,and MPC uses an anticipative action. Note that the same MAT ?le is used for generating the actual signals r, v in the simulation. % Reference and measured disturbance signals time=0:Ts:Tstop; trv=[time; zeros(size(0:Ts:5)),ones(size(5+Ts:Ts:Tstop)) zeros(size(0:Ts:10)),ones(size(10+Ts:Ts:Tstop))]; save rv trv The input and output trajectories are reported in Fig. 7. 4 Example 2 The script EX2NL.M prepares the parameters to simulate the Simulink model MPC EX2NL.MDL,de- picted in Fig. 8. 9 % MIMO (3x2) SYSTEM: MV=[1 2 3], MY=[1 2] % UD=[5 6] % Prepare parameters for mpc_ex2nl % Nonlinear model sys=’nl3x2’; % Linearize the model at (0,0) [A,B,C,D]=linmod2(sys); Ts=.2; %Sample time model=c2d(ss(A,B,C,D),Ts); %Convert to discrete time clear A B C D p=5; %Prediction horizon moves=2; %Free control moves % Define lower and upper bounds for u, delta u, y umin=[-2,-1,-1]; umax=[2,1,1]; dumin=[-Inf,-Inf,-Inf]; dumax=[Inf,Inf,Inf]; ymin=[-Inf,-Inf]; ymax=[Inf,Inf]; limits={umin,umax,dumin,dumax,ymin,ymax}; % Define weights for u, delta u, y weights={[0 0 0],[.1 .1 .1],[1 1]}; % No noisemodel, no covmat supplied implies that % output step disturbances will be added by default open(’mpc_ex2nl.mdl’) % Display the Simulink model % Run simulation Tfinal=8; sim(’mpc_ex2nl’,Tfinal) 10 The LTI system model used for the MPC block is obtained by converting to discrete time the Simulink nonlinear model NL3X2.MDL depicted in Fig. 9. A copy of NL3X2.MDL is also used as the plant to be controlled in Fig. 8. In order to take into account steady-state o?sets due to the nonlinearity, a step-like unmeasured disturbances is added on each output. The reference trajectoriesr 1 (t)andr 2 (t) are unit step, andr 1 (t) switches back to 0 at timet= 4. Constraints are present on the input variables u 1 (t), u 2 (t), u 3 (t). The input and output trajectories are reported in Fig. 10. 5 Upgrading from MPC Toolbox v1.0 ? In the MPC Toolbox v1.0, the variable uwt refers to weights on increments of manipulated variables. Here, the corresponding variable has been called DUWeight. 6 Installation The provided ?le archive must be extracted into a directory (e.g. MPCBLOCK/). This directory must be added to Matlab’s path. The archive contains the following ?les: mpc1.dll mpc2.dll mpcbl1.tif mpcbl2.tif mpcbl3.tif mpcbl4.tif mpclib.mdl mpcsfun.m Directory of \demos ex1.m ex10.m ex11.m ex12.m ex1a.m ex2nl.m ex3.m ex6.m ex7.m ex7delay.m ex7u.m ex9.m mpc_ex1.mdl mpc_ex10.mdl mpc_ex11.mdl 11 mpc_ex12.mdl mpc_ex1a.mdl mpc_ex2nl.mdl mpc_ex3.mdl mpc_ex6.mdl mpc_ex7.mdl mpc_ex7delay.mdl mpc_ex7u.mdl mpc_ex9.mdl nl3x2.mdl Directory of \manual mpcblock.pdf mpcblock.ps 12 (a)Output y(t)andreference r(t). (b)Manipulatedvariable u(t). (c)Measureddisturbance v(t). (d)Unmeasureddisturbance d(t). Figure 5: Trajectories obtained running EX1.M. 13 Figure 6: Simulink model MPC EX1A.MDL. 14 (a)Output y(t)andreference r(t). (b)Manipulatedvariable u(t). (c)Measureddisturbance v(t). (d)Unmeasureddisturbance d(t). Figure 7: Trajectories obtained running EX1A.M. 15 Figure 8: Simulink model MPC EX2NL.MDL. Figure 9: Simulink model NL3X2.MDL. 16 (a)Outputs. (b)Inputs. Figure 10: Trajectories obtained running EX2NL.M. 17