Patent Publication Number: US-7587253-B2

Title: Partial enumeration model predictive controller

Description:
STATEMENT OF GOVERNMENT RIGHTS 
   This invention was made with United States government support awarded by the following agency: NSF #CTS-0105360. The United States government has certain rights in this invention. 

   TECHNICAL FIELD 
   This invention relates to system controllers, and more particularly, to a partial enumeration model predictive controller. 
   BACKGROUND OF THE INVENTION 
   Systems can be defined generally as an object in which signals interact to produce an observable output signal. The system can be a physical entity, such as a chemical process, an electrical circuit, or an engine. The system can also be an abstract entity, such as the stock market or a financial system. One important type of system is the multiple input, multiple output (MIMO) system. As the name suggests, a MIMO system has multiple inputs and multiple outputs, in contrast to single input, single output (SISO) systems. MIMO systems are commonly controlled using model-based control methods, such as linear quadratic (LQ) control or model predictive control (MPC). 
   Model predictive control is the dominant advanced control method in process industries for control of large-scale MIMO processes. MPC systems use a model of the MIMO system and make control decisions based on the model. MPC systems are popular for complex systems because they can explicitly optimize the process, handle complex multivariate processes, and account for system constraints, such as valve motion limits. In spite of these benefits, existing MPC systems present problems that prevent widespread use. MPC systems run too slowly for large-scale processes, processes with long time delays, and processes with fast sampling times. MPC systems also run too slowly when the problem being solved has long control horizons. Specific problems that cause the MPC systems to run slowly include solution tables becoming so large as to prevent searching in a reasonable time and solution tables lacking the optimal solution. Performance of the MPC system often suffers when a suboptimal solution is used. 
   Previous attempts to use MPC systems for complex MIMO processes have used complete enumeration strategies, i.e., strategies which solve the MIMO system over the whole parameter region for all possible solutions. One example of a complete enumeration strategy for solving the MPC systems formulated as quadratic programs is a multi-parametric quadratic programming solver technique, which can use active set techniques or interior point quadratic programming solver techniques to compute the optimal solution for each current parameter. The number of computations involved in these complete enumeration strategies grows rapidly with the number of dimensions and length of the horizon, making the strategies slow to run and unsuitable for real time control of complex MIMO processes. 
   It would be desirable to have a partial enumeration model predictive controller that would overcome the above disadvantages. 
   SUMMARY OF THE INVENTION 
   One aspect of the present invention provides a method of predictive control for a multiple input, multiple output (MIMO) system, including providing a solution table with problem solutions to a model predictive control problem for the MIMO system over a partial parameter region, scanning the solution table for an optimal solution for current parameters, using the optimal solution to control the MIMO system when the optimal solution is in the solution table, and using an alternative solution to control the MIMO system when the optimal solution is not in the solution table. 
   Another aspect of the present invention provides a system of predictive control for a multiple input, multiple output (MIMO) system, including means for providing a solution table with problem solutions to a model predictive control problem for the MIMO system over a partial parameter region, means for scanning the solution table for an optimal solution for current parameters, means for using the optimal solution to control the MIMO system when the optimal solution is in the solution table, and means for using an alternative solution to control the MIMO system when the optimal solution is not in the solution table. 
   Another aspect of the present invention provides a computer readable medium storing a computer program for predictive control for a multiple input, multiple output (MIMO) system, including computer readable code for providing a solution table with problem solutions to a model predictive control problem for the MIMO system over a partial parameter region, computer readable code for scanning the solution table for an optimal solution for current parameters, computer readable code for using the optimal solution to control the MIMO system when the optimal solution is in the solution table, and computer readable code for using an alternative solution to control the MIMO system when the optimal solution is not in the solution table. 
   The foregoing and other features and advantages of the invention will become further apparent from the following detailed description of the presently preferred embodiments, read in conjunction with the accompanying drawings. The detailed description and drawings are merely illustrative of the invention, rather than limiting the scope of the invention being defined by the appended claims and equivalents thereof. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  is a block diagram of a partial enumeration model predictive controller made in accordance with the present invention; 
       FIG. 2  is a flow chart of a method of model predictive control for a MIMO system made in accordance with the present invention; 
       FIG. 3  is a hardware diagram of a MIMO model predictive controller made in accordance with the present invention; and 
       FIGS. 4 &amp; 5  are simulation results showing the performance of a MIMO model predictive controller made in accordance with the present invention. 
   

   DETAILED DESCRIPTION OF THE INVENTION 
     FIG. 1  is a block diagram of a partial enumeration model predictive controller made in accordance with the present invention. The partial enumeration model predictive controller  100  controls a system  102  having an input u k  and an measured output y k . The partial enumeration model predictive controller  100  includes a state and disturbance estimator  104  receiving the input u k  and the output y k , and generating an estimated state x*hat*k and a estimated disturbance d*hat*k; a constrained target calculator  106  receiving the estimated disturbance d*hat*k and a steady state target output z*bar, and generating a steady state target state x*bar*k and a steady state target input u*bar*k; and a constrained dynamic optimizer  108  receiving the estimated state x*hat*k, the estimated disturbance d*hat*k, the steady state target state x*bar*k, and the steady state target input u*bar*k, and generating an optimized input u k . Those skilled in the art will appreciate that there is an implicit time shift in the variables of  FIG. 1 . The u k  received at the state and disturbance estimator  104  is the previous input, i.e., u k−1 . The output y k , estimated state x*hat*k, and estimated disturbance d*hat*k depend on u k−1 . The steady state target state x*bar*k and steady state target input u*bar*k are computed and a present u k  computed, which will affect y k+1 . 
   The system  102  is any multiple input, multiple output (MIMO) system, having an input u k  and output y k . The system  102  can be linear or non-linear. The state of the system  102  is described with a state vector x k . The system can be a physical entity, such as a chemical process, an electrical circuit, or an engine, or an abstract entity, such as a financial system. A linear time-invariant model augmented with integrating states to remove offset for the system  102  can be written as:
 
 x   k+1   =Ax   k   +Bu   k   +B   d   d   k  
 
 d   k+1   =d   k  
 
 y   k   =Cx   k   +C   d   d   k ,  (1)
 
   where the state x is a real state vector, the input u is a real input vector, the output y is a real output vector, and the disturbance d is a real disturbance vector. A, B, C, B d , and C d  are fixed matrix model factors representing the system  102 . The sampling time k is a non-negative integer. The state x and input u is constrained without loss of generality by the following hard constraint:
 
 Du   k   −Gx   k   ≦h   (2)
 
   where D and G are matrices of appropriate dimensions. The hard constraint accounts for limitations on the input to the system  102 , such as input bound and input rate-of-change constraints. Examples of such constraints are valve travel limits in a chemical process or current limits in an electrical circuit. Those skilled in the art will appreciate that the constraints can be hard constraints, such as the hard constraint of this example, and/or soft constraints, such as output constraints. Those skilled in the art will also appreciate that the system  102  can be any number of MIMO systems. 
     FIG. 2  is a flow chart of a method of model predictive control for a MIMO system made in accordance with the present invention. The method is a partial enumeration method, i.e., problem solutions are provided over a partial parameter region rather than the including all problem solutions. 
   The method includes providing a solution table with problem solutions to a model predictive control problem for the MIMO system over a partial parameter region  210 , scanning the solution table for an optimal solution for current parameters  220 , using the optimal solution to control the MIMO system when the optimal solution is in the solution table  230 , and using an alternative solution to control the MIMO system when the optimal solution is not in the solution table  240 . The method can then continue scanning the solution table for the next optimal solution for the next sampling time  250 . The alternative solution can be a best suboptimal solution, a predetermined time current solution, or a simplified model current solution. The solution table can be provided by modeling the MIMO system with model factors, generating the model predictive control problem from the model factors, solving the model predictive control problem for the problem solutions offline, and populating the solution table with the problem solutions. 
   The method can further include calculating a current solution to the model predictive control problem for the MIMO system for the current parameters when the optimal solution is not in the solution table and adding the current solution to the solution table. An old problem solution, such as a least used solution or a least recently used solution, can be deleted when the current solution is added to maintain the solution table at a predetermined length. The addition to the solution table can also be a solution family related to the current solution. The method can further include determining a use score for each of the problem solutions over a predetermined number of sampling times, so that the scanning the solution table for an optimal solution includes scanning the problem solutions in order of the use score from most used to least used. The method can further include determining an age for each of the problem solutions, so that the scanning the solution table for an optimal solution includes scanning the problem solutions in order of the age from newest to oldest. 
   Providing a solution table with problem solutions to a model predictive control problem for the MIMO system over a partial parameter region  210  typically includes modeling the MIMO system with model factors, generating the model predictive control problem from the model factors, determining the partial parameter region for an expected operating state, solving the model predictive control problem offline for the problem solutions in the partial parameter region, and populating the solution table with the problem solutions. Modeling the MIMO system with model factors includes modeling the MIMO system as a linear time-invariant model with model factors as discussed above for  FIG. 1 . 
   Generating the model predictive control problem from the model factors includes formulating a target calculation that computes the steady state target state and the steady state target input. The estimator calculates estimated state x*hat*k|k and estimated disturbance d*hat*k|k, based on the current measured output y k . The model assumes that the estimated state x*hat*k|k and estimated disturbance d*hat*k|k are available at every sampling time. The model also assumes that it is possible to make a linear transformation z=C z y of the output y with known setpoints z*bar. The z is typically a subvector of the output y. 
   The target calculator computes the steady state target state and the steady state target input at each sampling time from the estimated state x*hat*k|k, estimated disturbance d*hat*k|k, and desired output z*bar. The target calculation computes the targets that may drive the controlled variables to their setpoints within the constraints, such as the hard constraints of Equation (2). The target calculation formulation can be selected as desired from a number of possible formulations, such as linear objective functions, quadratic objective functions, hard constraint functions with precedence rank ordering, soft constraint functions, multi-stage quadratic objective functions, or the like. 
   In this example, the target calculation formulation is a two-stage quadratic objective functions approach. A first quadratic program is written as: 
   
     
       
         
           
             
               
                 
                   min 
                   
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           u 
                           _ 
                         
                         k 
                       
                       - 
                       
                         
                           x 
                           _ 
                         
                         k 
                       
                     
                   
                 
                 ⁢ 
                 
                   
                     1 
                     2 
                   
                   ⁢ 
                   
                     
                       ( 
                       
                         
                           
                             u 
                             _ 
                           
                           k 
                         
                         - 
                         
                           u 
                           _ 
                         
                       
                       ) 
                     
                     ′ 
                   
                   ⁢ 
                   
                     
                       R 
                       s 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             u 
                             _ 
                           
                           k 
                         
                         - 
                         
                           u 
                           _ 
                         
                       
                       ) 
                     
                   
                 
               
             
             
               
                 ( 
                 
                   3 
                   ⁢ 
                   a 
                 
                 ) 
               
             
           
         
       
     
   
   subject to
 
   x     k   =A  x     k   +Bū   k   +B   d   {circumflex over (d)}   k|k   (3b)
 
   z =C   z ( C  x     k   +C   d   {circumflex over (d)}   k|k )  (3c)
 
 Dū   k   −G  x     k   ≦h   (3d)
 
   where R s  is a positive definite matrix, u*bar is the steady state desired input, and z*bar is the steady state desired output setpoint. The steady state desired input u*bar and the steady state desired output setpoint z*bar are desired values, whereas the results of the target calculation, steady state target input u*bar*k and steady state target state x*bar*k, are the targets. 
   In some cases, the first quadratic program of Equation (3) is infeasible, i.e., the estimated state, estimated disturbance, and steady state target output may not be able to provide a satisfactory steady state target state and the steady state target input satisfying the first quadratic program due to conflicting constraints and/or requirements. When the first quadratic program of Equation (3) is infeasible, a second quadratic program is provided to supply an offset, which allows the controlled variables z to track to the setpoint vector z*bar. 
   The second quadratic program aimed at minimizing the offset is written as: 
   
     
       
         
           
             
               
                 
                   min 
                   
                     
                       
                         
                           
                             u 
                             _ 
                           
                           k 
                         
                         - 
                         
                           
                             x 
                             _ 
                           
                           k 
                         
                       
                       , 
                       
                         
                           z 
                           _ 
                         
                         k 
                       
                     
                     ⁢ 
                     
                         
                     
                   
                 
                 ⁢ 
                 
                   
                     1 
                     2 
                   
                   ⁢ 
                   
                     { 
                     
                       
                         
                           
                             ( 
                             
                               
                                 
                                   z 
                                   _ 
                                 
                                 k 
                               
                               - 
                               
                                 z 
                                 _ 
                               
                             
                             ) 
                           
                           ′ 
                         
                         ⁢ 
                         
                           
                             Q 
                             s 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   z 
                                   _ 
                                 
                                 k 
                               
                               - 
                               
                                 z 
                                 _ 
                               
                             
                             ) 
                           
                         
                       
                       + 
                       
                         
                           
                             ( 
                             
                               
                                 
                                   u 
                                   _ 
                                 
                                 k 
                               
                               - 
                               
                                 u 
                                 _ 
                               
                             
                             ) 
                           
                           ′ 
                         
                         ⁢ 
                         
                           
                             R 
                             s 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   u 
                                   _ 
                                 
                                 k 
                               
                               - 
                               
                                 u 
                                 _ 
                               
                             
                             ) 
                           
                         
                       
                     
                     } 
                   
                 
               
             
             
               
                 ( 
                 
                   4 
                   ⁢ 
                   a 
                 
                 ) 
               
             
           
         
       
     
   
   subject to:
 
   x     k   =A  x     k   +Bū   k   +B   d   {circumflex over (d)}   k|k   (3b)
 
 Dū   k   −G  x     k   ≦h   (3d)
 
   z     k   =C   z ( C  x     k   +C   d   {circumflex over (d)}   k|k )  (4b)
 
   where Q is a positive definite matrix. 
   The steady state target input u*bar*k and the steady state target state x*bar*k can be calculated at each sampling time from the first or second quadratic programs above. For convenience of notation, the system variables are reparameterized with a stabilizing gain matrix K and deviation variables are defined as follows:
 
 w   j   ={circumflex over (x)}   k+j|k   −  x     k   , v   j =( u   k+j    −ū   k )− Kw   j .  (5)
 
   Using the deviation variables, an MPC optimization problem can be defined as: 
   
     
       
         
           
             
               
                 
                   
                     min 
                     
                       w 
                       , 
                       v 
                     
                   
                   ⁢ 
                   
                     
                       ∑ 
                       
                         j 
                         = 
                         0 
                       
                       
                         N 
                         - 
                         1 
                       
                     
                     ⁢ 
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         { 
                         
                           
                             
                               w 
                               j 
                               ′ 
                             
                             ⁢ 
                             
                               Qw 
                               j 
                             
                           
                           + 
                           
                             
                               v 
                               j 
                               ′ 
                             
                             ⁢ 
                             
                               Rv 
                               j 
                             
                           
                           + 
                           
                             2 
                             ⁢ 
                             
                               w 
                               j 
                               ′ 
                             
                             ⁢ 
                             
                               Mv 
                               j 
                             
                           
                         
                         } 
                       
                     
                   
                 
                 + 
                 
                   
                     1 
                     2 
                   
                   ⁢ 
                   
                     w 
                     N 
                     ′ 
                   
                   ⁢ 
                   
                     Pw 
                     N 
                   
                 
               
             
             
               
                 ( 
                 
                   6 
                   ⁢ 
                   a 
                 
                 ) 
               
             
           
         
       
     
   
   subject to:
 
 w   j+1   =A   K   w   j   +Bv   j   , j= 0, . . . ,  N− 1,  (6b)
 
 Dv   j   −G   K   w   j   ≦h −( Dū   k   −G  x     k ),  j= 0, . . . ,  N− 1,  (6c)
 
   where A K  is a strictly stable matrix defined as A K =A+BK; G K  is defined as G K =G−DK; Q and P are positive semidefinite matrices; and R is a positive definite matrix. 
   In this formulation of the MPC optimization problem, the parameters which may change with each solution at each sampling time are b s  and w 0 , where:
 
 b   s   =[ū′   k     x ′   k ]′.
 
   The other quantities Q, R, M, P, A K , B, D, G, and G K  can be considered constant, i.e., independent of time. The state variables
 
 w=[w′   1    . . . w′   N ]′
 
   can be eliminated by applying Equation (6b) to the MPC optimization problem, arriving at the following strictly convex quadratic program for the MPC optimization problem: 
   
     
       
         
           
             
               
                 
                   
                     min 
                     v 
                   
                   ⁢ 
                   
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       v 
                       ′ 
                     
                     ⁢ 
                     Hv 
                   
                 
                 + 
                 
                   
                     g 
                     ′ 
                   
                   ⁢ 
                   v 
                 
               
             
             
               
                 ( 
                 
                   7 
                   ⁢ 
                   a 
                 
                 ) 
               
             
           
         
       
     
   
   subject to:
 
 Av≧b   (7b)
 
   where H is a positive definite constant matrix; A is a constant matrix; and
 
 v=[v′   0    . . . v′   N−1 ]′
 
   The vectors g and b depend on the parameters b s  and w 0  as follows:
 
 g=G   w   w   0   (8a)
 
 b=  b +B   s   b   x   +B   w   w   0   (8b)
 
   where b*bar is a constant vector, and G w , B s , and B w  are constant matrices. 
   Assuming the optimal solution v* obtained from the formulation of the MPC optimization problem of Equation (6) or (7), the first input u k  is obtained from Equation (5) for injection into the system as:
 
 u   k   =ū   k   +Kw   0   +v*   0 ,  (9)
 
   Solving the model predictive control problem offline for the problem solutions in the partial parameter region includes solving the formulation of the MPC optimization problem of Equation (6) or (7) once at each sampling time. MPC optimization problems for consecutive sampling times differ only in the parameters b s  and w 0 . The unique solution to the strictly convex formulation of the MPC optimization problem of Equation (7) satisfies the following Karush-Kuhn-Tucker (KKT) optimality conditions for some choice of active constraint indices a , whose complement I denotes the inactive constraint indices:
 
 Hv*+g−A′   a λ* a =0  (10a)
 
 A   a   v*=b   a   (10b)
 
 A   i   v*≧b   i   (10c)
 
