**2.2 Variable, fractional-order linear time-invariant difference equations**

On the base of the Grünwald-Letnikov variable, fractional-order linear timeinvariant backward-difference the difference Eqs. (GL-VFOBE) for *i* ¼ 1, 2, ⋯, *p* and *j* ¼ 1, 2, ⋯, *q* representing discrete models of real dynamical systems or discrete control strategies are defined by the variable, fractional-order linear time-invariant difference equation (VFODE). *h*> 0 denotes the sampling time.

$$\sum\_{l=0}^{n\_i} a\_{i,l\_{k\_0}} \Delta\_k^{\left[\iota\_{il}(k)\right]} y(kh) = \sum\_{l=0}^{m\_i} b\_{i,l\_{k\_0}} \Delta\_k^{\left[\mu\_{il}(k)\right]} u(kh) \tag{7}$$

In the transfer functions defined by the one-sided **Z** transform one assumes zero

*l*¼0

*bik*0**A**½ � *<sup>μ</sup>i*ð Þ*<sup>k</sup>*

*aik*0**A**½ � *<sup>ν</sup>i*ð Þ*<sup>k</sup>*

*bik*0**A**½ � *<sup>μ</sup>i*ð Þ*<sup>k</sup>*

*<sup>k</sup>* is invertible, so for *k*<sup>0</sup> ¼ 0 one can write

*<sup>k</sup>* **u**ð Þ*k* (13)

*<sup>k</sup>* (14)

*<sup>k</sup>* (15)

*<sup>k</sup>* **u**ð Þ *kh* (16)

*<sup>N</sup>*½ � *<sup>μ</sup>P*ð Þ*<sup>k</sup> <sup>k</sup>***u**ð Þ *kh* (17)

*<sup>k</sup>* (18)

*<sup>k</sup>* **u**ð Þ *kh* (19)

*<sup>k</sup>* (20)

*<sup>k</sup> e kh* ð Þ (21)

initial conditions. Following this assumption equality (11) simplifies to

*<sup>k</sup>* **<sup>y</sup>**ð Þ¼ *<sup>k</sup>* <sup>X</sup>*mi*

**<sup>k</sup>** <sup>¼</sup> <sup>X</sup>*ni*

**<sup>k</sup>** <sup>¼</sup> <sup>X</sup>*mi*

*l*¼0

*l*¼0

*<sup>k</sup>* **<sup>y</sup>**ð Þ¼ *kh <sup>k</sup>*0*N*½ � *<sup>μ</sup>P*ð Þ*<sup>k</sup>*

h i�<sup>1</sup>

*<sup>k</sup>* <sup>¼</sup> <sup>0</sup>*D*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup> k* h i

**<sup>y</sup>**ð Þ¼ *kh* <sup>0</sup>*G*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup>* ,*μP*ð Þ*<sup>k</sup>*

**<sup>G</sup>***o*ð Þ¼ *kh* <sup>0</sup>*G*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup>* ,*μP*ð Þ*<sup>k</sup>*

fer function it is different by the real discrete variables. It relates discrete SISO

One considers a closed-loop system illustrated in **Figure 1**. Where a plant is

*<sup>k</sup> e kh* ð Þþ *KD*0Δ*<sup>ν</sup>*ð Þ*<sup>k</sup>*

systems by vectors and matrices related to its dimensions *k* þ 1∈ 0.

The classical PID controller output is desribed by three terms

*u kh* ð Þ¼ *KPe kh* ð Þþ *KI*0Δ�*μ*ð Þ*<sup>k</sup>*

*Remark 2.1*. Though the relation (19) looks similar to the classical discrete trans-

0

0 *N*½ � *<sup>μ</sup>*ð Þ*<sup>k</sup>*

*aik*0**A**½ � *<sup>ν</sup>i*ð Þ*<sup>k</sup>*

*Variable, Fractional-Order PID Controller Synthesis Novelty Method*

**k0D**½ � *<sup>ν</sup>***P**ð Þ **<sup>k</sup>**

**k0N**½ � *<sup>μ</sup>***P**ð Þ **<sup>k</sup>**

