Extended horizon adaptive block predictive controller with an efficient prediction system
FIELD OF THE INVENTION 
     This application relates generally to the control and prediction of multivariable dynamic systems. The present invention is particularly applicable to real time automatic control systems and apparatus and more specifically to a class of adaptive controllers utilizing dynamic system prediction techniques and will be described with particular reference thereto, although it will be appreciated that the invention has broader applications. 
     BACKGROUND OF THE INVENTION 
     High performance, computationally efficient real-time dynamic controller software and hardware are required for use in complex multi-input, multi-output, nonlinear, time-varying systems that are operating in challenging environments. In particular, there is a need for efficient and reliable controllers for dynamic systems having large numbers of system inputs and outputs. The need for robust, adaptive controllers arises from the fact that systems engineers always encounter modeling uncertainties and unmeasurable external perturbations when designing control systems. It is generally impossible to describe the dynamics of most real systems by precise mathematical models. Although very good system models can be obtained for some devices in artificial laboratory environments, accurate modeling is unrealistic in complex and challenging environments. 
     Existing adaptive controllers such as the adaptive Pole Placement controller and the Minimum Variance controllers are promising approaches to accomplish these control tasks. However, they suffer from two potentially crippling limitations: (1) computational complexity, which limits their feasibility in multivariable applications; and (2) sensitivity to the choice of the input-output delays and model order selection. There has been considerable research recently in the development of adaptive controllers that attempt to overcome these limitations. A major focus has been the development of extended horizon predictive control methods. 
     PREVIOUS PREDICTIVE CONTROLLERS 
     The simplist predictive controllers, such as the Minimum Variance and the Generalized Minimum Variance methods, take into account the fact that the dynamic system has an input-output delay D. Control input is chosen to make the system match some desired trajectory in D steps ahead. Recently, predictive controllers that consider time horizons beyond the system input-output delay (extended horizon predictive controllers) have been used in engineering applications. The motivation for doing this is two-fold. First, the input-output delays of a dynamic system are usually not known in advance, and if the time delays are incorrectly estimated, or the delays are time varying as the system progresses, then the system input-output stability can suffer. Second, for high performance controller designs with fast sampling, the resulting sampled dynamic system often has nonminimum phase zeros (having zeros outside the unit circle). In this case, when the choice of control is only based on the beginning of the system step response, the controller often does not perform well. 
     Inside all extended horizon predictive and adaptive predictive controllers is a `predictor` that estimates future values of certain quantities related to the system outputs, based on current and past values of inputs and outputs. If the prediction horizon k extends beyond the system delay, D, then the set of (k-D) future inputs (u(t),u(t+1), . . . , u(t+k-D)) must be assumed. Different predictive and adaptive predictive controllers make different assumptions (or place different constraints) on these `extra` inputs. The predicted system outputs are then used to compute the control inputs, based upon some criterion. Several extended horizon predictive controllers have been developed previously. They differ in how the set of future controls is chosen. These controllers include: the Extended Horizon Adaptive Control (EHC), the Receding Horizon Adaptive Control (RHC), the Control Advance Moving Average Controller (CAMAC), the Extended Prediction Self-Adaptive Control (EPSAC), the Generalized Predictive Controller (GPC), the Model Predictive Heuristic Control (MPHC), and the Dynamic Matrix Control (DMC). 
     The EHC predictive controller developed by Ydstie is based on an Auto-Regressive Moving-Average with auxiliary input (ARMAX) model description of the dynamic system, as described in Goodwin and Sin. The set of k future controls is chosen by minimization of the control effort, subject to the constraint that the k-step ahead predicted output is equal to some desired value. In the EHC, two implementation approaches have been used: extended horizon and receding horizon implementations. In the extended horizon implementation, the set of k future controls is implemented sequentially and updated every k samples. In the receding horizon implementation of the EHC (i.e., the RHC), only the first control in the set is applied to the system. At the next sampling instant, the whole set of k future controls is computed and again only the first one is used. In the extended horizon implementation, the EHC can stabilize both open loop unstable and nonminimum phase systems. The RHC, however, cannot stablize an open loop unstable system in general. 
     The CAMAC controller developed by Voss et al, differs from the EHC controller in the selection of the set of future controls. In the CAMAC controller, all k future controls are assumed to be constant. This choice assures offset-free tracking performance in steady state, even without an integral action. The CAMAC controller can also be used for both receding horizon and extended horizon implementations. Again, the CAMAC controller can fail if the open loop system is unstable. Both the EHC and the CAMAC controllers have considered the output at only one future point in time. 
     The GPC controller developed by Clarke et al extends this idea by allowing the incorporation of multi-step output predictions. The set of k future controls in the GPC is determined by the minimization of a quadratic cost function of predicted output tracking errors and controls, up to a horizon value into the future. The GPC is implemented in a receding horizon manner, as in the RHC controller. The GPC controller is effective for both open loop unstable and nonminimum phase systems. Because a Controlled Auto-Regressive Integrated Moving-Average (CARIMA) model is used to model the dynamic systems and to predict the output, the GPC controller always contains an integrator. The GPC controller includes the EPSAC controller as a special case. 
     The MPHC controller by Richalet et. al. makes output predictions, based on an approximated impulse response model. On the other hand, the DMC controller by Cutler and Ramaker are based on an approximated step response model. Therefore, both MPHC and DMC controllers can be only applied to open loop stable systems. However, the DMC controller can stablize a nonminimum phase system. 
     The block controller of the subject invention differs from all above-mentioned controllers in that a block input and block output dynamic relationship and a, block state structure are used. The subject multi-input, multioutput block controller is built upon the single input and single output predictive controller developed in Dong and Chizeck. In particular, it is applicable for the control of multi-input and multi-output dynamic systems. It includes all existing predictive controllers as special cases. Each of the above controllers can be represented as a different special case of the subject block controller if certain values of the controller tuning parameters are enforced. 
     PREVIOUS OUTPUT PREDICTION METHODS 
     Prediction of the outputs of the dynamic system can be accomplished based on a system model description. For the dynamic systems that can be linearized using a linear model the general ARMAX model can be used to describe the input and output relationship of the dynamic systems: 
     
         A(q.sup.-1)y(t)=B(q.sup.-1)u(t-1)+C(q.sup.-1)w(t)          (1) 
    
     where u(t) is the dynamic system input of dimension M at time t: 
     
         y(t)=(y.sub.1 (t)y.sub.2 (t), . . . , u.sub.M (t)).sup.T   ( 2) 
    
     y (t) is the dynamic system output of dimension N at time t: 
     
         y(t)=(y.sub.1 (t)y.sub.2 (t), . . , y.sub.N (t)).sup.T     ( 3) 
    
     and w(t) is the noise term of dimension N at time t. It includes the effects of all unmeasured disturbance signals (stochastic noises as well as load disturbances) and of linearization offsets: 
     
         w(t)=(w.sub.1 (t)w.sub.2 (t), . . . ,w.sub.M (t)).sup.T    ( 4) 
    
     Here q -1  denotes the backward shift operator operating as g -  y(t)=y(t-1). This notation is customary for the description of real-time adaptive controllers, and is described in Goodwin, et al., Adaptive Filtering Prediction and Control, 1984 A(q -4 )={A if  (q -1 )}(i=1,2, . . . , N,j=1,2, . . , N) is an N by N polynomial matrix with polynomial entries A if  (q -1 ): 
     
         A.sub.if (q.sup.-1)=a.sub.0.sup.if +a.sub.1.sup.if q.sup.-1 +. . . +a.sub.n.sbsb.aif q.sup.-n.sbsp.aif                       ( 5) 
    
     and B(q -1 )={B if  (q -1 )}(i=1,2, . . . ,N,j=1,2, . . ,M) is N by M polynomial matrix with polynomial entries B if  (q -1 ): 
     
         B.sub.if (q.sup.-1)=b.sub.0.sup.if +b.sub.1.sup.if q.sup.-1 +b.sub.n.sbsb.bif.sup.if q.sup.-n.sbsp.bif                ( 6) 
    
     For the noise filter C(q -1 ), it is assumed that the noise components for the i th  loop can be combined into a single colored noise sequence C if  (q -1 )w i  (t) where for i=1,2, . . . ,N: 
     
         C.sub.if (q.sup.-1)=1+c.sub.1.sup.if q.sup.-1 +. . . +c.sub.n.sbsb.cif.sup.if q.sup.-n.sbsp.cli                ( 7) 
    
     If this is not the case, we can make it so by premultiplying the adjoint matrix of C (q -1 ), so as to convert it into an equivalent ARMAX model where C(q -1 ) is diagonal. 
     If the parameters in (1) are considered to be unknown parameters, then their estimates can be obtained on the basis of measured system inputs and outputs, and adjusted in response to the changes in the system dynamics, using a parameter identification method such as an recursive least square method as described in Ljung and Soderstrom. The estimated parameters (and thus the prediction system and controller operation) will, in general, be time varying. This allows for the tracking of time varying parameters and progressive refitting of a linear model approximation of a nonlinear system. 
     For the control of certain type of dynamic systems, such as lightly damped systems, it is often difficult to achieve the desired response to the command inputs on the basis of predictions of the system output alone. In such situations, predictor filters can be used to change the closed loop response properties. Another use of predictor filters is to reduce the interaction in multivariable control applications. The filtered output to be predicted can, in general, be described as: 
     
         z(t).tbd.T(q.sup.-1)y(t)                                   (8a) 
    
     where predictor filter T(q -1 )={T y  (g -1 )} for i=1,2, . . . ,N;j=1,2, . . . ,N where each entry of T can be factored into a ratio of two polynomials by ##EQU1## 
     Where it is assumed that polynomial R ij  (q -1 ) is monic, (i.e., the first coefficient is unity) and that it is stable. 
     Linear model representation of the dynamic system, in conjunction with an identification tool, can suffer if the effect of linearization offsets are significant. To reduce the effect of linearization offsets in system modeling, an integrator can be added to an ARMAX model, resulting a CARIMA model. Another benefit of the use of integral action with the ARMAX model is that this allows automatic removal of offsets when incorporated in a predictive controller. 
     Prediction of dynamic system outputs also involves a variety of data filtering and update operations. Data filtering in the output prediction serves two main purposes. First, it can be used to enhance the input and output signals by removing the noise and disturbances. Second, data filtering can be used to minimize the mean square output prediction errors. 
     There are well known time series analysis techniques for real time. prediction of ARMAX or CARIMA models, in which the coefficients of a k-step ahead predictor are computed sequentially, via the solutions of a matrix Diophantine equation. Several forms of such predictor systems are described in Favier and Dubois, Borsson, Koivo and Clark et al. In prior art applications, a typical Diophantine equation procedure requires the solution of two polynomial matrix Diophantine equations. Specifically, the first matrix Diophantine equation 
     
         det(C(q.sup.-1))=adj(C(q.sup.-1))A(q.sup.-1) F(q.sup.-1)+q.sup.-k G(q.sup.-1)                                               (9) 
    
     is used for the computation of the output data coefficient matrix G(q -1 ). The second matrix Diophantine equation 
     
         F(q.sup.-1)adj(C(q.sup.-1))B(q.sup.-1)=det(C(q.sup.-1))F&#39;(q.sup.-1)+q.sup.-k G&#39;(q.sup.-1)                                             (10) 
    
     is used to obtain the input data coefficient matrices F&#39;(q -1 ) and G&#39;(q -1 ). 
     The solution to equation (9) involves the determination of (N 3  ×n)+(N 2  ×k) scalar quantities from the N coupled polynomial Diophantine equations. The solution to equation (10) amounts to the solution of N coupled polynomial scalar Diophantine equations involving the determination of (N 2  ×M×n)+(N×M×k) scalar quantities. Thus a total of (N 3  ×n)+(N 2  ×M×n)+(N 2  ×k)+(N×M×k) scalar quantities must be determined for each k-step predictor. For the large dimensional dynamic systems, this computation is very inefficient because of the large memory storage for the coefficient matrices, and because of the computational intensity in the solutions of those coefficient matrices. This is particularly true for the system having a large number of outputs and a high model. parameter order. 
     The primary computational task in solving (9) and (10) are the calculation of the predictor input and output coefficient matrices. A vectored notation for equations (9) and (10) reveals that the solution of equation (9) requires approximately (2×N 3  ×n×k)+(2×N 3  ×n 2 ) floating point operations (flops) and that the solution of equation (10) requires roughly (2×N 2  ×M×n 2 )+(2×N 2  ×M×n×k) flops. Secondary tasks are to perform multiplications of the predictor input and output coefficient matrices with the associated input and output data vectors. Approximately 2×N 3  ×n flops are needed for multiplying the output coefficient matrices with the output data vectors and 2×N 2  ×M×n flops are needed for multiplying the input coefficient matrices with the input data vectors. Thus, the standard Diophantine equation method requires about 2×N 2  ×(N+M)×(n+k+1)×n flops for accomplishing the computation of a k-step predictor. For example, to compute the 10-step predictor for a modest system with 10 inputs and 10 outputs and a model order of 10, it will requires about 840 thousand flops. 
     For controller methods where a bank of predictors (i.e., k-step predictors for i=1,2, . . . ,k) are needed, the Diophantine equation method must solve 2k matrix Diophantine equations (with each i-step predictor being associated with two equations). Although the computational intensity can be somewhat reduced by using a recursive solution method for the Diophantine equations, the input and output coefficient matrices must still be computed for each i-step predictor. Thus, a total of at least (N 3  n+N 2  Mn+N 2 )k scalar quantities must be determined when a bank of predictors is needed. Thus, recursive solution methods for the Diophantine equations for the above mentioned example requires about 21 kilobytes of memory. Therefore, techniques to efficiently compute the predictor and to minimize the storage requirements, such that the computation of the large number of predictor coefficients can be avoided and hence memory storage can be reduced, as is accomplished by the subject invention, are crucial for real time prediction of dynamic systems. This includes real time control and adaptive control applications. 
     Although the output predictor reported in Dong and Chizeck makes use of blocks of inputs and outputs, the subject prediction system incorporates a unique construction of the block state structure which makes the subject prediction system computationally efficient. 
     Accordingly, it is an object of the present invention to provide a unified method and apparatus for the implementation of extended horizon predictive controllers. 
     It is another object of the present invention to provide a method and apparatus for the output prediction of those multivariable dynamic systems which can be approximated by a linear model. 
     It is another object of the present invention to provide an improved prediction method for multivariable dynamic system outputs. 
     It is a further object of the present invention to streamline the predictor implementation and to obtain predictions in a computationally efficient way. 
     It is also an object of the present invention to provide an efficient predictor system for adaptive multivariable predictive controllers. 
     These objects of the invention will become clear to those having the benefit of this disclosure. 
     SUMMARY OF THE INVENTION 
     The method of the present invention includes three embodiments. The first embodiment is the general block predictive controller that computes, updates and implements a predictive control law. The second embodiment is the efficient prediction system that generates the output predictions for multivariable dynamic systems that can be approximated as a time varying linear model. The third embodiment contains a detailed implementation of the efficient prediction system. The subject invention unifies existing extended horizon predictive controllers. The subject invention contains a general method for obtaining output predictions, utilizing a block state structure. The subject invention contemplates a new and improved prediction system for multivariable dynamic system output prediction. The subject invention reduces run time and memory requirements. The computational savings of the subject invention, in terms of floating point operations, compared to the Diophantine equation method, is cubic in the number of the system outputs N, quadratic in the dynamic system order n, and linear in both number of the system inputs M and horizon k. Specifically, let: 
     N=number of outputs in the multivariable system 
     M=number of inputs in the multivariable system 
     n=maximum number of {A, B, C} polynomial terms in the ARMAX model 
     k=prediction horizon 
     Then for calculation of a single, fixed step k-step ahead predictor: using standard Diophantine equation type methods, approximate N 3  n+N 2  Mn+N 2  k+NMk scalar quantities must be computed, requiring about 2N 2  (N+M)(n+k+1)n floating point operations. For the subject predictor system herein, only N 2  k+Nn+NMK quantities must be computed. This results in a savings of (N 3  +N 2  M-N)n fewer scalar quantities to be stored. The number of flops required is 2N 2  [(N+M)(n+k)]n. Thus using the subject predictor system, a savings of 2N 2  (N+M)n flops is obtained. In situations where a whole block of predictors or a bank of predictors is required (for horizons up to k), the savings are even larger. Compared to Diophantine methods, the savings is 2N 2  n[(N+M)n (k-1)+(N+M-1)k] flops (requiring Nn[KN 2  +NM-1] fewer scalar quantities). The savings in flops, for the predictor portion alone of the adaptive controller, is cubic in system output dimension N, quadratic in the ARMAX model order n, and linear in both input dimension M and horizon k. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The subject invention will be described in conjunction with a sell of figures wherein the illustrations are for the purpose of illustrating the preferred embodiment of the invention only. 
     FIG. 1 the diagram of the structure of the Adaptive Block Predictive Controller. 
     FIG. 2 is the diagram outlining the procedures of obtaining a prediction of blocks of dynamic system outputs 
     FIG. 3 is the diagram of the Block Control System. 
     FIG. 4 is the diagram of the preferred Block Output Predictor. 
     FIG. 5 is the diagram of the Data Filtering and Updating Unit. 
     FIG. 6 is the flow chart of the Predictor Input Data Filtering and Updating Section. 
     FIG. 7 is the flow chart of the Predictor Output Data Filtering and Updating Section. 
     FIG. 8 is the diagram of the Block State Processing Unit. 
     FIG. 9 is the flow chart of the Partial Output Block State Processing Section. 
     FIG. 10 is the flow chart of the Partial Input Block State Processing Section. 
     FIG. 11 is the flow chart of the ADS Mapping Processing. 
     FIG. 12 is the diagram of the Diagonal-Exclusive Outer Product Processing. 
     FIG. 13 is the diagram of the Block Predictor Processing Unit. 
     FIG. 14 is the description of computation engines for the Adaptive Block Predictive Controller. 
     FIG. 15 Embedded Controller Block Diagram 
     FIG. 16 Single Board Computer and DSP Function Block Diagram. 
     FIG. 17 is a plot of run time savings versus the number of system outputs. 
     FIG. 18 is a plot of run time savings versus the model parameter order. 
     FIG. 19 is a plot of run time savings versus the number of system inputs. 
     FIG. 20 a plot of run time savings versus the prediction horizon of the controller. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS: 
     Now, embodiments of the present invention will be described in detail with reference to the drawings. 
     FIG. 1 illustrates the overall structure of extended horizon adaptive block predictive controller with an efficient prediction system. A block of potential system inputs (block controls) (u(t),u(t+1), . . . ,u(t+k-1)) is determined by Block Control System 1, based upon a block of the desired output trajectory (y*(t+1), y*(t+2), . . . , y*(t+k)) and i) a block of predicted outputs (y p  (t+1),y p  (t+2), . . . ,y p  (t+k)) when no filtering of the system outputs is used, or a block of predicted filtered outputs (z p  (t+1),z p  (t+2), . . . ,z p  (t+k)) when filtering of the system outputs is used, both being generated by Block Predictor System 2. The block of predicted outputs are based upon measurements of the system outputs of the Dynamic System 3, and either dynamic system parameters from 3, if known or identified parameters from the Identifier 4. Both the inputs and outputs of the Dynamic System 3 may be corrupted by noise. The Block Predictive Controller 5 consists of the Block Control System 1, the Block Predictor System 2, the Control Implementation 6 and the Control Update 7. The Block Predictive Controller 5 together with Identifier 4, comprises the Adaptive Block Predictive Controller 8. This Adaptive Block Predictive Controller 8 resides within a computational engine described later. For predictive control with horizon k, the Block Predictor System 2 must generate a one of the predicted outputs (y p  (t+1),y p  (t+2), . . . ,y p  (t+k)) or predicted filtered outputs (z p  (t+1),z p  (t+2), . . . ,z p  (t+k)) in order to choose block controls. If k=D, then u(t) is chosen at time t. If k≧D, then additional inputs u(t+1), . . . ,u(t+k-D) must also be chosen at time t. If D is not known, then u(t+1), . . . ,u(t+k-1) can be chosen. Two issues in this adaptive block predictive controller configuration that concern how measured information and parameter estimates are used, are the update method and implementation method. Control Implementation 6 refers to the implementation method. The whole block controls (u(t),u(t+1), . . . ,u(t+k-1)) that is computed by the Block Control System 1 at time t can be used. Then new block controls are obtained k steps later. This results in an extended horizon implementation. An alternative is to have only the first element in the block controls, u(t), actually delivered to the system. The other calculated elements of the block controls at time t are not used. Then at time t+1, the calculation process is repeated for new block controls, and only the first element is used. This is a receding horizon implementation method. It is also possible to use the first i&lt;k of the inputs. This is the i-step implementation method. The full extended horizon and the receding horizon implementation methods have been used by prior art predictive controllers. The i-step implementation method described here is novel. Control Update 7 refers to the update method. The calculation of the output predictions in the Block Predictor System 2 depends upon measured outputs, y, and parameter estimates. In block-by-block updates, new parameter estimates and output measurements are provided only every k time steps. In element-by-element updates, the output measurements and parameter estimates are updated at each time step. 
     FIG. 2 further describes the Block Predictor System 2 of FIG. 1. It comprises a Block State Structure Construction 10 and a Block Output Predictor 20. The Block Output Predictor 20 generates a block of k output predictions (y p  (t+1),y p  (t+2), . . ., y p  (t+k)) (when the filtering is used) or a block of k filtered output predictions (z p  (t+1),z p  (t+2), . . . ,z p  (t+k)) (when filtering is used), where k is the prediction horizon. The system inputs from the Control Implementation 6, the system outputs from the Dynamic System 3, and the parameters of the Dynamic System 5 when the parameters are known, or the parameter estimates from the Identifier 4 when the parameters of the Dynamic System 3 are unknown (adaptive case), are used for the Block State Structure Construction 10. This requires five steps. 
     Construction of the N by N polynomial matrix E(q -1 ) at time t as follows: ##EQU2## with polynomial R j  (q -1 ) for j=1,2, . . . ,N being the least common multiples of polynomials R if  (q -1 ) of (8b) for i=1, 2, . . . ,N. 
     Construction of the filtered outputs ψ(t) at time t as: ##EQU3## where polynomial matrix R(q -1 ) of dimension N by N is defined to be: 
     
         R(q.sup.-1)=diag{R.sub.1 (q.sup.-1),R.sub.2 (q.sup.-1), . . . , R.sub.N (q.sup.-1)}                                               (2.4) 
    
     Note that the R if  &#39;s in (8b) are not the ij th  element of this R(q -1 ). Rather, R(q -1 ) in (2.4) is defined from the R j  (q 1 ) of (2.2) which are the least common multiples of polynomials R if  (q -1 ) of (8b). Of course, setting (Rcq -1 ) equal to the identify matrix results in ψ(t) to be exactly equal to y(t) and hence no filtering is performed, 
     Construction of the Input Block U(t) of dimension Mk by 1, the Filtered Output Block Ψ(t) of dimension Nk by 1 and Noise Block W(t) of dimension Nk by 1 at time t as follows: 
     
         Ψ(t)=(ψ.sub.1 (t),Ψ.sub.2 (t), . . . , Ψ.sub.N (t)).sup.T(2.5) 
    
     
         U(t)=(U.sub.1 (t),U.sub.2 (t), . . . , U.sub.M (t)).sup.T  (2.6) 
    
     
         W(t)=(W.sub.1 (t),W.sub.2 (t), . . . , W.sub.N (t)).sup.T  (2.7) 
    
     
         Ψ.sub.i (t)=(ψ.sub.i (t+k),ψ.sub.i (t+k-1), . . . , ψ.sub.i (t+1)).sup.T                                              (2.8) 
    
     
         U.sub.i (t)=(u.sub.il (t+k-1),u.sub.i (t+k-2), . . . , u.sub.i (t)).sup.T(2.9) 
    
     
         W.sub.i (t)=(w.sub.i (t+k),w.sub.i (t+k-1), . . . , w.sub.i (t+1)).sup.T(2.10) 
    
     In (2.5), ψ i  (t) for i=1,2, . . . , N is the i th  filtered output block of size k. In (2.6), U,(t) for i=1, 2, . . . , M is the i th  input block of size k. In (2.3), W i  (t) for i=1, 2, . . . , N is the i th  loop noise block of size k. 
     Construction of the Block System that relates nonoverlapping blocks of ψ(t), U(t) and W(t) in terns of matrix polynomials E(d -1 ), B(d -1 ) and C(d -1 ) in the k-fold backward shift operator d -1  defined as q -k  : 
     
         E(d.sup.-1)ψ(t)=B(d.sup.-1)U(t)+C(d.sup.-1)W(t)        (2.11) 
    
     
         E(d.sup.-1)={E.sub.i,j (d.sup.-1)) (i=1,2, . . . ,N;j=1,2, . . . ,N)(2.12) 
    
     
         B(d.sup.-1)={B.sub.i,j (d.sup.-1 }(i=1,2, . . . , N;j=1,2, . . . , M)(2.13) 
    
     
         C(d.sup.-1)=diag{C.sub.11 (d.sup.-1),C.sub.22 (d.sup.-1), . . . , C.sub.NN (d.sup.-1)                                                (2.14) 
    
     Here the polynomial submatrices E if  (d -1 ), b if  (d -1 )and C if  (d -1 ) are all constructed in a special circulant structure, known as the dynamic shift Toeplitz matrix structure. For example, the submatrix E if  (d -1  ) for k=3 case can be written as: ##EQU4## where E i   if  (d -1 ) for l=1,2,3 is a polynomial in d -1 . ##EQU5## Here n is the parameter order defined as max ##EQU6## The notation ceiling in (2.17) denotes the largest integer in the fraction of n/k. 
     Note that the matrix polynomial E if  (d -1 ) is completely determined by its first row. All of the submatrices in the above block system representation have this special structure. They are constant along all diagonals, and have a factor of d -1  for all elements below the principal diagonal. Matrices of this circulant (&#34;shift Toeplitz&#34;) type are commutative. transformation of the Block System (2.5), (2.6) and (2.7) into the following Block State Structure using the commutivity property of dynamic shift Toeplitz matrices. 
     
         ψ(t)=MARKOV×U(t)+REC×X(t)+W(t)             (2.18) 
    
     where 
     •MARKOV is the Nk by Mk Impulse Response Matrix constructed as: 
     
         MARKOV={MARKOV.sub.if }(i=1, . . . ,N;j=1, . . . ,M)       (2.19) 
    
     Here submatrix MARKOV if  is a k by k upper-triangular Toeplitz matrix with vector (g l   if , . . . ,g k   if ) being its first row and {g l   if  } l=1   k  are the first k impulse response parameters from the input A us(t) (defined as u i  (t)-u i  (t-1)) to the output ψ i  (t) 
     •REC is the Nk by nk Block State Mapping Matrix constructed as: 
     
         REC={REC.sub.if }(i=1, . . . ,N;j=1, . . . ,N)             (2.20) 
    
     Here submatrix REC u  is a k by n matrix defined as: ##EQU7## Here s l   if  for l&gt;0 is defined by the inverse of autoregressive polynomial matrix E (q -1 ): ##EQU8## with s l   if  =0 for l&lt;0 or, equivalently, by the inverse of the lowest degree coefficient matrix in E(d -1 ). 
     •X(t) is the Block State of dimension N by n, constructed as: ##EQU9## 
     Here Γ if  and His-are all anti-diagonal symmetric matrices of dimension n by n, defined as 
     
         Γ.sup.if =Δ.sub.L (c.sub.l.sup.if, . . . , c.sub.n.sup.if)Δ.sup.U (e.sub.0.sup.if, . . ., e.sub.n-1.sup.if) -Δ.sub.L (e.sub.l.sup.if, . . . ,e.sub.n.sup.if)Δ.sup.U (c.sub.0.sup.if, . . . , c.sub.n-1.sup.if)                (2.24) 
    
     
         II.sup.if =Δ.sub.L (b.sub.l.sup.if, . . . ,b.sub.n.sup.if)Δ.sup.U (C.sub.0.sup.if, . . . ,c.sub.n-1.sup.if) -Δ.sub.L (c.sub.l.sup.if, . . . , c.sub.n.sup.if)Δ.sup.U (b.sub.0.sup.if, . . . ,b.sub.n-1.sup.if)                 (2.25) 
    
     with Δ L , (a 1 , . . . ,a n ) denoting an n by n lower triangular Toeplitz :matrix with (a 1 ,. . . ,a n ) as its first row and Δ U  (a l , . . . ,a n ) denoting an n by n upper triangular Toeplitz matrix with (a l , . . . ,a n ) as its last row. 
     •Y f  (t) of dimension N by N by n is the Predictor Output Data at time t, constructed as follows: 
     
         Y.sup.f (t)=(Y.sub.1.sup.f (f),Y.sub.2.sup.f (t), . . . , Y.sub.N.sup.f (t))(2.26) 
    
     Here submatrix Y i   f  (t) for i=1,2, . . . ,N of dimension N by n contains the filtered i th  output: 
     
         Y.sub.i.sup.f (t)=(Y.sub.i.sup.f1 (t),Y.sub.i.sup.f2 (t), . . . , Y.sub.i.sup.fN (t))                                       (2.27) 
    
     The vectors Y i   fj  (t) for j=1,2, . . . ,N are obtained as follows: ##EQU10## 
     •U f  (t) of dimension N by M by n is the Predictor Input Data at time t: 
     
         U.sup.f (t)=(U.sub.1.sup.f (t),U.sub.2.sup.f (t), . . . ,U.sub.M.sup.f (t))(2.30) 
    
     Here submatrix U i   f  (t) for i=1,2, . . . ,M is of dimension N by n. It contains the filtered i th  input: 
     
         U.sub.i.sup.f (t)=(U.sub.i.sup.f1 (t),U.sub.i.sup.f2 (t), . . . , U.sub.i.sup.fN (t))                                       (2.31) 
    
     The vectors U i   fj  (t) for j=1,2, . . . ,N are obtained as follows: ##EQU11## 
     Block Output Prediction 20 uses the block state structure of (2.18) to generates the block predictor Z p  (t) as follows. First, block predictor Ψ p  (t) is obtained by minimizing the conditional expected mean squared prediction errors E{||Ψ(t)-Ψ p  (t)|| 2  }, conditioned on all inputs and outputs up to and including time t: 
     
         Ψ.sub.p (t)=MARKOV×U(t)30 REC×X(t)         (2.34) 
    
     Note that the block filtered system output prediction Ψ p  (t) is the product of the Impulse Response Matrix MARKOV with the Block Control U(t) in summation with the product of the Block State Mapping Matrix-REC with the Block State X(t), which is independent of the prediction horizon k. 
     The second step in the method of block 20 is to filter Ψ P  (t) using the Nk by Nk polynomial matrix filter P(d -1 ): 
     
         Z.sub.P (t)=P(d.sup.-1)Ψ.sub.P (t)                     (2.35) 
    
     with d -1  Ψ P  (t)=Ψ(t-i) for i&lt;0, where 
     
         P(d.sup.-1)={P.sub.if (d.sup.-1)}(i=1,2, . . . ,N; j=1,2, . . . , N)(2.36) 
    
     Here the polynomial submatrix P if  (d -1 ) is a dynamic shift Toeplitz matrix with its first row given by (P 1   if  (d -1 ), P 2   if  (d -1 ), . . . , P k   if  (d -1 )) (since this dynamic, shift Toeplitz matrix is completely determined by its first row) where P i   if  (d -1 ) for l=1,2, . . . , is a polynomial in d -1  defined as: ##EQU12## The block predictor Z p  (t) obtained by (2.35) can be also derived by directly minimizing the cost function of E{||Z(t)-Z p  (t)||} (conditioned by all past inputs and outputs up to and including time t). 
     FIG. 3 further describes the Block Control System 1. The Block Control Law 100 is described by the following equation: 
     
         R(d.sup.-1)U(t)=-Q(d.sup.-1)Z.sub.P (t)+M(d.sup.-1)Y.sub.r (t)(3.1) 
    
     where Y r  (t) is a given block desired output defined as: 
     
         Y.sub.r (t)=(Y.sub.1r.sup.T (t),Y.sub.2r.sup.T (t), . . . ,Y.sub.Nr.sup.T (t)).sup.T                                                (3.2) 
    
     with 
     
         Y.sub.ir =(y.sub.i.sup.o (t+1),y.sub.i.sup.o (t+2), . . . ,y.sub.i.sup.o (t+k))                                                    (3.3) 
    
     •Z p  (t) is the block output (or block filtered output) prediction at time t with d -1  Zp(t)=Z(t-i) for i&gt;0. 
     •R(d -1 )={R if  (d -1 )}(i=1,2, . . . ,N;j=1,2, . . . ,N)is the kM by kM polynomial matrix with polynomial submatrix R if  (d -1 ) being defined by 
     
         R.sub.if (d.sup.-1)=R.sub.0.sup.if +R.sub.1.sup.if d.sup.-1 +. . . +R.sub.n.sup.if, d.sup.-n,                                (3.4) 
    
     Here R 1   if  (l=1,2, . . . n,) is a k by k constant matrix. 
     •Q(d -1 )={Q if  (d -1 )}(i=1,2,. . . , M;j=1,2,. . . N)is the kM by kN polynomial matrix with polynomial submatrix Q if  (d -1 ) being defined by 
     
         Q.sub.if (d.sup.-1)=Q.sub.0.sup.if +Q.sub.1.sup.if d.sup.-1 +. . . +Q.sub.n.sbsb.q.sup.if d.sup.-n.sbsp.q                    (3.5) 
    
     Here Q 1   if  (l=1,2, . . . ,n) is a k by k constant matrix. 
     •and M(d -1 )={M if  (d -1 )}(i=1,2, . . . ,M;j=1,2, . . . ,N)is the kM by kN polynomial matrix with polynomial submatrix M if  (d -1 ) being defined by 
     
         M.sub.if (d.sup.-1)=M.sub.0.sup.if +M.sub.1.sup.if d.sup.-1 +. . . +M.sub.n.sbsb.q.sup.if d.sup.-n.sbsp.m                    (3.6) 
    
     Here M 1   if  (l=1,2, . . . , n m ) is a k by k constant matrix. 
     This block control equation subsumes most existing extended horizon adaptive predictive controllers including the MHPC controllers, the DMC controller, the EHC controller, the CAMAC controller, the GPC controller, and the EPSAC controller. In Dong and Chizeck, several of the above-mentioned controllers are translated into a single input, single output block structure. The subject block controller described here is more general. In particular, it is applicable to the multi-input and multi-output case. The multivariable version of the GPC controller is given below. This is a typical example of the most commonly used controllers. Let k 1 ,k 2 , . . . ,k N  be the prediction horizons corresponding to system outputs y 1  (t),y 2  (t), . . . ,y n  (t), respectively. Let nu 1 ,nu 2  , . . . ,nu M  be the control horizons corresponding to system inputs u 1  (t),u 2  (t), . . . ,u M  (t), respectively. Let k=max{k 1 ,k 2 , . . . ,k N  }. Let minimum prediction horizons for outputs y 1  (t),y 2  (t), . . . ,y N  (t) be l 1 ,l 2 , . . . ,l N , respectively. The multivariable GPC controller can be derived from the subject invention controller with the following choices of R(d -1 ), Q(d -1 ) and M(d -1 ): 
     
         R(d.sup.-1)=Diag(λ.sub.1 D(d.sup.-1), λ.sub.2 D(d.sup.-1), . . . , λ.sub.M D(d.sup.-1)}                         (3.7) 
    
     where D(d -1 ) is a k by k dynamic shift Toeplitz matrix with the first row given by (1,-1, 0, . . . , 0) and λ i  is a nonnegative control weighting factor for the i th  input u i  (t) for i=1,2, . . . ,M. Q(d -1 ) and M(d 4 ) are chosen as: 
     
         Q(d.sup.-1)=M(d.sup.-1)=(NU)(MARKOV.sup.T)(P.sub.0.sup.T)(Λ)(BPS)(3.8) 
    
     Here NU is a Mk by Mk block diagonal matrix: 
     
         NU=diag{NU.sub.1, NU.sub.2, . . . , NU.sub.M }             (3.9) 
    
     with diagonal submatrix NU i  for i=1,2, . . . ,M being ##EQU13## MARKOV is the Nk by Mk Impulse Response Matrix: 
     
         MARKOV={MARKOV.sub.if }(i=1, . . . ,N;j=1, . . . . ,M)     (3.11) 
    
     Here submatrix MARKOV if  is a k by k upper-triangular Toeplitz matrix with vector (g 1   if ,. . . . ,g k   if ) being its first row and {g 1   if  } l=1   k  , are the first k impulse response parameters from the input Δu j  (t) to the output ψ i  (t) defined in (2.3). 
     P 0  is Nk by Nk constant matrix: 
     
         P.sub.0 ={P.sub.0.sup.if }(i=1,2, . . . , N;j=1,2, l . . . , N)(3.12) 
    
     with submatrix P 0   if  is the k by k coefficient matrix of the lowest degree term of the dynamic shift Toeplitz matrix P if  (d -1 ) of dimension k by k, with its first row given by (P 1   if  (d -1 ), P 2   if  (d -1 ),. . . , P k   if  (d -1 )) (since a dynamic shift Toeplitz matrix is completely determined by its first row) where P l   if  (d -1 ) for l=1,2, . . . ,k is a polynomial in d -1  defined as: 
     
         P.sub.1.sup.if (d.sup.-1)=P.sub.l-1.sup.if +P.sub.l+k-1.sup.if d.sup.-1 +. . . +P.sub.l+kh-1.sup.if d.sup.-h (l=1,2,. . . , k)       (3.13) 
    
     Here ##EQU14## and A=diag{A 1 ,A 2 , . . . ,A N  } is a Ark by Nk block diagonal matrix with submatrix A i  (i=1,2, . . . ,N) of dimension k by k being a diagonal matrix defined by 
     
         Λ.sub.i =ρ.sub.i I                              (3.16) 
    
     where ρ i  for i=1,2, . . . ,N is a posive scalar and I is the k by k identity matrix. BPS is the Nk by Ark matrix defined as 
     
         BPS=diag {BPS.sub.1, BPS.sub.2, . . . , BPS.sub.N }        (3.17) 
    
     where BPS, (i=1,2, . . . ,N) is the k by k matrix: ##EQU15## Since the GPC controller includes the EPSAC controller as a special case, the above construction for the GPC controller also includes the EPSAC controller. 
     FIG. 4 is the block diagram of the preferred block output predictor. It consists of three main units: the Data Filtering and Updating Unit A, the Predictor State Processing Unit B, the Block Predictor Processing Unit C. The preferred predictor system is organized to process the output y(t) or the filtered output z(t) predictions as follows. First, the Data Filtering and Updating Unit A filters inputs and outputs obtained from the dynamic system 3 at time t and then uses these filtered inputs and outputs to update the Predictor Input Data and the Predictor Output Data. Second, Unit B computes 2Block State X(t) using the updated Predictor Input Data and Predictor Output Data from Unit A. Next, the Block Predictor Processing Unit C computes the block output predictions Y p  (t) or the block filtered output predictions Ψ p  (t), recursively in horizon, using the Block State from Unit B. Finally, Unit C obtains the desired output predictions Z p  (t) by filtering Ψ p  (t) through the matrix polynomial filter P(q -1 ). The detailed operations of the preferred predictor system are described below in conjunction with FIG. 5 to FIG. 13. 
     FIG. 5 is the diagram of the Data Filtering and Updating Unit A. It contains two sections. Section A100 performs input filtering and updates the Predictor Input Data, U f , from U f  (t-1) at time t-1 to U f  (t) of the present time, using newly obtained filtered system inputs. Here the t is the current time index. Section A200 performs output filtering and updates the Predictor Output Data, Y f , from Y f  (t-1) at time t-1 to Y f  (t) at the present time t, using newly obtained filtered system outputs. The Predictor Input Data U f  (t) and Predictor Output Data Y f  (t) are then used to compute the Block State X(t) in Unit B. The Predictor Input Data and the Predictor Output Data need to be initialized at time t=0 and thereafter are filtered and updated sequentially. 
     FIG. 6 is the flow chart for the Predictor Input Data Filtering and Updating Section. The Predictor Input Data in A101 is constructed as a 3-dimensional array. At time t-1 (before being updated), U f  can be written as U f  (t-1)=(U 1   f  (t-1), U 2   f  (t-1), . . . ,U M   f  (t-1)) where each U i   f  (t-1) for i=1,2, . . . ,M is a 2-dimensional array comprising the following vectors U i   f  (t-1)=(U i   fi  (t-1),U i   f2  (t-1), . . . ,U i   fN  (t-1)). The vector U i   fjp  (t-1) can be further expressed component-wise as U i   fj  (t-1)=(u i   fj  (t-1), u i   fj  (t-2), . . . , u i   fj  (t-n)) T  for the value of j=1,2, . . . ,N that contains n most recently filtered i th  input. The Predictor Input Data filtering computes the filtered i th  input u i   fj  (t) for i=1,2, . . . ,M and j=1,2, . . . ,N in the following two steps. First, A102 performs a vector dot product of vector U i   fj  (t-1) from A101 and filter polynomial C if  (q -1 ) from A103. The resulting value is saved as S in A104. Second, in A102 the filtered input u i   fj  (t) is computed as: u i   fj  (t)←u i  (t)-u i  (t-1)-S. A105 performs updating operations for the Predictor Input Data U f , from U f  (t-1) to U f  (t), by updating each individual vector U i   fj  (t-1) i=1,2, . . . ,M and j=1,2, . . . , N. The updated Predictor Input Data U f  (t) is then used for the block state processing. 
     The Predictor Output Data filtering and updating operations are performed similarly to those of the Predictor Input Data filtering and updating. FIG. 7 is the flow chart for the Predictor Output Data Filtering and Updating Section. The Predictor Output Data in A201 is constructed as a 3-dimensional array. At time t-1 (before being updated), Y f  can be written as Y f  (t-1)=(Y 1   f  (t-1), Y 2   f  (t-1), . . . ,Y M   f  (t-1)) where each Y 1   f  (t-1) for the value of i=1,2, . . . ,N is a 2-dimensional array comprising the following vectors Y i   f  (t-1)=(Y i   f1  (t-1), Y i   f2  (t-1), . . . , Y i   fN  (t-1)). The vector Y i   fj  (t-1) can be further expressed component-wise as Y i   fj  (t-1)=(y i   fj  (t-1),y i   fj  (t=2), . . . , y i .sup. fi (t-n)) T  for the j=1,2, . . . , N that contains n most recently filtered i th  output. The Predictor Output Data filtering operations compute the filtered i th  output y i   fj  (t) for i=1,2, . . . ,N and j=1,2, . . . . ,N in the following two steps. First, A202 performs a vector dot product of vector Y i   fj  (t-1) from A201 and filter polynomial C jj  (q -1 )R j  (q -1 ) from A203. The resulting dot product value is saved in S of A204. Second, in A202 the filtered output y i   fj  (t) is computed as: y i   fj  (t)ιy i  (t)-S. A105 performs updating operations for the Predictor Output Data Y f  from Y f  (t-1) of t-1 to Y f  (t) of the present time by updating each individual vector Y i   fj  (t-1) i=1,2, . . . ,N and j=1,2, . . . , N. The updated Predictor Output Data Y.sup. f (t) is then used for the block state processing. 
     FIG. 8 outlines the processes of the Block State Processing Unit B. Referring to the figure, the Block State Processing Unit contains the Partial Output State Processing Section B100 and the Partial Input State Processing Section B200. The Block State Processing Unit B determines the block state at time t as the sum of the Partial Input State and Partial Output State. The Block State is used in Block Predictor Processing Unit C. The detailed operations of the Block State Processing Unit are described below with reference to FIG. 9 to FIG. 14. 
     FIG. 9 is the flow chart for the method of computing the (N×n) 2-dimensional Partial Output State array X y  (t)=(X y  (1,t),X y  (2,t), . . . ,X y  (N,t)) T . Each vector X y  (i,t) for i=1,2, . . . . ,N is of length n. Referring to the figure, the computation of the Partial Output State Xy(t) is accomplished as follows. For each i, the vector Xy(i,O is first initialized at zero. The ADS Mapping Processing B101 then performs an anti-diagonal symmetric (ADS) mapping process for polynomials C jj  (q -1 ) and E if  (q -1 ), in order to obtain the n×n anti-diagonal symmetric matrix Γ ij . The detailed operations of the ADS Mapping Processing will be described below. Next, the Matrix-Vector Product B102 carries out the multiplication of matrix Γ if  with vector Y i   fj  (t) and saves the resulting vector as S. Finally, the Vector Summation B103 performs the vector accumulated summation: Xy(i,t)=Xy(i,O+S for j=1,2, . . . . ,N. 
     FIG. 10 is the flow chart for the method of computing the (N×n) 2-dimensional Partial Input State X u  (t)=(X u  (1,t),X u  (2,t), . . . ,X u  (N,t)) r  with each X u  (i,t) being a vector of length n for i=1,2, . . . . ,N. The procedures for the computation of the Partial Input State are analogous to those for the computation of the Partial Output State. Referring to the figure, the computation of the Partial Input State Xu(t) is accomplished as follows: for each i, the vector Xu(i,t) is first initialized at zero. The ADS Mapping B201 then performs an anti-diagonal symmetric mapping operation for polynomials B ji  (q -1 ) and C ii  (q -1 ), in order to obtain the n×n anti-diagonal symmetric array II if . Next, the Matrix-Vector Product B202 carries out the multiplication of matrix II if  with vector U i   fi  (t) and saves the resulting vector as S. Finally, the Vector Summation B203 performs the vector accumulated summation: Xu(i,t)=Xu(i,t)+S for j=1,2, . . . . ,M. 
     FIG. 11 is the diagram of the detailed construction and operation of ADS mapping of two polynomials. Referring to the figure, Block ADS300 performs the diagonal-exclusive outer product on polynomials X of degree n in ADS100 and Y of degree n in ADS200. As shown in FIG. 12, the diagonal-exclusive out-product matrix resulting from the operations of ADS300 is divided into two upper left triangular matrices P1 as in ADS301 and P2 as in ADS302. Referring back to FIG. 11, ADS400 determines the upper left triangular part of the matrix in the following steps. First, P2 is subtracted P1. Second, P1 is updated according to: P1(i+j+1,j+1)←P1(i+j+1,j+1) +P1(i,O) for j≦(n-i-2)/2 and i=1,2, . . . ,n. Finally, elements in the lower left triangular part of P1 are assigned with values determined as the followings: P1(i+j,j+1)←P(i,j) for j&gt;(n-i-2)/2 for i=1,2, . . . ,n. The updated P1 resulting from the three steps gives the ADS matrix of polynomials x and y. Note that the elements in the lower anti-diagonal part need not be computed because P1 is an anti-diagonal symmetric matrix. 
     FIG. 13 is the block diagram of the Block Predictor Processing unit which computes the block predictor Ψ p  (t) recursively in the prediction horizon. Since the block state X(t) is invariant with the horizon k and since the block of future controls is assumed to be known in advance, the recursive computation of the Ψ P  (t) amounts to the recursire computation of the Impulse Response Matrix MARKOV and the Block State Mapping Matrix REC. C100 first constructs N vectors H 1   i  of dimension 1 by Nn for i=1,2, . . . ,N obtained by concatenation of the l th  rows of submatrices REC il  to REC iN  (in order). Here l runs from to k. When l increases from l to l+1, only n new values need to be computed for each vector H l+1   i , for i=1,2, . . . ,N, because all submatrices in E(d -1 ) are dynamic shift Toeplitz matrices. The relationship between vector H l   i  and H l+1   i  can be described as: 
     
         H.sub.l+1.sup.i =H.sub.l.sup.i C.sub.e                     (13.1) 
    
     where C e  ={C e   if  } for (i=1,2, . . . , N;j=1,2, . . . ,N) with the submatrices C e   if  being n by n matrices defined as ##EQU16## C200 computes the impulse response parameter g l   if  for i=1,2, . . . . ,N and j=1,2, . . . ,M via vectors H l-1   i  as 
     
         g.sub.l.sup.if =H.sub.l-1.sup.i V.sub.i                    (13.4) 
    
     where V j  and B j  for j=1,2, . . . . ,M are the Nn by 1 vectors: ##EQU17## Here B if   c  (for i=1,2, . . . . ,N and j=1,2, . . . ,M) denotes the n by 1 vector that contains all of the coefficients of the polynomial B if  (q -1 ) except b o   if . Similarly E y   c  (for i=1,2, . . . . ,N and j=1,2, . . . ,N) denotes the n by 1 vector that contains all the coefficients of the polynomial E if  (q -1 ) except e o   if  : 
     
         B.sub.if.sup.c =(b.sub.n.sup.if,b.sub.n-1.sup.if, . . . ,b.sub.l.sup.if,).sup.T                                   (13.7) 
    
     
         E.sub.if.sup.c =(e.sub.n.sup.if,e.sub.n-1.sup.if, . . . ,e.sub.i.sup.if,).sup.T                                   (13.8) 
    
     Thus, C300 determines/-step ahead predictor ψ p   i  (t+l) for the i th  output y i  (t) for i=1,2, . . . . ,N as ##EQU18## Equations (13.1), (13.4) and (13.9) provide a recursive-in-horizon method for the computation of block predictor Ψ P  (t) for increasing horizon l, from 1 through k, with the initial conditions for H 1   i  for i=1,2, . . . ,N given by 
     
         H.sub.l.sup.i =(REC.sub.1.sup.1i,REC.sub.1.sup.2i, . . . , REC.sub.1.sup.Ni)(13.10) 
    
     with REC 1   if  =(0,0, . . . 1) if i=j and otherwise REC j   if  =(0,0, . . . 0). 
     Finally, C400 performs filtering operation on Ψ p  (t) to obtain the block predictor Z p  (t). 
     FIG. 14 is the description of the computational engine which is used in conjunction with the extended horizon adaptive block predictive controller and efficient predictor system. The computational engine is comprised of computing means 50 and memory means 60. 
     When required due to the nature of the dynamic system input or output signals, dynamic system signal output processing means 70 and dynamic system signal input processing means 80 may also be included in the computational engine. Analytical computational resources for the computational engine are provided by the computing means 50, which is composed of various registers; accumulators; and timing, synchronization, instruction and address decoding circuitry. The memory means 60 may comprise either or both nonvolatile and volatile memory storage. In the prefered embodiment, non-volatile memory storage is used to contain both low level operating system information and high level application software, and volatile memory is used to store data and other temporary intermediate information. The computing means 50 executes operational code which is stored in both memory types. 
     When required by the type of dynamic system input and output signals, the dynamic system output signal processing means 70 may contain digital data buffers, digital level shifters, digital and/or analog communications interface circuitry, analog to digital conversion, analog signal conditioning and protection circuitry. The dynamic system output signal processing means 70 is used to convert real time dynamic system output signals to a form compatible with the computing means 50. The dynamic system input signal processing means 80 may also contain digital data buffers, digital level shifters, digital and/or analog communications interface circuitry, digital to analog signal conditioning and protection circuitry. This allows the dynanic system input signal processing means 80 to provide either digital or analog inputs to the dynamic system. 
     There are many computational engines which can work in conjunction with the extended horizon adaptive block predictive controller with efficient prediction system. These include 1) microcontrollers and 2) microprocessors with different memory adaptations. New products come out frequently which improve computational speed, improve memory size or speed, add peripheral features, and/or lower cost. Due to the wide variety of computational engines, computing means and memory means which are currently available and rapidly becoming available, each application may have different requirements which will call for a different computational engine configuration. Conversely, old computational engines, and their computing means and memory means, are rapidly becoming obsolete. 
     FIG. 15 shows the preferred embodiment of the Embedded Controller. The Embedded Controller is comprised of a Single Board Computer and Digital Signal Processing (DSP) Function Card, Electrically Programmable Read Only Memory (EPROM) Card, Analog Input Card and Analog Output Card. Analytical computational resources for the Embedded Controller are provided by the Single Board Computer and Digital Signal Processing Function Card which is composed of various registers, accumulators, and timing, synchronization, instruction and address decoding circuitry. The memory portion of the Embedded Controller consists of both non-volatile and volatile memory types. Non-volatile memory storage contained on the EPROM Card is used to contain both low level operating system information and high level application software. While volatile memory contained on the Single Board Computer and Digital Signal Processing Function Card is used to store data and other temporary intermediate information. The Single Board Computer and Digital Signal Processing Function Card executes operational code which is stored in both memory types. The Analog Input Signal Processing component of the Embedded Controller contains digital data buffers, analog to digital conversion, analog signal conditioning and protection circuitry. It is used to convert real time control input signals to a form compatible with the Single Board Computer and Digital Signal Processing Function Card. The Analog Output Signal Processor also has digital data buffers, analog signal conditioning and protection circuitry. However, the Analog Output Signal Processor&#39;s primary function is digital to analog conversion which makes the digital output of the Single Board Computer and Digital Signal Processing Function Card compatible with the real time actuators of the control system. If required, both the Input and Output Signal Processors can condition and function with direct digital information. 
     FIG. 16 is the Single Board Computer and Digital Signal Processing Function Card Block Diagram. The salient feature of the preferred embodiment Embedded Controller is the very highly integrated Single Board Computer and Digital Signal Processing Function Card manufactured by Spectrum Signal Processing Inc. The Single Board Computer and Digital Signal Processing Function Card integrates all of the logic functions from the original IBM AT mother board while including other enhancements such as DSP functions. This highly integrated card facilitates the substantial reduction in controller size, weight, power consumption, and cost. The core of the Single Board Computer and Digital Signal Processing Function Card is a Intel 80386 microprocessor with a full 32 bit instruction set, registers, internal memory organization, all addressing modes, real and protected environment modes, interrupt priorities, and halt and shut down cycles. An external Intel 80387 math coprocessor is used to enhance the floating point execution times. A Texas Instrument TMS320C31 DSP chip is also include, d on the Single Board Computer and Digital Signal Processing Function Card. Said DSP chip features 32 bit instruction and data words, 24 address lines, and 40 bit floating point operations. Said DSp chip also includes on chip peripherals such as a Direct Memory Access (DMA) controller, instruction cache, serial port, two 32 bit timers, two general purpose external flags and four external interrupts. These features provide up to 33 million floating point operations per second (MFLOPS) of DSP performance. Internal memory management functions on the Single Board Computer and Digital Signal Processing Function Card facilitate direct connection to memory to the Intel 80386. An internal DRAM control unit provides direct memory interface of up to 32MB of physical memory. A shadow RAM control feature and a ROM chip select function provides direct interface to a Built In Operating System (BIOS) EPROM. The general AT Core Peripherals such as: an interrupt unit which accepts requests from peripherals and resolves priority levels, a Counter/Timer that provides binary or BCD counting functions, and a Real Time Clock which combines a complete time of day clock with alarm and calendar functions are all included on the Single Board Computer and Digital Signal Processing Function Card. Control inputs and outputs are interfaced to the Single Board Computer and Digital Signal Processing Function Card using a high speed DSP-Link Bus and use the AT ISA Bus only for power. Control inputs are processed using a AT half size Analog Input Card. The card contains analog signal conditioning, analog to digital converters, digital signal input buffers, and analog and digital input signal protection circuitry. It is used to convert real time control input signals to a form compatible with the DSP-Link Bus. The control outputs are processed using a AT half size Analog Output Card. It contains digital to analog converts, analog power amplifiers, digital line drives, special actuator drive circuits, and output protection circuitry. It is used to convert the DSP-Link Bus digital data into analog output signals required for the real time control of the actuators of the control system. 
     Other Embedded Controllers can be based on the computing means which use the Intel 80486, Intel Pentium, Intel i960, R3000, R4000, R6000, Motorola 68HC16, and many others which are familiar to those skilled in the art of computing means. 
     FIG. 17 is a plot of the run time savings the subject predictor system over the Diophantine equation method, as the number of outputs varies. The measured run time for each prediction method (and the difference between them) is shown. The dynamic system has 15 inputs and a model parameter order of 10. The number of outputs is allowed to vary from 1 to 12. Both predictor systems compute predicted outputs for each horizon from k=1 to k=20. The run time savings of the subject prediction system increases cubically with N. As demonstrated in this trial, a dramatic savings in computation time is obtained using the subject predictor system, for large values of N. Thus for dynamic systems having a large number of outputs, significant computational efficiencies are obtainable using the subject predictor system for extended horizon adaptive predictive control. 
     FIG. 18 is a plot of the the run time savings obtained by the use of the subject predictor system, relative to the Diophantine equation method, as the model parameter order n varies. The measured run time for each prediction method (and the difference between them) is shown. Here the dynamic system has 5 inputs and 5 outputs. Both predictor systems perform the task of computing predicted outputs, for each horizon from k=1 to k=30. The model parameter order n varies from n=1 to n=20. The run time savings of the subject prediction system increases quadratically with the model parameter order n. Thus for dynamic systems having a large model parameter order n, significant computational efficiencies are obtainable using the subject predictor system for extended horizon adaptive predictive control. 
     FIG. 19 is a plot of the the run time savings obtained by the subject predictor system, relative to the Diophantine equation method, as number of inputs varies. The measured run time for each prediction method (and the difference between, them) is shown. The dynamic system has 12 outputs and a model parameter order of 10. The number of inputs is varied from 1 to 15. Both predictor systems perform the same task of computing predicted outputs, for each horizon from k=1 to k=20. The run time savings of the subject prediction system increases linearly with M. 
     FIG. 20 is a plot of the run time savings obtained by the subject predictor system, relative to the Diophantine equation method, as the horizon k varies. The dynamic system in this example has 5 inputs and 5 outputs, and a model parameter order of 10. Both predictor systems perform the same task of computing predicted outputs for each horizon, as k is varied from 1 to 30. The run time savings of the subject prediction system increases linearly with prediction horizon k. 
     To further demonstrate the computational efficacy of the predictor system, we embedded the subject system in a controller testbed which uses a 80286 microprocessor and a 80287 math coprocessor. The system considered has two inputs and two outputs with a model order of 2. Here, it is assumed that the system auto-regression dynamic and predictor filter matrices are all diagonal. For this system, a typical execution time on the controller testbed with the subject predictor system is 0.057 seconds. 
     There are numerous applications of extended horizon adaptive block predictive controller with an efficient prediction system in the control of dynamic systems. One system requiring low cost components is the control of the automobiles. The extended horizon adaptive block predictive controller with an efficient prediction system ,can be used to control acceleration, braking and steering of automobiles operating in close proximity to one another at high speed on highways as well as maintenance of the distances between moving vehicles. It can also be used in active suspension systems to smooth out the ride on automobiles or for controlling the suspension systems on magnetically levitated vehicles. It can also be used for engine and transmission control, including control of fuel emissions, engine power, speed, and response of the vehicle to driver .commands. The same applications apply to gasoline, fossil, and electrically powered trucks, transport vehicles, buses and utility vehicles. This includes systems using rail, track and wheels. Other vehicular applications include control of ships, including engine, sail configuration, stabilizer, and hydrofoil control. 
     The speeds and feeds of milling, grinding and cutting machines which process complex parts, particularly those with non homogeneous compositions, can be more effectively controlled with the extended horizon adaptive block predictive controller with an efficient prediction system. Part accuracy can be increased and power consumption can be decreased. This includes fabric cutting. 
     A major application of the subject invention is the control of aircraft, including active control of flight control surfaces, engines and other ancillary systems. 
     Multi-stage steam turbines and jet engine turbines which have inter-stage extractions can have their performance improved with the subject invention. Missiles can be more efficiently targeted and spacecraft pointing and docking accuracy can also be improved with the subject invention. The invention is also particularly useful for robot hand/eye coordination. These are but a few examples of the applications of the subject invention to real time control of dynamic systems. The invention can also be used for simulations of such systems.