λ* a ≧0,  (10d).
 
   Assuming that the matrix A a  has full row rank, i.e., that the active constraints are linearly independent, Equations (10a) and (10b) provide a linear system with a square nonsingular coefficient matrix as follows: 
   
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           H 
                         
                         
                           
                             - 
                             
                               A 
                               a 
                               ′ 
                             
                           
                         
                       
                       
                         
                           
                             A 
                             a 
                           
                         
                         
                           0 
                         
                       
                     
                     ] 
                   
                   ⁡ 
                   
                     [ 
                     
                       
                         
                           
                             v 
                             * 
                           
                         
                       
                       
                         
                           
                             λ 
                             a 
                             * 
                           
                         
                       
                     
                     ] 
                   
                 
                 = 
                 
                   
                     [ 
                     
                       
                         
                           
                             - 
                             g 
                           
                         
                       
                       
                         
                           
                             b 
                             a 
                           
                         
                       
                     
                     ] 
                   
                   . 
                 
               
             
             
               
                 ( 
                 11 
                 ) 
               
             
           
         
       
     
   
   Equation (11) uniquely determines variables v* and λ a  *. 
   The dependence of the formulation of the MPC optimization problem solely on the parameters b s  and w 0  can be used to express the solution explicitly in terms of the parameters b s  and w 0 . The parameter region of (b s , w 0 ) can be partitioned into polyhedral regions within which the parameters are valid. Substituting Equation (8) into the linear system of Equation (11) provides: 
   
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         
                           
                             H 
                           
                           
                             
                               - 
                               
                                 A 
                                 a 
                                 ′ 
                               
                             
                           
                         
                         
                           
                             
                               A 
                               a 
                             
                           
                           
                             0 
                           
                         
                       
                       ] 
                     
                     ⁡ 
                     
                       [ 
                       
                         
                           
                             
                               v 
                               * 
                             
                           
                         
                         
                           
                             
                               λ 
                               a 
                               * 
                             
                           
                         
                       
                       ] 
                     
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             0 
                           
                         
                         
                           
                             
                               
                                   
                               
                               ⁢ 
                               
                                 
                                   b 
                                   _ 
                                 
                                 a 
                               
                             
                           
                         
                       
                       ] 
                     
                     + 
                     
                       
                         [ 
                         
                           
                             
                               0 
                             
                           
                           
                             
                               
                                 B 
                                 as 
                               
                             
                           
                         
                         ] 
                       
                       ⁢ 
                       
                         b 
                         s 
                       
                     
                     + 
                     
                       
                         [ 
                         
                           
                             
                               
                                 - 
                                 
                                   G 
                                   w 
                                 
                               
                             
                           
                           
                             
                               
                                 B 
                                 aw 
                               
                             
                           
                         
                         ] 
                       
                       ⁢ 
                       
                         w 
                         0 
                       
                     
                   
                 
                 , 
               
             
             
               
                 ( 
                 12 
                 ) 
               
             
           
         
       
     
   
   where b a  *bar, B as , and B aw  are row subvectors and submatrices corresponding to the active set a of problem solutions. Continuing with the assumption that the matrix A a  has full row rank, the solution of Equation (12) depends linearly on the parameters b s  and w 0 , so that:
 
 v*=K   as   b   s   +K   aw   w   0   +c   v ,  (13a)
 