*<sup>k</sup>*0*D*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup>*

**<sup>y</sup>**ð Þ¼ *kh* <sup>0</sup>*D*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup> <sup>k</sup>*

0*G*½ � *<sup>ν</sup>P*ð Þ*<sup>k</sup>*

one gets similar to the transfer function description

Under assumption (12) *<sup>k</sup>*0*D*½ � *<sup>ν</sup>*ð Þ*<sup>k</sup>*

X*ni l*¼0

*DOI: http://dx.doi.org/10.5772/intechopen.95232*

Defining matrices

one gets

Denoting

or for simplicity

*2.4.1 VFO\_PID*

**5**

**2.4 VFO linear system description**

described by (19) where *e kh* ð Þ and *u kh* ð Þ.

where *mi* ≤*ni*, *νni*,*<sup>l</sup>* ð Þ*k* ≥*νni*,*l*�<sup>1</sup> ð Þ*k* ≥ ⋯*νi*,1ð Þ*k* ≥*νi*,0ð Þ¼ *k* 0, *μmi*,*<sup>l</sup>* ð Þ*<sup>k</sup>* <sup>≥</sup>*μmi*,*l*�1�1ð Þ*<sup>k</sup>* ≥ ⋯ *μi*,1ð Þ*k* ≥*μi*,0ð Þ*k* ≥0, *ai*,*<sup>l</sup>* and *bi*,*<sup>l</sup>* are constant coefficients for *l* ¼ 0, 1, ⋯, *ni* and *l* ¼ 0, 1, ⋯, *mi*, respectively. It is assumed that *a*0,*n*<sup>0</sup> ¼ 1.

According to the notation (5) Eq. (7) takes the form

$$\sum\_{l=0}^{n\_i} a\_{i,lk\_0} \mathbf{A}\_k^{\left[\mu\_{il}(k)\right]} \mathbf{y}(k) = \sum\_{l=0}^{m\_i} a\_{i,lk\_0} \mathbf{A}\_k^{\left[\mu\_{il}(k)\right]} \mathbf{u}(k) \tag{8}$$

The vector **u***j*ð Þ*k* satisfies the condition **u** *<sup>j</sup>*ð Þ¼ *k* **0***<sup>k</sup>* for *k*< *k*0. In the general solution of (8) to the assumed **<sup>u</sup>** *<sup>j</sup>*ð Þ*<sup>k</sup>* and initial conditions vector **<sup>y</sup>***<sup>i</sup>*,*k*0�<sup>1</sup> <sup>¼</sup> *yi*,*k*0�<sup>1</sup> *yi*,*k*0�<sup>2</sup> <sup>⋯</sup> h i<sup>T</sup> (T denotes the transposition) must be taken into account with �∞ ¼ *k*<sup>0</sup> <sup>0</sup> <0≤*k*<sup>0</sup> ≤*k*. Then, the infinite number of initial conditions (8) are formed in the following vector

$$\mathbf{y}\_{i,k\_0-1} = \begin{bmatrix} \mathbf{y}\_{i,k\_0-1} \\ \mathbf{y}\_{i,k\_0-2} \\ \vdots \end{bmatrix} \tag{9}$$

and the combined Eq. (8) is of the form

$$\begin{aligned} \left[\sum\_{l=0}^{n\_i} a\_{ij, lk\_0} \mathbf{A}\_k^{\left[\nu\_{i,l}(k)\right]} \quad \sum\_{l=0}^{n\_i} a\_{i, l-\infty} \mathbf{A}\_{k\_0-1}^{\left[\nu\_{i,l}(k)\right]}\right] \\ \times \left[\begin{matrix} \mathbf{y}\_i(kh) \\ \mathbf{y}\_{i,k\_0-1} \end{matrix}\right] = \sum\_{l=0}^{m\_i} b\_{ij, lk\_0} \mathbf{A}\_k^{\left[\mu\_{i,l}(k)\right]} \mathbf{u}(kh) \end{aligned} \tag{10}$$

or after simple transformation

