3.2.2. Recursive Least Squares (RLS) Algorithm

Given an input samples set {*u*(1), *<sup>u</sup>*(2),..., *u*(*N*)} and a desired response set {*d*(1), *d*(2),..., *d*(*N*)}, the output *y*(*n*) could be computed by linear filter method as follows,

$$y(n) = \sum\_{k=0}^{M} w\_k u(n-k), n = 0, 1, 2, \dots \tag{56}$$

Minimizing the sum of the error squares, we have,

$$\begin{aligned} \varepsilon(n) &= \varepsilon(\omega\_1(n), \omega\_2(n), \dots, \omega\_{M-1}(n)) \\ &= \sum\_{i=i\_1}^n \beta(n, i) \times \left[ \varepsilon(i)^2 \right] \\ &= \sum\_{i=i\_1}^n \beta(n, i) \times \left[ d(i) - \sum\_{k=0}^{M-1} \omega\_k(n) \mu(i - k) \right]^2 \end{aligned} \tag{57}$$

where the error *e*(*i*) is *e*(*i*) = *d*(*i*) − *y*(*i*) = *d*(*i*) − *M*−1 ∑ *k*=0 *<sup>ω</sup>k*(*n*)*u*(*<sup>i</sup>* − *k*), and the *β*(*<sup>n</sup>*, *i*) is the forgetting factor and 0 < *β*(*<sup>n</sup>*, *i*) ≤ 1, *i* = 1, 2, ... , *n*. To simplify the writing format, taking the form as *β*(*<sup>n</sup>*, *i*) = *<sup>λ</sup><sup>n</sup>*−*i*, *i* = 1, 2, . . . , *n*. Thus, the sum of error squares can be rewritten as,

$$\begin{aligned} \varepsilon(n) &= \varepsilon(\omega\_1(n), \omega\_2(n), \dots, \omega\_{M-1}(n)) \\ &= \sum\_{i=i\_1}^n \lambda^{n-i} \times \left[ d(i) - \sum\_{k=0}^{M-1} \omega\_k(n) u(i-k) \right]^2 \end{aligned} \tag{58}$$

Defining the following formulas,

$$u'(i) = \sqrt{\lambda^{n-i}} u(i); \quad d'(i) = \sqrt{\lambda^{n-i}} d(i) \tag{59}$$

Then the Equation (58) can be rewritten as,

$$\begin{split} \varepsilon(n) &= \sum\_{i=i\_1}^{n} \lambda^{n-i} \times \left[ d(i) - \sum\_{k=0}^{M-1} \omega\_k(n) u(i-k) \right]^2 \\ &= \sum\_{i=i\_1}^{n} \left[ d'(i) - \sum\_{k=0}^{M-1} \omega\_k(n) u'(i-k) \right]^2 \end{split} \tag{60}$$

which is the standard least squares (LS) criterion, the solution of the LS can be obtained as,

$$\begin{split} \omega(n) &= \left(\sum\_{i=i\_1}^{n} u'(i)u'(i)^T\right)^{-1} \sum\_{i=i\_1}^{n} u'(i)d'(i) \\ &= \left(\sum\_{i=i\_1}^{n} \lambda^{n-i} u(i)u(i)^T\right)^{-1} \sum\_{i=i\_1}^{n} \lambda^{n-i} u(i)d(i) \\ &= \left[\Phi(n)\right]^{-1} \Psi(n) \end{split} \tag{61}$$

where parameters <sup>Φ</sup>(*n*) = *n* ∑ *i*=1 *<sup>λ</sup><sup>n</sup>*−*iu*(*i*)*u*(*i*)*<sup>T</sup>* and *ψ*(*n*) = *n*∑ *i*=1 *<sup>λ</sup><sup>n</sup>*−*iu*(*i*)*d*(*i*). Based on *ω*(*n*) = [Φ(*n*)]−<sup>1</sup>*ψ*(*n*), thus the *ω*(*n*) at time *n* − 1 could be obtained, i.e.,

$$
\omega(n-1) = \left[\Phi(n-1)\right]^{-1} \psi(n-1) \tag{62}
$$

Therefore, the variables <sup>Φ</sup>(*n*) and *ψ*(*n*) can be rewritten using <sup>Φ</sup>(*n* − 1) and *ψ*(*n* − <sup>1</sup>), i.e.,

