Abstract:
A method and system for attenuating the effects of unknown, unmeasurable and time-varying exogenous disturbances on multiple-input multiple-output dynamical systems are described. The disturbance rejection system is characterized in terms of an ARMARKOV or predictive model controller. The parameters of this controller are revised in real time at discrete time steps so as to generate an input to the dynamical system that attenuates the effect of the exogenous disturbance on any chosen set of measured outputs of the dynamical system. The method for revising the controller parameters involves the steps of defining a novel retrospective cost function based on windows of past data, calculating a gradient that is based on this cost function, and using an implementable adaptive step size that brings the controller parameters closer to optimal controller parameters after each revision. The method and system are applicable to active noise and vibration control and reject single-tone, multi-tone, sine sweeping and broadband disturbances in acoustic spaces.

Description:
This invention was made with government support under Grant #F49620-95-1-0019 awarded by the Air Force Office of Scientific Research. The government has certain rights in the invention. 
    
    
     BACKGROUND OF THE INVENTION 
     The field of the invention relates to the problem of rejecting exogenous disturbances acting on dynamical systems (or “plants”). In particular, the invention pertains to a method and system for noise and vibration suppression that does not require measurement of the actual disturbance. 
     Heretofore, several methods required a priori knowledge of the spectral characteristics of the disturbance in addition to models of all four paths in the plant including actuators and sensors shown in FIG. 1, that is, G zw  (the “primary path”), G zu  (the “secondary path”), G yw  (the “reference path”) and G yu  (the “feedback path”). These methods are described in reference  5  of the appended listing of references. Other methods used for active noise control described in references  1 ,  2 ,  6  through  12 , do not require knowledge about the disturbance but often require a direct measurement of it, and require an FIR (finite impulse response) or IIR (infinite impulse response) model of G zu . These methods use instantaneous measurements for adaptation and do not accurately characterize the effect of the control over a window of time. 
     Despite the need for a method and system that can adapt based on retrospective information obtained from sensors to account for the effect of the system and method over a window of time, none was known. Thus, there was the need for a method and system using a retrospective performance evaluation in a special heretofore unknown form. A need also existed to determine an explicit step size or well-defined distance based upon the retrospective performance evaluation. 
     The disclosed method and system of this invention is applicable to a wide class of disturbance rejection problems, including but not limited to active noise and vibration control. Other applications include command-tracking in which the command is viewed as a disturbance signal whose effect is rejected in the output error signal. 
     The present method does not require knowledge of the disturbance spectrum nor a measurement of it, and only requires the numerator of the ARMARKOV model G zu  denoted by the Toeplitz matrix B zu . 
     The present method uses ARMARKOV models to describe the plant including sensors and actuators as well as the disturbance rejection controller. These models are described below. 
     SUMMARY OF THE INVENTION 
     It is an object of the invention to provide a method and system that evaluates performance based upon past data and determines an explicit step size or distance for adaptation, for differentiation from existing methods. 
     In contrast to the prior art, the method and system for achieving such rejection comprises of a set of sensors which measure the outputs of the plant for which the effect of the disturbance is to be minimized, an optional additional set of sensors which measure other outputs of the plant, converters that digitize analog signals from the sensors, a microprocessor capable of storing data from the converters and performing the calculations described in the method herein, converters that create command signals from the results of the calculations of the microprocessor, and actuators that act on the plant based on these command signals. A graphical representation of the system according to the invention is given in FIG.  1 . 
     The plant with sensors and actuators comprise the four block unit in FIG. 1, while the microprocessor implementing the method described herein, or “controller”, is the lower block marked G c . The exogenous disturbance is denoted by w(k), the outputs of the system for which the effects of this disturbance are to be minimized or “performance measurements” by z(k), the additional measured outputs of the system or “feedback measurements” by y(k) and the signals generated by the controller or “control signals” by u(k). The paths from the two sets of inputs to the plant, namely, w(k) and u(k), to the two sets of outputs of the plant, namely, z(k) and y(k) are denoted by G zw , G zu , G yw  and G yu . 
     For a more complete understanding of the present invention, reference is made to the following detailed description when read in conjunction with the accompanying drawings wherein like reference characters refer to like elements throughout the several views, in which: 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 illustrates a graphical representation of the system according to the invention; 
     FIG. 2 illustrates a geometrical interpretation of the method according to the invention; 
     FIG. 3 illustrates a graphical representation of an experimental set-up for the invention; 
     FIG. 4 illustrates the results of the system in active mode compared with the results of the system in inactive mode for a single-tone disturbance at 139.65 Hz; 
     FIG. 5 illustrates the results of the system in active mode compared with the results of the system in inactive mode for a dual-tone disturbance at 139.74 Hz and 160.4 Hz; 
     FIG. 6 illustrates the results of the system in active mode compared with the results of the system in inactive mode for band-limited white noise; and 
     FIG. 7 illustrates the results of the system in active mode compared with the results of the system in inactive mode with AM radio disturbance. 
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENT 
     To begin, we describe the ARMARKOV model of the nth-order discrete-time finite-dimensional linear time-invariant system given by 
     
       
           x ( k+ 1)= Ax ( k )+ Bu ( k ),  (1) 
       
     
     k=0, 1, 2, . . . , 
     
       
           y ( k )= Cx ( k )+ Du ( k ),  (2) 
       
     
     where A, B, C and D are real matrices of appropriate size, u(k) is of size m u  and y(k) is of size l y , and whose Markov parameters H j  of size l y ×m u  are defined as                  H   j          =   Δ        D     ,       for                 j     =     -   1       ,           (   3   )                                =   Δ            CA   j        B       ,                  for                 j     ≥   0.               (   4   )                                
     This system (1), (2) may be alternatively described by the auto-regressive moving average (ARMA) representation given by 
     
       
           y ( k )=− a   1   y ( k− 1)− . . . − a   n   y ( k−n )+ B   0   u ( k )+ . . . + B   n   u ( k−n ),  (5) 
       
     
     or the μ-ARMARKOV (ARMA+Markov) model or μ step ahead predictor model                      y        (   k   )       =                    ∑     j   =   1     n            -     α   j            y        (     k   -   μ   -   j   +   1     )           +                                  ∑     j   =   1     μ            H     j   -   1            u        (     k   -   j   +   1     )           +                                  ∑     j   =   1     n            B   j          u        (     k   -   μ   -   j   +   1     )           ,                   (   6   )                                
     where α j  are scalars and B j  are of size l y ×m u , j=1, . . . , n. We note that in the special case μ=1, the ARMARKOV form (6) is the same as the ARMA form. 
     Now, let p denote the data window length and define the extended measurement vector Y(k) of size l y p and the ARMARKOV regressor vector Φ yu (k) of size l y (p+n−1)+m u (μ+p+n−1) by                  Y        (   k   )            =   Δ          [           y        (   k   )               ⋮             y        (     k   -   p   +   1     )             ]       ,                    Φ   yu          (   k   )            =   Δ            [           y        (     k   -   μ     )               ⋮             y        (     k   -   μ   -   p   -   n   +   2     )                 u        (   k   )               ⋮             u        (     k   -   μ   -   p   -   n   +   2     )             ]     .               (   7   )                                
     Using (6), Y(k) and Φ yu (k) are related by 
     
       
           Y ( k )= W   yu Φ yu ( k ),  (8) 
       
     
     where the block-Toeplitz ARMARKOV weight matrix W yu  of size pl y ×[l y (p+n−1)+m u (μ+p+n−1)] is defined by                W   yu            =   Δ            [             -     α   1            I     I   y             …           -     α   n            I     t   y               0     t   y           ⋯         0     I   y             H     -   1           ⋯         H     μ   -   2             B   1         ⋯         B   n           0       t   y     ×     m   u             ⋯         0       t   y     ×     m   u                   0     t   y           ⋰                   ⋰       ⋰       ⋮         0       t   y     ×     m   u             ⋰                   ⋰       ⋰                   ⋰       ⋰       ⋮           ⋮       ⋰       ⋰                   ⋰         0     t   y           ⋮       ⋰       ⋰                   ⋰       ⋰                   ⋰         0       t   y     ×     m   u                   0     t   y           ⋯         0     t   y               -     α   1            I     t   y             ⋯           -     α   n            I     t   y               0       t   y     ×     m   u             ⋯         0       t   y     ×     m   u               H     -   1           ⋯         H     μ   -   2             B   1         …         B   n           ]                .               (   9   )                                
     We now develop the ARMARKOV/Toeplitz model of the two vector input, two vector output plant with sensors and actuators whose inputs are the disturbance w(k) and the control u(k), and whose outputs are the feedback measurement y(k) and the performance measurement z(k) as shown in FIG.  1 . The ARMARKOV form of the plant with actuators and sensors is                      z        (   k   )       =                    ∑     j   =   1     n            -     α   j            z        (     k   -   μ   -   j   +   1     )           +       ∑     j   =   1     μ            H     zw   ,     j   -   2              w        (     k   -   j   +   1     )           +                                  ∑     j   =   1     n            B     zw   ,   j            w        (     k   -   μ   -   j   +   1     )           +       ∑     j   =   1     μ            H     zu   ,     j   -   2              u        (     k   -   j   +   1     )           +                                  ∑     j   =   1     n            B     zu   ,   j            u        (     k   -   μ   -   j   +   1     )           ,                   (   10   )                       y        (   k   )       =                    ∑     j   =   1     n            -     α   j            y        (     k   -   μ   -   j   +   1     )           +       ∑     j   =   1     μ            H     yw   ,     j   -   2              w        (     k   -   j   +   1     )           +                                  ∑     j   =   1     n            B     yw   ,   j            w        (     k   -   μ   -   j   +   1     )           +       ∑     j   =   1     μ            H     yu   ,     j   -   2              u        (     k   -   j   +   1     )           +                                  ∑     j   =   1     n            B     yu   ,   j            u        (     k   -   μ   -   j   +   1     )           ,                   (   11   )                                
     where α j  are scalars, B zw,j  and H zw,j  are of size l z ×m w , B zu,j  and H zu,j  are of size l z ×m u , B yw,j  and H yw,j  are of size l y ×m w , and B yu,j  and H yu,j  are of size l y ×m u . 
     Next, define the extended performance measurement vector Z(k), the extended feedback measurement vector Y(k) and the extended control vector U(k) by                  Z        (   k   )            =   Δ          [           z        (   k   )               ⋮             z        (     k   -   p   +   1     )             ]       ,       Y        (   k   )            =   Δ          [           y        (   k   )               ⋮             y        (     k   -   p   +   1     )             ]       ,                  U        (   k   )            =   Δ          [           u        (   k   )               ⋮             u        (     k   -   p   +   1     )             ]       ,           (   12   )                                
     where the controller window size p c  is given by μ+n+p−1, and the ARMARKOV regressor vectors Φ zw (k) and Φ yw (k) are defined by                    Φ   zw          (   k   )            =   Δ          [           z        (     k   -   μ     )               ⋮             z        (     k   -   μ   -   p   -   n   +   2     )                 w        (   k   )               ⋮             w        (     k   -   μ   -   p   -   n   +   2     )             ]       ,                    Φ   yw          (   k   )            =   Δ            [           y        (     k   -   μ     )               ⋮             y        (     k   -   μ   -   p   -   n   +   2     )                 w        (   k   )               ⋮             w        (     k   -   μ   -   p   -   n   +   2     )             ]     .               (   13   )                                
     Furthermore, define the block-Toeplitz ARMARKOV weight matrices W zw  of size pl z ×[(n+p−1)l z +(μ+n+p−1)m w ] and W yw  of size pl y ×[(n+p−1)l y +(μ+n+p−1)m w ] by                  B   zu          =   Δ          [             -     α   1            I     l   z             …           -     α   n            I     l   z               0     l   z           ⋯         0     l   z             H     zw   ,     -   1             ⋯         H     zw   ,     μ   -   2               B     zw   ,   1           ⋯         B     zw   ,   n             0       l   z     ×     m   w             ⋯         0       l   z     ×     m   w                   0     l   z           ⋰                   ⋰       ⋰       ⋮         0       l   z     ×     m   w             ⋰                   ⋰       ⋰                   ⋰       ⋰       ⋮           ⋮       ⋰       ⋰                   ⋰         0     l   z           ⋮       ⋰       ⋰                   ⋰       ⋰                   ⋰         0       l   z     ×     m   w                   0     l   z           ⋯         0     l   z               -     α   1            I     l   z             ⋯           -     α   n            I     l   z               0       l   z     ×     m   w             ⋯         0       l   z     ×     m   w               H     zw   ,     -   1             ⋯         H     zw   ,     μ   -   2               B     zw   ,   1           …         B     zw   ,   n             ]       ,           (   14   )                   W   yw          =   Δ          [             -     α   1            I     l   y             …           -     α   n            I     l   y               0     l   y           ⋯         0     l   y             H     yw   ,     -   1             ⋯         H     yw   ,     μ   -   2               B     yw   ,   1           ⋯         B     yw   ,   n             0       l   y     ×     m   w             ⋯         0       l   y     ×     m   w                   0     l   y           ⋰                   ⋰       ⋰       ⋮         0       l   y     ×     m   w             ⋰                   ⋰       ⋰                   ⋰       ⋰       ⋮           ⋮       ⋰       ⋰                   ⋰         0     l   y           ⋮       ⋰       ⋰                   ⋰       ⋰                   ⋰         0       l   y     ×     m   w                   0     l   y           ⋯         0     l   y               -     α   1            I     l   y             ⋯           -     α   n            I     l   y               0       l   y     ×     m   w             ⋯         0       l   z     ×     m   w               H     yw   ,     -   1             ⋯         H     yw   ,     μ   -   2               B     yw   ,   1           …         B     yw   ,   n             ]       ,           (   15   )                                
     and the block-Toeplitz ARMARKOV control matrices B zu  of size pl z ×p c m u  and B yu  of size pl y ×p c l u  by                  B   zu          =   Δ          [           H     zu   ,     -   1             ⋯         H     zu   ,     μ   -   2               B     zu   ,   1           ⋯         B     zu   ,   1             0       l   z     ×     m   u             ⋯         0       l   z     ×     m   u                   0       l   z     ×     m   u             ⋰                   ⋰       ⋰                   ⋰       ⋰       ⋮           ⋮       ⋰       ⋰                   ⋰       ⋰                   ⋰         0       l   z     ×     m   u                   0       l   z     ×     m   u             ⋯         0       l   z     ×     m   u               H     zu   ,     -   1             ⋯         H     zu   ,     μ   -   2               B     zu   ,   1           ⋯         B     zu   ,   n             ]       ,           (   16   )                 B   yu            =   Δ            [           H     yu   ,     -   1             ⋯         H     yu   ,     μ   -   2               B     yu   ,   1           ⋯         B     yu   ,   1             0       l   y     ×     m   u             ⋯         0       l   y     ×     m   u                   0       l   y     ×     m   u             ⋰                   ⋰       ⋰                   ⋰       ⋰       ⋮           ⋮       ⋰       ⋰                   ⋰       ⋰                   ⋰         0       l   z     ×     m   u                   0       l   y     ×     m   u             ⋯         0       l   y     ×     m   u               H     yu   ,     -   1             ⋯         H     yu   ,     μ   -   2               B     yu   ,   1           ⋯         B     yu   ,   n             ]     .               (   17   )                                
     Then (10) and (11) can be written in the form 
     
       
           Z ( k )= W   zw Φ zw ( k )+ B   zu   U ( k ),  (18) 
       
     
     
       
           Y ( k )= W   yw Φ yw ( k )+ B   yu   U ( k ),  (19) 
       
     
     which is the ARMARKOV/Toeplitz model of the augmented plant. 
     Next, we formulate an adaptive disturbance rejection feedback algorithm for the system represented by (18) and (19). We use a strictly proper controller G c  in ARMARKOV form of order n c  with μ c  Markov parameters, so that, analogous to (6), the control input u(k) is given by                      u        (   k   )       =                    ∑     j   =   1       n   c              -       α     c   ,   j            (   k   )              u        (     k   -     μ   c     -   j   +   1     )           +                                  ∑     j   =   1         μ   c     -   1                H     c   ,     j   -   1              (   k   )            y        (     k   -   j   +   1     )           +                                  ∑     j   =   1       n   c                B     c   ,   j            (   k   )            y        (     k   -     μ   c     -   j   +   1     )           ,                   (   20   )                                
     where the controller Markov parameter H c,j  is of size m u ×l y . Next, define the controller parameter block vector θ(k) by                θ        (   k   )            =   Δ          [             -       α     c   ,   1            (   k   )              I     m   u             ⋯           -       α     c   ,     n   c              (   k   )              I     m   u                 H     c   ,   0            (   k   )           ⋯           H     c   ,       μ   c     -   2              (   k   )               B     c   ,   1            (   k   )           ⋯               B     c   ,     n   c              (   k   )       ]     ,                     (   21   )                                
     where θ(k) is of size m u ×[n c m u +(n c +μ c −1)l y ]. Now from (12) and (20) it follows that u(k) and U(k) are given by 
     
       
           u ( k )=θ( k ) R   1 Φ uy ( k )  (22) 
       
     
     and                  U        (   k   )       =       ∑     i   =   1       p   c              L   i          θ        (     k   -   i   +   1     )            R   i            Φ   uy          (   k   )             ,           (   23   )                                
     where                    Φ   uy          (   k   )            =   Δ          [           u        (     k   -     μ   c       )               ⋮             u        (     k   -     μ   c     -     n   c     -     p   c     +   2     )                 y        (     k   -   1     )               ⋮             y        (     k   -     μ   c     -     n   c     -     p   c     +   2     )             ]       ,           (   24   )                                
     and where                L   i          =   Δ          [           0       (     i   -   1     )          m   u     ×     m   u                   I     m   u                 0       (       p   c     -   i     )          m   u     ×     m   u               ]             (   25   )                                
     is of size p c m u ×m u  with I m     u    denoting the identity matrix of size m u , and                R   i          =   Δ          [           0     q1   ×     (     i   -   1     )          m   u               I     q1   ×   q1             0     q1   ×     (       p   c     -   i     )          m   u               0     q1   ×     (     i   -   1     )          l   y               0     q1   ×   q2             0     q1   ×     (       p   c     -   i     )          l   y                   0     q2   ×     (     i   -   1     )          m   u               0     q2   ×   q1             0     q2   ×     (       p   c     -   i     )          m   u               0     q1   ×     (     i   -   1     )          l   y               I     q2   ×   q2             0     q2   ×     (       p   c     -   i     )          l   y               ]             (   26   )                                
     is of size [n c m u +(n c +μ c −1)l y ]×[(n c +p c −1)m u +(n c +μ c +p c −2)l y ], with q 1 n c m u  and q 2 (n c +μ c −1)l y . Thus, from (18) and (23) we obtain                Z        (   k   )       =         W   zw            Φ   zw          (   k   )         +       B   zu            ∑     i   =   1       p   c              L   i          θ        (     k   -   i   +   1     )            R   i              Φ   uy          (   k   )       .                     (   27   )                                
     Next, we describe the update law for the controller parameter block vector θ(k). To do this, we define a retrospective performance cost function that evaluates the performance of the controller obtained from the current value of θ(k) based upon the measurements of the system during the previous p c  steps. Therefore, we define the estimated performance {circumflex over (Z)}(k) by                    Z   ^          (   k   )            =   Δ              W   zw            Φ   zw          (   k   )         +       B   zu            ∑     i   =   1       p   c              L   i          θ        (   k   )            R   i            Φ   uy          (   k   )                 ,           (   28   )                                
     which has the same form as (27) but with θ(k−i+1) replaced by the current controller parameter block vector θ(k). Using (28) we define the retrospective performance cost function 
     
       
           J ( k )=½ {circumflex over (Z)}   T ( k ) {circumflex over (Z)} ( k ),  (29) 
       
     
     with “T” denoting the transpose of a vector. Next, the gradient of J(k) with respect to θ(k) is given by                  ∂     J        (   k   )           ∂     θ        (   k   )           =       ∑     i   =   1       p   c              L   i   T          B   zu   T            Z   ^          (   k   )              Φ   uy   T          (   k   )              R   i   T     .                 (   30   )                                
     Since w(k) is not available, which implies that Φ zw (k) is unknown, {circumflex over (Z)}(k) cannot be calculated from (28). However, it follows from (18) and (28) that                    Z   ^          (   k   )       =       Z        (   k   )       -       B   zu          (       U        (   k   )       -       ∑     i   =   1       p   c              L   i          θ        (   k   )            R   i            Φ   uy          (   k   )             )           ,           (   31   )                                
     which can be used to evaluate (30). 
     The gradient (30) is used in the update law                  θ        (     k   +   1     )       =       θ        (   k   )       -       η        (   k   )              ∂     J        (   k   )           ∂     θ        (   k   )                 ,           (   32   )                                
     where η(k) is the adaptive step size. To determine the adaptive step size η(k), we assume that there is a controller parameter block vector θ* that minimizes J(k) for all k. The method does not need to know θ*. 
     Now, we define the desired performance                      Z   ^     *          (   k   )            =   Δ              W   zw            Φ   zw          (   k   )         +       B   zu            ∑     i   =   1       p   c              L   i          θ   *          R   i            Φ   uy          (   k   )                 ,           (   33   )                                
     and the performance error                ɛ        (   k   )            =   Δ                Z   ^     *          (   k   )       -         Z   ^          (   k   )       .               (   34   )                                
     Our goal is to determine η(k) such that θ(k) moves closer to θ* after each update. For convenience, we define the optimal adaptive step size                    η   opt          (   k   )            =   Δ                   ɛ        (   k   )            2   2                ∂     J        (   k   )           ∂     θ        (   k   )                F   2         ,           (   35   )                                
     where ∥ ∥ F  denotes the matrix Frobenius norm and ∥ ∥ 2  denotes the vector Euclidean norm. It is shown in reference 16 that 
     
       
         ∥θ*−θ( k+ 1)∥ F ≦∥θ*−θ( k )∥ F   (36) 
       
     
     if and only if η(k) satisfies 
     
       
         0&lt;η( k )&lt;2η opt ( k ).  (37) 
       
     
     Furthermore, η(k)=η opt (k) minimizes ∥θ*−θ(k+1)∥ F −∥θ*−θ(k)∥ F . 
     A geometrical interpretation of the procedure detailed above is now presented. Using FIG. 2 for reference, the objective of the algorithm is to move the controller parameter block vector θ(k) closer to the optimal controller parameter block vector θ*. The direction in which to move is the negative of the gradient          ∂     J        (   k   )           ∂     θ        (   k   )                                
     which is obtained from the retrospective performance cost function. The distance to move at each time step is determined by the adaptive step size η(k). It is shown that the step size η opt (k) moves θ(k) to the point closest to θ* along the negative gradient direction, that is, to a point such that the vectors E(k+1) and        -       ∂     J        (   k   )           ∂     θ        (   k   )                                  
     are at right angles. 
     In practice, η opt (k) is not computable since ε(k) is not available from sensor measurements. The crucial innovative feature of the method in this invention is the use of an implementable adaptive step size which can be calculated from available data and is guaranteed to be within the range (37) that mathematically demonstrates that θ(k) moves closer to θ*. Three such step sizes are given below:                    η   1          (   k   )            =   Δ          1       (     [       ∑     i   =   1       p   c                       R   i            Φ   uy          (   k   )              2            σ   _          (       B   zu          L   i       )           ]     )     2         ,           (   38   )                     η   2          (   k   )            =   Δ          1                  Φ   uy          (   k   )            2   2          [       ∑     i   =   1       p   c              σ   _          (       B   zu          L   i       )         ]       2         ,           (   39   )                   n   3          (   k   )            =   Δ          1       p   c                 B   zu          F   2                   Φ   uy          (   k   )            2   2                 (   40   )                                
     where {overscore (σ)}(B zu L i ) denotes the maximum singular value of the matrix B zu L i . Note that if B zu  is known, then η 1 (k), η 2 (k) or η 3 (k) can be calculated and used to implement (32). Other implementable adaptive step sizes satisfying (37) may be obtained. The step sizes η 1 (k), η 2 (k) and η 3 (k) satisfy 
     
       
         0&lt;η 3 ≦η 2 ( k )≦η 1 ( k )≦η opt ( k ),  (41) 
       
     
     and thus satisfy (37). 
     The steps involved in implementing the adaptive algorithm are as follows: 
     0. Obtaining the matrix B zu  using the identification algorithm of reference 14, 15 or by calculating from an ARMA or state space representation of G zu . 
     1. Calculating the control signal u(k) from the controller parameter block vector θ(k) and the vector Φ uy (k) using (20). 
     2. Using the signals u(k), z(k) and y(k) updating the estimated performance vector {circumflex over (Z)}(k) as defined in (31). 
     3. Calculating the retrospective gradient          ∂     J        (   k   )           ∂     θ        (   k   )                                
     using (30). 
     4. Calculating an implementable adaptive step size such as η 1 (k), η 2 (k) or η 3 (k) from (38), (39) or (40). 
     5. Revising the controller parameter block vector θ(k) using (32). 
     6. Updating Φ uy (k+1) as defined in (24). 
     Steps 1 through 5 are performed at each time step k. 
     Experimental demonstration of the ARMARKOV adaptive disturbance algorithm for active noise control is performed on an acoustic duct of circular cross section. The duct is 80 inches long and has a diameter of 4 inches. The disturbance speaker (w) is located at one end of the duct and the measurement sensor (y), a microphone, is located 4 inches in from the same end of the duct. The performance sensor (z), a microphone, is positioned 6 inches in from the other end. Alternative sensors for vibration control are accelerometers and piezo-electric sensors. The control actuator (u), a speaker, is placed 16 inches in from that end of the duct. A servovalve for flow modulation of compressed air is another form of actuation for noise control while proof mass actuators can be used for vibration control. The signals from the two microphones are amplified by a dbx 760x microphone preamplifier while the control signal is amplified by an Alesis RA-100 amplifier. Both speakers are Radio Shack 6 inch woofers. A graphical representation of the experimental set-up is shown in FIG.  3 . 
     The algorithm is tested on four types of disturbances, namely, a single-tone disturbance (139.65 Hz), a two-tone disturbance (135.74 Hz and 160.4 Hz), band-limited white noise (up to 390 Hz) and AM radio noise. The algorithm uses n=4 and μ=12 for the secondary path matrix B zu , and n c =2, μ c =10 and p=2 for control. The controller is implemented on a dSPACE ds1102 real time board running a TMS320C30 DSP processor at a sampling frequency of 800 Hz. The microphone signals are processed through an Ithaco DL 4302 low pass filter that rolls off at 315 Hz. The tonal and band-limited white noise disturbances are generated by a Stanford Research Systems 770 FFT network analyzer and amplified by an Optimus STA-825 stereo receiver. 
     FIG. 4 shows the acoustic response with the disturbance rejection system inactive (“open-loop”) and with the disturbance rejection system active (“closed-loop”) with a single-tone disturbance. Disturbance attenuation of more than 40 dB is achieved with convergence in about 1 second. The system and method provide the same level of attenuation by adaptation when the frequency of the disturbance tone is changed, as in sine sweeps, while the system is active. 
     FIG. 5 shows the open-loop and closed-loop performance with a two-tone disturbance. In this case, disturbance attenuation of more than 35 dB is observed. FIG. 6 shows the open-loop and closed-loop magnitude plots of the transfer function from disturbance to performance with a white noise disturbance, and noise suppression of up to 15 dB is observed over a frequency range from 0 to 300 Hz. Finally, FIG. 7 shows the open-loop and closed-loop frequency response with an AM radio disturbance. Noise reduction levels of up to 40 dB are observed over the frequency range 0 to 300 Hz. 
     In contrast and improvement to the prior art, the present method has three innovative features. The first is the use of ARMARKOV/Toeplitz structures for describing both the plant and controller. While these structures have been used for predictive and neural net control as described in references 3, 4 and 13, the present method uses them in retrospective fashion to obtain a controller update law that learns from past data. 
     The second innovation is the definition of the retrospective cost function and calculation of the gradient with respect to this cost function. In the prior art, instantaneous or predicted cost functions are used. 
     The third innovation is the use of an implementable adaptive step size for the controller update which guarantees that the controller parameters move closer to the unknown optimal controller parameters at each time step. 
     Having described the invention, many modifications thereto will become apparent to those skilled in the art to which it pertains without deviation from the spirit of the invention as defined in the appended claims. 
     REFERENCES 
     [ 1 ] G. E. Warnaka, Active Noise Control, U.S. Pat. No. 5,511,127, April 1996. 
     [ 2 ] M. K. Ferguson, S. C. Southward and M. C. Heath, Active Noise and Vibration Cancellation System, U.S. Pat. No. 5,619,581, April 1997. 
     [ 3 ] D. W. Cawlfield, Model Predictive Control Apparatus and Method, U.S. Pat. No. 5,519,605, May 1996. 
     [ 4 ] D. C. Hyland and J. N. Juang, Series Parallel Approach to Identification of Dynamic Systems, U.S. Pat. No. 5,680,513, October 1997. 
     [ 5 ] B. Francis,  A Course in H∞ Control Theory , New York: Springer-Verlag, 1987. 
     [ 6 ] P. A. Nelson and S. J. Elliot,  Active Control of Sound.  New York: Academic Press, 1992. 
     [ 7 ] S. M. Kuo and D. R. Morgan,  Active Noise Control Systems.  New York: Wiley, 1996. 
     [ 8 ] L. A. Sievers and A. H. von Flotow, “Comparison and Extensions of Control Methods for Narrow Band Disturbance Rejection,”  IEEE Transactions on Signal Processing,  Vol. 40, pp. 2377-2391, 1992. 
     [ 9 ] W. Messner and M. Bodson, “Design of Adaptive Feedforward Alogrithms Using Internal Model Equivalence,”  International Journal of Adaptive Control and Signal Processing,  Vol. 9, pp. 199-212, 1995. 
     [ 10 ] S. J. Elliot, I. M. Stothers and P. A. Nelson. “A Multiple Error LMS Algorithm and its Applications to the Active Control of Sound and Vibration,”  IEEE Transactions on Acoustics, Speech and Signal Processing,  Vol. ASSP-35, pp. 1423-1434, 1987. 
     [ 11 ] B. Widrow and E. Walach,  Adaptive Inverse Control.  NJ: Prentice Hall, 1996. 
     [ 12 ] F. Jiang, H. Tsuji, N. Ojiro, H. Ohmori and A. Sano, “Adaptation for Active Noise Control,”  IEEE Control Systems Magazine,  pp. 36-47, December 1997. 
     [ 13 ] K. J. Astrom and B. Wittenmark,  Adaptive Control,  second edition, Reading Mass.: Addison-Wesley, 1995. 
     [ 14 ] J. C. Akers and D. S. Bernstein, “Time-Domain Identification Using ARMARKOV/Toeplitz Models,”  Proc. Amer. Contr. Conf.,  pp. 191-195, Albuquerque, N.Mex., June 1997. 
     [ 15 ] J. C. Akers and D. S. Bernstein, “ARMARKOV Least-Squares Identification,”  Proc. Amer. Contr. Conf.,  pp. 186-190, Albuquerque, N.Mex., June 1997. 
     [ 16 ] R. Venugopal and D. S. Bernstein, “Adaptive Disturbance Rejection Using ARMARKOV System Representations,”  Proc. IEEE Conf. Dec. Contr.,  pp. 1884-1889, San Diego, Calif., December 1997.