$$\begin{split} \sum\_{l=0}^{n\_i} a\_{i,lk\_0} \mathbf{A}\_k^{\left[\nu\_{il}(k)\right]} \mathbf{y}(kh) &= \sum\_{l=0}^{m\_i} b\_{i,lk\_0} \mathbf{A}\_k^{\left[\mu\_{il}(k)\right]} \mathbf{u}(kh) \\ &- \sum\_{l=0}^{n\_i} a\_{i,l-\alpha} \mathbf{A}\_{k\_0-1}^{\left[\nu\_{il}(k)\right]} \mathbf{y}\_{i,k\_0-1} \end{split} \tag{11}$$

#### **2.3 Main assumptions**

To preserve the VFOBDE order one assumes that

$$\mathbf{1} + \sum\_{l=0}^{n\_i - 1} a\_{i,l} \neq \mathbf{0} \text{ for } i = \mathbf{1}, 2, \cdots, p \tag{12}$$

*Variable, Fractional-Order PID Controller Synthesis Novelty Method DOI: http://dx.doi.org/10.5772/intechopen.95232*

In the transfer functions defined by the one-sided **Z** transform one assumes zero initial conditions. Following this assumption equality (11) simplifies to

$$\sum\_{l=0}^{n\_i} a\_{ik\_0} \mathbf{A}\_k^{[\boldsymbol{\mu}\_i(k)]} \mathbf{y}(k) = \sum\_{l=0}^{m\_i} b\_{ik\_0} \mathbf{A}\_k^{[\boldsymbol{\mu}\_i(k)]} \mathbf{u}(k) \tag{13}$$

Defining matrices

$$\mathbf{a\_{k\_0}} \mathbf{D\_k^{[\nu\_\mathbf{P}(k)]}} = \sum\_{l=0}^{n\_i} a\_{ik\_0} \mathbf{A\_k^{[\nu\_i(k)]}} \tag{14}$$

$$\mathbf{k\_{0}}\mathbf{N\_{k}^{[\mu\_{\rm P}(k)]}} = \sum\_{l=0}^{m\_{i}} b\_{ik\_{0}} \mathbf{A}\_{k}^{[\mu\_{i}(k)]} \tag{15}$$

one gets

**2.2 Variable, fractional-order linear time-invariant difference equations**

*Control Based on PID Framework - The Mutual Promotion of Control and Identification…*

*<sup>k</sup> y kh* ð Þ¼ <sup>X</sup>*mi*

difference equation (VFODE). *h*> 0 denotes the sampling time.

X*ni l*¼0 *ai*,*<sup>l</sup> GL <sup>k</sup>*<sup>0</sup> <sup>Δ</sup>½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup>*

0, 1, ⋯, *mi*, respectively. It is assumed that *a*0,*n*<sup>0</sup> ¼ 1. According to the notation (5) Eq. (7) takes the form

*ai*,*lk*0**A**½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup>*

X*ni l*¼0

and the combined Eq. (8) is of the form

or after simple transformation

**2.3 Main assumptions**

**4**

X*ni l*¼0

P*ni l*¼0

> � **y***i* ð Þ *kh* **<sup>y</sup>***<sup>i</sup>*,*k*0�<sup>1</sup>

*ai*,*lk*0**A**½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup>*

To preserve the VFOBDE order one assumes that

1 þ *n* X*i*�1 *l*¼0

where *mi* ≤*ni*, *νni*,*<sup>l</sup>*

On the base of the Grünwald-Letnikov variable, fractional-order linear timeinvariant backward-difference the difference Eqs. (GL-VFOBE) for *i* ¼ 1, 2, ⋯, *p* and *j* ¼ 1, 2, ⋯, *q* representing discrete models of real dynamical systems or discrete control strategies are defined by the variable, fractional-order linear time-invariant

> *l*¼0 *bi*,*<sup>l</sup> GL <sup>k</sup>*<sup>0</sup> <sup>Δ</sup>½ � *<sup>μ</sup>i*,*l*ð Þ*<sup>k</sup>*

*l*¼0

*yi*,*k*0�<sup>1</sup> *yi*,*k*0�<sup>2</sup> ⋮

3 7