λ* a    =L   as   b   s   +L   aw   w   0   +c   λ ,  (13b)
 
   where K as , K aw , L as , L aw , c v  and c λ  depend on the active set a of problem solutions and problem data, such as the system matrices A, B, C, the weight matrices Q, R, M, P, and so forth, but not on the parameters b s  and w 0 . The region of validity of the solutions to Equation (13) can be determined from the unused optimality conditions of Equations (10c) and (10d). Substituting Equation (13) into the optimality conditions of Equations (10c) and (10d) provides explicit tests for the parameters as follows:
 
 L   as   b   s   +L   aw   w   0   +c   λ ≧0,  (14a)
 
 A   l   K   as   b   s   +A   l   K   aw   w   0 +( A   l   c   v   −b   l )≧0,  (14b).
 
   Those skilled in the art will appreciate that any technique with redundant constraint removal can be used to eliminate redundant inequalities in Equation (14). One example of such an algorithm is described in Tondel, Johansen, and Bemporad,  Automatica , vol. 39, p. 489-497, 2003. 
   Solving the model predictive control problem for the problem solutions offline further includes calculating and storing the coefficient matrices and vectors of Equations (13) and (14) for a selected subset of possible active sets a of problem solutions that are valid active sets for some choice of the parameters b s  and w 0 . In one embodiment, it is not necessary to compute (13b). The parameters b s  and w 0  can be selected for the initial anticipated operating state for the system, although the selection is not critical as new solutions are calculated and stored during regular operation when the system moves to a different operating state. The selected subset of possible active sets a of problem solutions is determined by solving the MPC problem for a training period over a desired number of sampling times to simulate MPC operation, while adding random disturbances and setpoint changes to force the MPC problem into the region of the parameters b s  and w 0  for the initial anticipated operating state. Populating the solution table with the problem solutions includes storing the coefficient matrices and vectors of Equations (13) and (14) for a selected subset of possible active sets a of problem solutions in the solution table. Typically, the solution table is small and of fixed length, so the selected subset is small. The size of the solution table can be selected to be large enough to fit the range of optimal active sets a of problem solutions likely to be encountered for the current operating state. This maintains the number of times that the optimal solution is not found in the solution table to a reasonable level. Exemplary solution tables accommodate a selected subset for 50 to 200 possible active sets a of problem solutions. Those skilled in the art will appreciate that the number of possible active sets depends on the control horizon N and on the system dimensions n,m. These parameters can be used to compute an upper bound to the number of possible active sets, although this bound is often very conservative. In another embodiment, the solution table is allowed to grow with addition of new optimal active sets of problem solutions for solutions not initially stored in the solution table. 
   Scanning the solution table for an optimal solution for current parameters  220  includes determining the parameters b s  and w 0  for the current operating state and checking the solution table for an active set a of problem solutions corresponding to the current parameters. The checking can be performed by using matrix-vector products to see whether the inequalities of Equation (14) hold for the active set a of problem solutions. The optimal v* can be calculated from Equation (13) for the identified active set a problem solution. When the optimal solution is in the solution table, the identified optimal solution is used to control the MIMO system  230 . 
   When the optimal solution is not in the solution table, an alternative solution is used to control the MIMO system  240 . Using an alternative solution when the optimal solution is not in the solution table to control the MIMO system  240  includes identifying an alternative solution such as a best suboptimal solution, a predetermined time current solution, a simplified model current solution, or the like. The alternative v* can be calculated from Equation (13) for the identified active set a problem solution. 
   The best suboptimal solution can be determined by scanning the solution table for an active set a problem solution that ensures the solution of the quadratic program is feasible, while violating the fewest possible inequalities for the system. Feasibility is assured by determining that the parameter test of Equation (14b) is satisfied, ensuring the solution to the strictly convex quadratic program for the MPC optimization problem of Equation (7) is feasible. Violation of the fewest possible inequalities is assured by determining that the active set a of problem solutions violates the fewest possible inequalities the parameter test of Equation (14a). Violations of the parameter test of Equation (14a) indicate that components of the Lagrange multiplier vector λ a  have the wrong sign and that a better suboptimal solution can be obtained by releasing each corresponding constraint from its boundary. The suboptimal v* and λ a  * can be calculated from Equation (13) for the identified active set a problem solution. 
   A predetermined time current solution can be calculated by solving the MPC optimization problem of Equation (6) or (7) online for the optimal solution with the parameters b s  and w 0  for the current operating state. To be useable as the alternative solution to the MPC problem, the calculation of the new solution can be completed within a predetermined time, such as within a sampling time. When the new solution cannot be completed within a predetermined time, a different alternative solution is used. However, the calculation for the full problem solution is also started and the corresponding optimal solution is inserted into the table whenever the calculation is completed. 
   A simplified model current solution can be calculated by solving a simplified version of the MPC optimization problem of Equation (6) or (7) online to find a simplified model current solution. The simplified version of the MPC optimization problem can have a shortened time horizon to find a shortened time horizon solution, additional constraints on the final state to find an additionally constrained solution, and/or other modeling simplifications to shorten the calculation time. Alternatively, a fall back solution can be obtained from the optimal solution computed at some previous decision time point. 
   The solution table can be managed to maintain a current selected subset of possible active sets a of problem solutions for the partial parameter region of the present operating condition. Maintaining a current selected subset reduces calculations at each sampling time since the optimal solution is likely to be found in the solution table. When the optimal solution is not in the solution table, a current solution to the model predictive control problem for the MIMO system for the current parameters can be calculated and added to the solution table, so that the solution is available at future sampling times. In one embodiment, the size of the solution table is allowed to increase. In another embodiment, an old problem solution can be deleted when the current solution is added to maintain the solution table at a predetermined length. Maintaining the solution table at a predetermined length assures that the scanning time does not become too large and that the storage space is conserved. The old problem solution that is deleted can be a least used solution or a least recently used solution. The least used solution can be determined by maintaining a use score for each of the problem solutions indicating how often the particular solution has been used. The least recently used solution can be determined by maintaining an age score for each of the problem solutions indicating how recently the particular solution has been used. 
   In one embodiment, the use score for each of the problem solutions can be determined over a predetermined number of sampling times and the scanning of the solution table for an optimal solution can proceed by scanning the problem solutions in order of the use score from most used to least used. In another embodiment, the age score for each of the problem solutions can be determined and the scanning of the solution table for an optimal solution can proceed by scanning the problem solutions in order of the age score from newest to oldest. In another embodiment when the optimal solution is not in the solution table, a solution family related to the current solution to the model predictive control problem for the MIMO system for the current parameters can be calculated and added to the solution table, so that the solution family is available at future sampling times. 
     FIG. 3  is a hardware diagram of a MIMO model predictive controller made in accordance with the present invention. The controller  300  includes a detector  302  providing a measured output from the system  304  to an online processor  306 , which supplies an input to the system  304 . The storage  308  stores problem solutions calculated by the offline processor  310  for use by the online processor  306 . 
   The detector  302  is any detector capable of measuring the output from the system  304 , such as fluid level, fluid flow, electrical voltage, electrical current, or any other measurable output. The offline processor  310  is any general purpose computer, dedicated computer, personal computer, or the like, including a computer readable medium storing a computer program for calculating the problem solutions as discussed for  FIG. 2  above. The storage  308  is any computer readable medium, such as read only memory, random access memory, a memory device, a compact disc, a floppy disk, or the like, able to store the problem solutions and computer readable code storing a computer program for controlling operations of the online processor  306 . The online processor  306  is any personal computer, microcontroller, microprocessor, or the like, programmed to calculate the optimal solution as discussed for  FIG. 2  above. In various embodiments, the detector  302 , online processor  306 , storage  308 , offline processor  310 , and/or any combinations or subcombinations thereof are imprinted on a single chip, such as an Application Specific Integrated Circuit (ASIC). 
   In an alternative embodiment, both the offline processor  310  and the online processor  306  are different operations on a single computer. In another embodiment, the offline processor  310  is operational to update, adjust, and optimize the problem solutions while the online processor  306  is operating. In yet another embodiment, the offline processor  310  is inactive or removed from the controller  300  when the online processor  306  is operating. 
   The method of model predictive control for a MIMO system is illustrated by the following two exemplary simulations: a co-polymerization reactor and a crude distillation unit. The simulations were run on a 1.2 GHz personal computer with a Debian Linux operating system using the Octave high-level interactive language for numerical computations, available at www.octave.org. The function cputime was used to measure the computational time required. The partial enumeration method results were compared to the qpsol method results for the same problems from a commercial active-set solver qpsol used to solve the quadratic programs. 
   The first exemplary simulation is for a co-polymerization reactor for methyl methacrylate and vinyl acetate. The simulation uses a discrete state-space model with eighteen states, five inputs, and four outputs, which is a medium sized system. The inputs are the methyl methacrylate flow rate, vinyl acetate flow rate, initiator flow rate, transfer agent flow rate, and jacket temperature. The outputs are polymer production rate, mole fraction methyl methacrylate in the polymer, average polymer molecular weight, and reactor temperature. The normalized constraints on the inputs are:
 