$$\begin{aligned} \Phi(n) &= \sum\_{i=1}^{n} \lambda^{n-i} u(i) u(i)^T = \lambda \sum\_{i=1}^{n-1} \lambda^{n-1-i} u(i) u(i)^T + u(n) u(n)^T \\ &= \lambda \Phi(n-1) + u(n) u(n)^T \end{aligned} \tag{63}$$

and

$$\begin{aligned} \psi(n) &= \sum\_{i=1}^{n} \lambda^{n-i} u(i) d(i) = \lambda \sum\_{i=1}^{n-1} \lambda^{n-1-i} u(i) d(i) + u(n) d(n) \\ &= \lambda \psi(n-1) + u(n) d(n) \end{aligned} \tag{64}$$

The matrix inversion formula of <sup>Φ</sup>−<sup>1</sup>(*n*) is as follows,

$$\Phi^{-1}(n) = \lambda^{-1}\Phi^{-1}(n-1) - \frac{\lambda^{-2}\Phi^{-1}(n-1)u(n)u^T(n)\Phi^T(n-1)}{1 + \lambda^{-1}u^T(n)\Phi^{-1}(n-1)u(n)}\tag{65}$$

Denoting

$$P(n) = \Phi^{-1}(n)\tag{66}$$

and

$$k(n) = \frac{\lambda^{-1} P(n-1) u(n)}{1 + \lambda^{-1} u^T(n) P(n-1) u(n)} = P(n) u(n) \tag{67}$$

Therefore,

$$P(n) = \lambda^{-1} P(n-1) - \lambda^{-1} k(n) u^T(n) P(n-1) \tag{68}$$

The main time-update equation *ω*(*n*) can be derived as,

$$\begin{aligned} \omega(n) &= [\Phi(n)]^{-1} \psi(n) = P(n)\psi(n) \\ &= P(n)(\lambda \psi(n-1) + u(n)d(n)) \\ &= P(n)(\lambda \Phi(n-1)\omega(n-1) + u(n)d(n)) \\ &= P(n)\left\{ [\Phi(n) - u(n)u^T(n)]\omega(n-1) + u(n)d(n) \right\} \\ &= \omega(n-1) - P(n)u(n)u^T(n)\omega(n-1) + P(n)u(n)d(n) \\ &= \omega(n-1) + P(n)u(n)\left[d(n) - u^T(n)\omega(n-1)\right] \\ &= \omega(n-1) + P(n)u(n)a(n) \\ &= \omega(n-1) + k(n)a(n) \end{aligned} \tag{69}$$

where *α*(*n*) = *d*(*n*) − *u<sup>T</sup>*(*n*)*ω*(*n* − 1) is the innovation process. In conclusion, the prediction processes of the ARMA combined with RLS algorithm are summarized as follows,

*Step 1: Algorithm initialization*. The forgetting factor 0 < *λ* ≤ 1, the predicted steps *n*.

*Step 2: Initial parameter calculation*. Calculate the initial order *p* and *q* of ARMA (*p*, *q*) model using Akaike information criterion (AIC) [50] based on historical data.

*Step 3: Coefficients extraction.* Extract the coefficients *φ*1, *φ*2, ... , *φp* and *θ*1, *θ*2, ... , *θq* as the input *u*(*i*) of RLS algorithm, i.e., *u*(*i*) = 1, *φ*1, *φ*2,..., *φ<sup>p</sup>*, *θ*1, *θ*2,..., *<sup>θ</sup>q*, *i* = 1, 2, . . . , *p* + *q*.

Step 4: Prediction iteration based on RLS algorithm. Repeat the following iterations,

$$k(i) = \frac{\lambda^{-1}P(i-1)u(i)}{1 + \lambda^{-1}u^T(i)P(i-1)u(i)};$$

$$a(i) = d(i) - u^T(i)\omega(i-1);$$

$$\omega(i) = \omega(i-1) + k(i)a(i);$$

$$P(i) = \lambda^{-1}P(i-1) - \lambda^{-1}k(i)u^T(i)P(i-1);$$

$$y(i) = y(i-1) + k \times \alpha(i);$$

*Step 5: Judgment of the end conditions*. If the iterative step is satisfied, then output predicted signal *y*(*n*), otherwise, *i* = *i* + 1, and go to step (4).