*ai*,*l*�∞**A**½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup> k*0�1

*bij*,*lk*0**A**½ � *<sup>μ</sup>i*,*l*ð Þ*<sup>k</sup>*

*bi*,*lk*0**A**½ � *<sup>μ</sup>i*,*l*ð Þ*<sup>k</sup>*

*ai*,*l*�<sup>∞</sup>**A**½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup>*

*<sup>k</sup>* **u**ð Þ *kh*

*<sup>k</sup>* **u**ð Þ *kh*

*<sup>k</sup>*0�<sup>1</sup> **<sup>y</sup>***<sup>i</sup>*,*k*0�<sup>1</sup>

*ai*,*<sup>l</sup>* 6¼ 0 for *i* ¼ 1, 2, ⋯, *p* (12)

The vector **u***j*ð Þ*k* satisfies the condition **u** *<sup>j</sup>*ð Þ¼ *k* **0***<sup>k</sup>* for *k*< *k*0. In the general solution

*ai*,*lk*0**A**½ � *<sup>μ</sup>i*,*l*ð Þ*<sup>k</sup>*

ð Þ*k* ≥*νni*,*l*�<sup>1</sup> ð Þ*k* ≥ ⋯*νi*,1ð Þ*k* ≥*νi*,0ð Þ¼ *k* 0, *μmi*,*<sup>l</sup>*

*μi*,1ð Þ*k* ≥*μi*,0ð Þ*k* ≥0, *ai*,*<sup>l</sup>* and *bi*,*<sup>l</sup>* are constant coefficients for *l* ¼ 0, 1, ⋯, *ni* and *l* ¼

*<sup>k</sup>* **<sup>y</sup>**ð Þ¼ *<sup>k</sup>* <sup>X</sup>*mi*

of (8) to the assumed **<sup>u</sup>** *<sup>j</sup>*ð Þ*<sup>k</sup>* and initial conditions vector **<sup>y</sup>***<sup>i</sup>*,*k*0�<sup>1</sup> <sup>¼</sup> *yi*,*k*0�<sup>1</sup> *yi*,*k*0�<sup>2</sup> <sup>⋯</sup>

Then, the infinite number of initial conditions (8) are formed in the following vector

2 6 4

P*ni l*¼0

<sup>¼</sup> <sup>X</sup>*mi l*¼0

*l*¼0

*<sup>k</sup>* **<sup>y</sup>**ð Þ¼ *kh* <sup>X</sup>*mi*

� X*ni l*¼0

� �

(T denotes the transposition) must be taken into account with �∞ ¼ *k*<sup>0</sup>

*aij*,*lk*0**A**½ � *<sup>ν</sup>i*,*l*ð Þ*<sup>k</sup> k*

" #

**<sup>y</sup>***<sup>i</sup>*,*k*0�<sup>1</sup> <sup>¼</sup>

*<sup>k</sup> u kh* ð Þ (7)

*<sup>k</sup>* **u**ð Þ*k* (8)

<sup>5</sup> (9)

ð Þ*<sup>k</sup>* <sup>≥</sup>*μmi*,*l*�1�1ð Þ*<sup>k</sup>* ≥ ⋯

h i<sup>T</sup>

<sup>0</sup> <0≤*k*<sup>0</sup> ≤*k*.

(10)

(11)

$$\_{k\_0}D\_k^{\left[\mu\_P(k)\right]} \mathbf{y}(kh) = \_{k\_0} \mathbf{N}\_k^{\left[\mu\_P(k)\right]} \mathbf{u}(kh) \tag{16}$$

Under assumption (12) *<sup>k</sup>*0*D*½ � *<sup>ν</sup>*ð Þ*<sup>k</sup> <sup>k</sup>* is invertible, so for *k*<sup>0</sup> ¼ 0 one can write

$$\mathbf{y}(kh) = \left[ \mathbf{o} \mathbf{D}^{[\nu\_{\mathcal{P}}(k)]\_{k}} \right]^{-1} \limits\_{\mathbf{0}} \mathbf{N}^{[\mu\_{\mathcal{P}}(k)]\_{k}} \mathbf{u}(kh) \tag{17}$$