[−0.18 −0.18 −0.02 −0.02 −0.18 ]′≦u   k ≦[0.02 0.02 0.18 0.18 0.02]′.
 
   The nominal online model was modified with the addition of unmeasured random bounded input step disturbances and unmeasured random output noise to simulate an actual plant. An (A, B, C) controller design model was augmented with an input disturbance model, using B d =[B 1 , B 2 , B 3 , B 5 ] with B i  as the ith column of B, C d =0, and a steady-state Kalman filter designed for the augmented system. The model predictive controller was designed with Q=C′C, R=I, M=0, P as the corresponding Riccati penalty, and a time horizon of N=50. Under this qpsol model, the MPC optimization problem of Equation (7) has 250 decision variables, 500 constraints, and 3 250  (or 1.9×10 119 ) problem solutions over the total parameter region. 
   The partial enumeration method was simulated using 50, 100, and 200 problem solutions in the solution table and were labeled as TAB50, TAB100, and TAB200, respectively. The solution table was scanned in order of use score from most used to least used, with the use score being determined over a training period simulating 10 days of operation. When the optimal solution was not found in the solution table, a simplified model current solution with a shortened time horizon of N=5 was calculated as an alternative solution. The simulation runs were for 22 days (or 31681 sampling times at one sample per minute) with ten random input disturbances and random output noise. 
     FIGS. 4 &amp; 5  are simulation results showing the performance of a MIMO model predictive controller made in accordance with the present invention.  FIG. 4  is a graph of cumulative frequency versus number of entries scanned for the simulation results showing the performance of a MIMO model predictive controller made in accordance with the present invention. Results are shown for each of the cases TAB50, TAB100, and TAB200. The graph demonstrates that the partial enumeration method using a limited number of problem solutions in the partial parameter region is sufficient to provide an optimal solution in most cases. In this example as shown in  FIG. 4 , 80 percent of the sampling times scanned twenty problem solutions or fewer before finding the optimal solution. Thus, a solution table with 50 problem solutions was sufficient to provide an optimal solution for a great majority of the sampling times. The optimal solution was found 91.8 percent of the time in a solution table with 50 problem solutions and 94.2 percent of the time in a solution table with 200 problem solutions. 
   Table 1 summarizes the performance of the partial enumeration method and the qpsol method. The table provides the maximum and average (mean) CPU time required to compute the optimal solution or the alternative solution when the optimal solution is not in the solution table. For the partial enumeration method cases TAB50, TAB100, and TAB200, the table also provides the total number of different active sets of problem solutions scanned during the simulation and a suboptimality rating. The suboptimality rating is defined as the ratio of the achieved to optimal closed-loop objective function, with 1.0 indicating optimal. 
   
     
       
         
             
             
             
             
             
           
             
                 
               TABLE 1 
             
             
                 
                 
             
             
                 
               TAB50 
               TAB100 
               TAB200 
               qpsol 
             
             
                 
                 
             
           
          
             
                 
             
          
         
         
             
             
             
             
             
          
             
               Max time (sec.) 
               0.144 
               0.278 
               0.463 
               1.56 
             
             
               Mean time (sec.) 
               0.0075 
               0.0098 
               0.0129 
               1.02 
             
             
               Active sets visited 
               771 
               814 
               877 
               . . . 
             
             
               Suboptimality 
               1.0090 
               1.0088 
               1.0094 
               . . . 
             
             
                 
             
          
         
       
     
   
   From Table 1, the partial enumeration method typically runs 100 times faster than the qpsol method as indicated by the average CPU time, with negligible suboptimality. Even the slowest maximum time for the partial enumeration method is three times faster than the qpsol method. Use of the simplified model current solution with a shortened time horizon of N=5 as an alternative solution in about 8 percent of the sampling times had little effect on optimality. The suboptimality was less than 1 percent, while the suboptimality for the shortened time horizon case of N=5 used alone would be 5.2 percent. 
     FIG. 5  provides graphs of relative solution table differences for the TAB50 simulation versus time for the simulation results showing the performance of a MIMO model predictive controller made in accordance with the present invention. The relative solution table differences are defined as: 
   
     
       
         
           
             
               
                 
                   
                     
                       R 
                       d 
                     
                     ⁡ 
                     
                       ( 
                       
                         k 
                         , 
                         0 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       D 
                       ⁡ 
                       
                         ( 
                         
                           k 
                           , 
                           0 
                         
                         ) 
                       
                     
                     
                       N 
                       ⁡ 
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                 
                 , 
                 
                   
                     
                       R 
                       d 
                     
                     ⁡ 
                     
                       ( 
                       
                         k 
                         , 
                         
                           k 
                           - 
                           N 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       D 
                       ⁡ 
                       
                         ( 
                         
                           k 
                           , 
                           
                             k 
                             - 
                             N 
                           
                         
                         ) 
                       
                     
                     
                       N 
                       ⁡ 
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                 
               
             
             
               
                 ( 
                 15 
                 ) 
               
             
           
         
       
     
   
   where D(k,j) is the number of problem solutions in the solution table at sampling time k that were not in the solution table at sampling time j, N(k) is the total number of problem solutions in the solution table at sampling time k. In this simulation N(k) is constant because the solution table was previously populated at its maximum size. As indicated by the graph of R d (k,0), the initial 50 problem solutions in the solution table were the same for the first 30,000 minutes of the simulation, then the system moved into a different operating state so that new problem solutions were calculated and entered in the solution table. As indicated by the graph of R d (k, k−N), the number of problem solutions recently added to the solution table within the time horizon N increases and decreases when the system changes to different operating states. The relative solution table difference R d (k, k−N) returns to zero as the problem solutions in the solution table become suited to the different operating state, with different problem solutions acquiring higher use scores, new problem solutions being added to the solution table, and least used problem solutions being removed from the solution table. 
   The second exemplary simulation is for a crude distillation unit. The simulation uses a model with 252 states, 32 inputs, and 90 outputs, which is a large system. The inputs and outputs are normalized, with the inputs subject to the following constraints:
 
 u   min   ≦u   k   ≦u   max  
 
   with u max −u min =2·1 32 . Four outputs for the quality of the crude distillation unit side products have desired setpoints varying within [−0.1, 0.1]. 
   To simulate an actual plant, the nominal online model used by the model predictive controller was modified with the addition of unmeasured random bounded disturbances on crude composition, fuel gas quality, and steam header pressure. An (A, B, C) controller design model was augmented with a mixed input-output disturbance model (B d =[B 1 , B 2 , B 3 , B 5 ] with B i  as the ith column of B, C d ≠0) designed to assure detectability, and a steady-state Kalman filter designed for the augmented system. The model predictive controller was designed with Q=C′C, R=0.5 I, M=0, P as the corresponding Riccati penalty, and a time horizon of N=25. 
   Table 2 summarizes the results for a crude distillation unit. 
   
     
       
         
             
             
             
             
             
           
             
                 
               TABLE 2 
             
             
                 
                 
             
             
                 
               TAB50 
               TAB100 
               TAB200 
               qpsol 
             
             
                 
                 
             
           
          
             
                 
             
          
         
         
             
             
             
             
             
          
             
               Max time (sec.) 
               1.315 
               2.493 
               3.131 
               52.6 
             
             
               Mean time (sec.) 
               0.056 
               0.071 
               0.097 
               25.7 
             
             
               Suboptimality 
               1.0001 
               1.0001 
               1.0001 
             
             
                 
             
          
         
       
     
   
   While the embodiments of the invention disclosed herein are presently considered to be preferred, various changes and modifications can be made without departing from the spirit and scope of the invention. The scope of the invention is indicated in the appended claims, and all changes that come within the meaning and range of equivalents are intended to be embraced therein.