Denoting

$${}\_{0}G\_{k}^{[\nu\_{P}(k)]} = \left[ {}\_{0}D\_{k}^{[\nu\_{P}(k)]} \right]\_{0}N\_{k}^{[\mu(k)]} \tag{18}$$

one gets similar to the transfer function description

$$\mathbf{y}(kh) = {}\_0\mathbf{G}\_k^{[\mu\_P(k), \mu\_P(k)]} \mathbf{u}(kh) \tag{19}$$

or for simplicity

$$\mathbf{G}\_o(kh) = \_0\mathbf{G}\_k^{[\nu\_P(k), \mu\_P(k)]} \tag{20}$$

*Remark 2.1*. Though the relation (19) looks similar to the classical discrete transfer function it is different by the real discrete variables. It relates discrete SISO systems by vectors and matrices related to its dimensions *k* þ 1∈ 0.

### **2.4 VFO linear system description**

One considers a closed-loop system illustrated in **Figure 1**. Where a plant is described by (19) where *e kh* ð Þ and *u kh* ð Þ.

### *2.4.1 VFO\_PID*

The classical PID controller output is desribed by three terms

$$u(kh) = K\_P e(kh) + K\_{I0} \Delta\_k^{-\mu(k)} e(kh) + K\_{D0} \Delta\_k^{\nu(k)} e(kh) \tag{21}$$

• **y***o*ð Þ *kh* - a plant output signal vector,

*DOI: http://dx.doi.org/10.5772/intechopen.95232*

• **e**ð Þ *kh* - a closed-loop system error signal,

A system error is evaluated by the formula

Eq. (29) simplifies to

**7**

• **y**ð Þ *kh* - a closed-loop system output signal vector,

*Variable, Fractional-Order PID Controller Synthesis Novelty Method*

**<sup>e</sup>**ð Þ¼ *kh* ½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* **<sup>H</sup>**ð Þ *kh* �<sup>1</sup>

**3. Variable, fractional-order PID controller synthesis**

�½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* **<sup>H</sup>**ð Þ *kh* �<sup>1</sup>

In the synthesis of the classical PID controller there are three parameters to evaluate. Namely, *K*,*KI*,*KD* known as the proportional, integral and differential gains. In the fractional-order PID controllers there are two additional parameters: the differentiation order *ν*ð Þ *kh* ∈ <sup>þ</sup> and the integration one �*μ*ð Þ *kh* ∈ þ. In the variable, fractionalorder PID controller the mentioned orders are generalized to functions. This means that there are three constant coefficients and two discrete variable functions to find

In the rejection of the external disturbation one can assume that **r**ð Þ¼ *kh* **0** so

Usually the sensor matrix **H**ð Þ *kh* is treated as constant, by assumption that sensors do not introduce its own dynamics to the system. Hence, **H**ð Þ¼ *kh* **H** ¼ const. It may

The optimal parameters (29) are evaluated due to the assumed optymality criterion. The most popular is so called ISE one (Integral of the Squared Error) or in

> *k* X*max i*¼0

In the proposed VFOPID controller synthesis method with partially intuitive and supported by closed-loop systems synthesis experience the classical optimisation due to the performance criterion (32) is performed. The pre-defined differentiation

<sup>¼</sup> **<sup>d</sup>**ð Þ *kh <sup>T</sup>*½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* �*<sup>T</sup>*½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* �<sup>1</sup>

*e ih* ð Þ<sup>2</sup>

*<sup>h</sup>* <sup>¼</sup> **<sup>e</sup>**ð Þ *kh* <sup>T</sup>

*ν*ð Þ *kh* ≥0 (34)

**<sup>e</sup>**ð Þ¼� *kh* ½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* **<sup>H</sup>**ð Þ *kh* �<sup>1</sup>

be assumed that **H** ¼ *h*0**1***<sup>k</sup>* or further, for *h*<sup>0</sup> ¼ 1, formula (30) takes a form

the discrete-system case: Sum of the Squared Error (SSE).

*SSE KP*,*KI* ½ ,*KD*, *ν*ð Þ *kh* , *μ*ð Þ *kh* � ¼

*SSE KP*, *KI* ½ � ,*KD*, *ν*ð Þ *kh* , *μ*ð Þ *kh*

and integration order functions orders are as follows

Substitution of (31) into (32) gives

**<sup>e</sup>**ð Þ¼� *kh* ½ � **<sup>1</sup>***<sup>k</sup>* <sup>þ</sup> **<sup>G</sup>***o*ð Þ *kh* **<sup>C</sup>**ð Þ *kh* �<sup>1</sup>

**r**ð Þ *kh*

*KP*,*KI*,*KD*, *ν*ð Þ *kh* , *μ*ð Þ *kh* (29)

**<sup>H</sup>**ð Þ *kh* **<sup>d</sup>**ð Þ *kh* (28)

**H**ð Þ *kh* **d**ð Þ *kh* (30)

**d**ð Þ *kh* (31)

**e**ð Þ *kh h* (32)

**<sup>d</sup>**ð Þ *kh* (33)

**Figure 1.** *Closed-loop system.*

and in the convention proposed above as

$$\mathbf{u}(kh) = K\_P \mathbf{1}\_k \mathbf{e}(kh) + K\_{D0} \mathbf{G}\_k^{[\nu\_C(k)]} \mathbf{e}(kh) + K\_{I0} \mathbf{G}\_k^{[-\mu\_C(k)]} \mathbf{e}(kh) \tag{22}$$

which may be expressed as

$$\mathbf{u}(kh) = \left[K\_P \mathbf{1}\_k + K\_{D0} \mathbf{G}\_k^{[\iota\_C(k)]} + K\_{I0} \mathbf{G}\_k^{[-\mu\_C(k)]}\right] \mathbf{e}(kh) \tag{23}$$

where *νC*ð Þ*k* , *μC*ð Þ*k* ≥0 and controlling and error signals are denoted as **u**ð Þ*k* and **e**ð Þ*k* , respectively. Then, denoting.

*Remark 2.2*. The plant may be described by classical integer order, fractional or even variable, fractional - order difference equations. The matrix - vector description used makes it possible.

$$\mathbf{C}(kh) = K\_P \mathbf{1}\_k + K\_{D0} \mathbf{G}\_k^{[\boldsymbol{\nu}\_C(k)]} + K\_{I0} \mathbf{G}\_k^{[-\boldsymbol{\mu}\_C(k)]} \tag{24}$$

one gets a VFOPID controller transfer function-like description

$$\mathbf{u}(kh) = \mathbf{C}(kh)\mathbf{e}(kh) \tag{25}$$

To simplify the description one assumes a sensor matrix as

$$\mathbf{H}(kh) = \mathbf{1}\_k\tag{26}$$

The closed-loop system is presented in **Figure 1** from which one gets the following relations

$$\begin{split} \mathbf{y}(kh) &= \left[\mathbf{1}\_{k} + \mathbf{G}\_{o}(kh)\mathbf{C}(kh)\mathbf{H}(kh)\right]^{-1}\mathbf{G}\_{o}(kh)\mathbf{C}(kh)\mathbf{r}(kh) \\ &+ \left[\mathbf{1}\_{k} + \mathbf{G}\_{o}(kh)\mathbf{C}(kh)\mathbf{H}(kh)\right]^{-1}\mathbf{d}(kh) \end{split} \tag{27}$$

where


*Variable, Fractional-Order PID Controller Synthesis Novelty Method DOI: http://dx.doi.org/10.5772/intechopen.95232*


A system error is evaluated by the formula

$$\begin{aligned} \mathbf{e}(kh) &= \left[\mathbf{1}\_k + \mathbf{G}\_o(kh)\mathbf{C}(kh)\mathbf{H}(kh)\right]^{-1}\mathbf{r}(kh) \\ &- \left[\mathbf{1}\_k + \mathbf{G}\_o(kh)\mathbf{C}(kh)\mathbf{H}(kh)\right]^{-1}\mathbf{H}(kh)\mathbf{d}(kh) \end{aligned} \tag{28}$$
