Suboptimal Regulation of a Class of Bilinear Interconnected Systems with Finite-Time Sliding Planning Horizons

Abstract

This paper focuses on the suboptimization of a class of multivariable discrete-time bilinear systems consisting of interconnected bilinear subsystems with respect to a linear quadratic optimal regulation criterion which involves the use of state weighting terms only. Conditions which ensure the controllability of the overall system are given as a previous requirement for optimization. Three transformations of variables are made on the system equations in order to implement the scheme on an equivalent linear system. This leads to an equivalent representation of the used quadratic performance index that involves the appearance of quadratic weighting terms related to both transformed input and state variables. In this way, a Riccati-matrix sequence, allowing the synthesis of a standard feedback control law, is obtained. Finally, the proposed control scheme is tested on realistic examples.

1. Introduction

Dynamic bilinear systems have received great attention by researchers in the last decades, from the classical works of Anderson and Moore [1], Feldbaum [2], Tarn [3], or Tarn et al. [4], to the recent works of Al-Baiyat [5], Kotta et al. [6] or Garrido et al. [7]. They form a transitional class between the linear and the general nonlinear systems. The importance of such systems is well-known in diverse areas in which bilinear systems appear naturally for a variety of dynamic processes such as those taking place in electrical systems [8], industrial processes [9], power plants [10], chemical processes [11], nuclear fusion [12], biomedical applications [13], information management [14], mechanical systems [15], and aerospace and avionics [16]. which can not be satisfactorily modeled under the assumption of linearity. On the other hand, the design of controllers for bilinear systems has been an area of major research during the recent years. This growing interest in practical applications requires the development of suitable algorithms for control problems associated with bilinear systems. Many of the results obtained rely on the optimization theory with an appropriate performance index, as it is the case of the proposed control scheme.

Well-known works as Goka et al. [17] and Tarn et al. [18] stated the controllability conditions for standard classes of discrete bilinear systems. In the first paper, an equivalent system description was derived with the equivalent input being dependent on the products of the state and the original input. Since the equivalent feed-forward loop is of a linear nature, the analysis becomes greatly simplified. In the second paper, the uniform controllability of such systems using bounded input was studied. On the other hand, it is well known that the systems are often interconnected and, in many cases, several dynamic subsystems can be distinguished for analysis purposes. This sometimes takes place in computer communication, transportation networks, control and power stations, and so on. There exists important literature dealing with both the associate multivariable and the large-scale problems [19–25]. For instance, standard adaptive control techniques were applied in [26] to compensate the undesirable deviations of the process parameters from their nominal values, being the overall process modeled as a bilinear system prior to the application of the adaptive scheme.

This paper reports suboptimal optimization techniques which are applied to bilinear models. Such models can be considered as direct extensions to the linear continuous interconnected systems stated by the work of Ramakrishna and Viswanadham [27]. A class of invariant discrete-time multivariable bilinear systems with interconnection subsystems is studied. First, the system is shown to be equivalently described by the linear feed-forward multivariable structure with multiplicative inputs including a deterministic disturbance vector. Also, other interpretations are stated. As a previous requirement for optimization, controllability results for the overall systems performance are given by extending those ones given in [17, 18, 27–31]. This is achieved through the above equivalent linear system with the use of centralized control methods. A central coordinator is supposed to be available for the local controllers to supply information to each control station from the remaining ones. That control technique is applied to the optimization of the aforementioned system class. The used performance index involves state-weighting terms only. An explicit solution of Riccati type (which is unusual for such an optimization criterion) is found out as the suboptimal solution by using manipulations on the input/state variables of the problem statement [1, 32]. The importance of this strategy arises from problems where constrains on the input rather than input weights are introduced in the optimization criterion, what leads to a feedback-type control law.

The paper is organized as follows. In Section 2, the problem statement for the class of the interconnected invariant discrete-time bilinear systems is given. Also, controllability necessary conditions obtained by the decomposition techniques are studied. Section 3 is devoted to the development of suboptimal strategies for the regulation criterion, which involves quadratic terms of the state variables. Three transformations of variables are used to turn the loss function into another one, which involves state and input quadratic terms leading to a Riccati-type suboptimal solution. The control strategies may be considered as extensions of those ones used in [1, 29] or [25], for the linear continuous problem. The formulation is accomplished from a centralized viewpoint. In Section 4, two simulated examples are presented. The first one consists of a simulated numerical example; and the second one consists of a real regulation problem of two amplidynes acting on a DC-motor. In the second example, it is not necessary to consider constant current or voltage as it is usual in DC-motor control, with their product—namely, the electric par—a bilinear term describing the resulting dynamical system. Finally, conclusions end the paper. Appendix A is related to the meaning explanation of certain equivalent interpretations of the overall system structure. Appendix B concerns with the development of the transformations of variables involved in the redefinition of the performance index. By convenience, equations and results from the appendices are sometimes invoked in the main body of the paper order not to repeat mathematical material.

2. On the System Structure and Its Controllability Properties

2.1. Controllability of a Class of Bilinear Systems

In [17], the following typical class of bilinear systems was considered: 𝑥(𝑘+1)=𝐴+𝑢(𝑘)𝐵𝑥(𝑘)+𝐶𝑢(𝑘);𝑘=0,1,2,…,𝑥(⋅)∈𝑅^𝑛;𝑢(⋅)∈𝑅.(2.1) Equation (2.1) can equivalently be described, if rank (𝐵,𝐶)=1 and then rank (𝐵)=1 if 𝐶 and ℎ are both nonzero, with 𝐵=𝐶ℎ^𝑇 (𝐶,ℎ being unique 𝑛-vectors within a scalar factor), by the linear feed-forward system with multiplicative input: 𝑥(𝑘+1)=𝐴𝑥(𝑘)+𝐶𝜈(𝑘);𝜈(𝑘)≡⎧⎪⎨⎪⎩𝑢(𝑘)ℎ^𝑇𝑥(𝑘);if𝐶=0,𝑢(𝑘)[ℎ^𝑇𝑥(𝑘)+1];if𝐶≠0(2.2) with 𝑘=0,1,2,…; and where 𝑥(𝑘) and 𝑢(𝑘) are the state and the control at time 𝑘, respectively; and 𝐴 and 𝐵 are real constant coefficient matrices of compatible dimensions.

Note that for arbitrary ranks of 𝐵 or (𝐵,𝐶) in (2.1), the decomposition 𝐵=∑^𝑝_𝑖=1𝑒_𝑖𝑏^𝑇_𝑖 can be used, with 𝑒_𝑖 and 𝑏^𝑇_𝑖   being, respectively, the 𝑖th unity vector (i.e., the 𝑖th component is unity and the remaining ones are zero) and the 𝑖th row vector of 𝐵. The controllability of the bilinear system was related in that paper to that of the linear system (2.2) as follows.

Theorem 2.1. (summary of results in [17]). For the equivalent systems (2.1)-(2.2), the following controllability results hold.
(i)If the bilinear system (2.1) is controllable, then the pair (𝐴,𝐶) in (2.2) is controllable.(ii)If (𝐴^𝑇,ℎ) is not controllable, then the homogeneous bilinear system is not controllable.(iii)If (𝐴,𝐶) and (𝐴^𝑇,ℎ) are both controllable pairs and if ℎ^𝑇𝐶, ℎ^𝑇𝐴⁻¹𝐶≠0, then the homogenous bilinear system is controllable in 𝑅^𝑛 outside the origin.(iv)The inhomogeneous bilinear system is not controllable in 𝑅^𝑛 if(a)rank [controllability matrix of (𝐴^𝑇,ℎ)] = rank [controllability matrix of (𝐴^𝑇,ℎ);1]; (1≡𝑛-vector with every element being unity),(b)∑^𝑛_𝑖=1𝑎_𝑖=1, with a_(⋅) being the elements of the last row of the 𝐴-matrix in the canonical state-space representation of system (2.2).(v)The inhomogeneous bilinear system is controllable in 𝑅^𝑛, if (𝐴,𝐶) is a controllable pair, ℎ^𝑇𝐶, ℎ^𝑇𝐴⁻¹𝐶≠0, and at least one of the conditions (iv-a) or (iv-b) is not fulfilled.
Although these results are restricted to a special class of discrete bilinear systems, many natural discrete bilinear systems [10, 15, 33–35] do satisfy these assumptions. More general controllability conditions for bilinear systems can be found in [18, 36].

Remark 2.2. More general results were derived in [30] for the case in which 𝐴⁻¹ does not exist. Such a paper studies the controllability conditions for three cases of violation of the conditions given in [37] (namely, (𝐴,𝐶) and (𝐴^𝑇,ℎ) are controllable pairs and 𝑘₀=g.c.d{𝑖/ℎ^𝑇𝐴^𝑖−1𝐶≠0,0<𝑖≤𝑛²}=1) (g.c.d denotes the greatest common divisor). Such cases include the following:
(a)rank ⌈(𝐴,𝐶)=𝑛, rank Ω(ℎ^𝑇,𝐴)<𝑛;(b)rank ⌈(𝐴,𝐶)<𝑛, rank Ω(ℎ^𝑇,𝐴)=𝑛; and(c)rank ⌈(𝐴;𝐶)<𝑛, rank Ω(ℎ^𝑇,𝐴)<𝑛, with ⌈(𝐴,𝐶)≡[𝐶,𝐴𝐶,…,𝐴^𝑛−1𝐶] and Ω(ℎ^𝑇,𝐴)≡[ℎ,𝐴^𝑇ℎ,

2.2. Models for a Class of Bilinear Systems with Interconnected Subsystems

In this paper, the following multivariable invariant discrete-time structure of bilinear systems (…,𝐴^{(𝑛−1)𝑇}ℎ]^𝑇.) of homogeneous type is considered 𝑆_𝑖 with the interconnection bilinear subsystems (𝑆_𝑖∶𝑥_𝑖(𝑘+1)=𝐴_𝑖+𝑢_𝑖(𝑘)𝐵_𝑖𝑥_𝑖(𝑘)+𝐷_𝑖𝑧_𝑖(𝑘),(2.3)𝑦_𝑖(𝑘)=𝐶_𝑖𝑥_𝑖(𝑘);𝑖∈𝐼;𝐼=1,2,…,𝑝;𝑘=0,1,2,…(2.4)), which include linear and bilinear coupling terms, given by 𝐻_i where 𝐻_𝑖∶𝑧_𝑖(𝑘+1)=𝑀_𝑖+𝑢_𝑖(𝑘)𝑃_𝑖𝑧_𝑖(𝑘)+^𝑝_𝑗=1𝐿_𝑖𝑗+𝑢_𝑖(𝑘)𝑀_𝑖𝑗𝑦_𝑗(𝑘);𝑖∈𝐼;𝑘=0,1,2,…,(2.5) are, respectively, the state vector of the 𝑥_𝑖(𝑘)∈𝑅^𝑛𝑖,𝑧_𝑖(𝑘)∈𝑅^𝑎𝑖,𝑦_𝑖(𝑘)∈𝑅^𝑟𝑖,𝑢_𝑖(𝑘)∈𝑅-subsystems, the state vectors of the interaction 𝑆_𝑖-subsystems, the corresponding output and input vectors at time 𝐻_𝑖, and the different coefficient matrices of appropriate dimensions (see Figure 1).

Figure 1 

The class of multivariable discrete invariant-time bilinear subsystems with the interaction bilinear subsystems.

The linear continuous models of [27] for Figure 1 are more general than those ones above since additional measurements from the interconnection subsystems are used to generate the inputs to the plant. Structures of this type arise from practical systems such as countercurrent heat exchangers having weak nonlinearities modeled as bilinear.

Now, the extended state vector of the overall system is defined as 𝑘 all 𝑥(𝑘)≡𝑥^𝑇₁(𝑘),𝑥^𝑇₂(𝑘),…,𝑥^𝑇_𝑝(𝑘)^𝑇∈𝑅^𝑛+𝑎,where𝑥_𝑖(𝑘)≡𝑥^𝑇_𝑖(𝑘),𝑧^𝑖𝑇(𝑘)^𝑇∈𝑅^{𝑛𝑖+𝑎𝑖},(2.6) with 𝑖∈𝐼 and the extended output vector is defined as 𝑛≡∑^𝑝_𝑗=1𝑛_𝑗,𝑎≡∑^𝑝_𝑗=1𝑎_𝑗, with 𝑦(𝑘)≡[𝑦^𝑇₁(𝑘),𝑦^𝑇₂(𝑘),…,𝑦^𝑇_𝑝(𝑘)]^𝑇∈𝑅^𝑟, Thus, the following multivariable structure, equivalent to (2.3)–(2.5), is deduced: 𝑟≡∑^𝑝_𝑗=1𝑟_𝑗. where the matrices of parameters have the following structures: 𝑥(𝑘+1)=𝐴𝑥(𝑘)+^𝑝_𝑖=1𝑢_𝑖(𝑘)𝐵_𝑖𝑥(𝑘);𝑘=0,1,2,3,…,(2.7)𝑦(𝑘)=𝐶𝑥(𝑘);𝑘=0,1,2,3,…,(2.8) for all 𝐴=𝐴^𝑇₁,𝐴^𝑇₂,…,𝐴^𝑇_𝑝^𝑇;𝐴_𝑖=𝐴_𝑖1,𝐴_𝑖2,…,𝐴_𝑖𝑖,…,𝐴_𝑖𝑝𝐴_𝑖𝑖=[𝐴_𝑖𝐷_𝑖𝐿_𝑖𝑖𝐶_𝑖𝑀_𝑖]𝐴_𝑖𝑗=[00𝐿_𝑖𝑗𝐶_𝑗0](for𝑗≠𝑖)𝐵_𝑖=0^𝑖−1,…,0,𝐵^𝑇_𝑖,0^{𝑝−1},…,0^𝑇;𝐵_𝑖=𝐵_𝑖1,𝐵_𝑖2,…,𝐵_𝑖𝑝𝐵_𝑖𝑖=[𝐵_𝑖0𝑀_𝑖𝑖𝐶_𝑖𝑃_𝑖]𝐵_𝑖𝑗=[00𝑀_𝑖𝑗𝐶_𝑗0](for𝑗≠𝑖)𝐶=𝐶^𝑇₁,𝐶^𝑇₂,…,𝐶^𝑇_𝑝;𝐶_𝑖=0^𝑖−1,…,0,𝐶_𝑖,0,0^{𝑝−1},…,0(2.9)

Remark 2.3. In general, the number of subsystems 𝑗,𝑖∈𝐼. and the number of interaction subsystems 𝑆_(⋅) are not equal. Thus, assume that 𝐻_(⋅) and 𝑝₁ are, respectively, the number of subsystems 𝑝₂(𝑝₁≠𝑝₂) and the number of subsystems 𝑆_𝑡,𝑡∈𝐼₁={1,2,3,…,𝑝₁}, Then, an extended state vector 𝐻_𝑠,𝑠∈𝐼₂={1,2,…,𝑝₂}. can be obtained similarly, if 𝜉(⋅) and 𝑝=max(𝑝₁𝑝₂)
Taking into account the dimensions and the structural zeros of the 𝐼=𝐼₁𝑈𝐼₂.-matrices, they can be decomposed into a sum of dyads as 𝐵_𝑖 where 𝐵_𝑖=^{𝑛𝑖+𝑎𝑖}_𝑗=1𝑒_{𝑗+𝑙𝑖−1}𝑏^𝑇_{𝑗+𝑙𝑖−1}(𝑖);forall𝑖∈𝐼,(2.10) is the unity 𝑒_𝑗 vector, 𝑗th is the 𝑏^𝑇_𝑗(𝑖) vector of the 𝑗th-matrix, 𝐵_𝑖, and the lower subscripts 𝑖∈𝐼 are related to the first nonstructural zero rows of the 1_𝑗,all𝑗∈𝐼,-matrix. If they have additional (nonstructural) rows being zero, the corresponding 𝐵_𝑖-vectors can be neglected in (2.10). From (2.7), it is clear that the rows of the 𝑒_(⋅)-matrices have the following structures: 𝐵_i𝑏^𝑇_𝑗(𝑖)≡⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩[0,…,^{(𝑖−1)block}0,𝑒^𝑇_{𝑗−𝑙𝑖+1}𝐵_𝑖,0,…,^{(𝑝−𝑖)blocks}0];if1_𝑖≤𝑗<1_𝑖+𝑛_𝑖,[𝑒^𝑇_{𝑗−𝑙𝑖−𝑛𝑖+1}𝑀_𝑖1𝐶₁,0,…,𝑒^𝑇_{𝑗−𝑙𝑖.𝑛𝑖+1}𝑀_{𝑖(𝑖−1)}𝐶_𝑖−1,^{(𝑖−1)blocks}0,𝑒^𝑇_{𝑗−𝑙𝑖−𝑛𝑖+1}𝑃_𝑖,_{(𝑝−1)blocks}𝑒^𝑇_{𝑗−𝑙𝑖−𝑛𝑖+1}𝑀_𝑖(𝑖+1)𝐶_𝑖+1,0,…,𝑒^𝑇_{𝑗−𝑙𝑖−𝑛𝑖+1}𝑀_𝑖𝑝𝐶_𝑖,0];if𝑙_𝑖+𝑛_𝑖𝑖≤𝑗<1_𝑖+𝑛_𝑖+𝑎_𝑖,(2.11) Then one deduces from (2.11) that forall𝑖∈𝐼.𝑏^𝑇_{𝑗+1𝑖−1}(𝑖)𝑥≡⎧⎪⎪⎨⎪⎪⎩𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘);if1≤𝑗≤𝑛_𝑖,𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖+1}𝑦_𝑙(𝑘);if𝑛_𝑖+1≤𝑗≤𝑛_𝑖+𝑎_𝑖,0;if𝑛_𝑖+𝑎_𝑖+1≤𝑗≤𝑛+𝑎,(2.12), when with𝑘=0,1,2,3,…, 𝑏^𝑇_𝑗(𝑖), 𝑝^𝑇_𝑗(𝑖) denote, respectively, the 𝑚^𝑇_𝑗(𝑖,1) vector row of the matrices 𝑗th From (2.7) through (2.12) one obtains 𝐵_𝑖,𝑃_𝑖,and𝑀_𝑖,forall𝐼,1∈𝐼. where the sequences of scalars 𝑥(𝑘+1)=𝐴𝑥(𝑘)+^𝑝_𝑖=1^{𝑛𝑖+𝑎𝑒}_𝑗=1𝑒_{𝑗+𝑙𝑖−1}𝜈_𝑗𝑖(𝑘);𝑘=0,1,2,…,(2.13) are defined by 𝜈_𝑗𝑖(𝑘)𝜈_𝑗𝑖(𝑘)≡⎧⎪⎨⎪⎩𝑢_𝑖(𝑘)𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘);if1≤𝑗≤𝑛_𝑖,𝑢_𝑖(𝑘)𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘);if𝑛_𝑖+1≤𝑗≤𝑛_𝑖+𝑎_𝑖,(2.14) are considered as equivalent (primary) inputs for the equivalent linear systems (2.13) to be determined prior to the generation of the forall𝑖∈𝐼,𝑘=0,1,2,…,-sequence in the design implementation.

2.3. Obtaining the 𝑢_(⋅)(⋅)-Sequences from the 𝑢_𝑖(𝑘)-Sequences

Equations (2.13) and (2.14) constitute a feed-forward linear equivalent system with multiplicative inputs. But the 𝛾_𝑗𝑖(𝑘)𝑝-inputs cannot, in general, be determined from the 𝑢_𝑖(𝑘) equivalent (𝑛+𝑎)>𝑝-inputs. Then, it is necessary to take a set of p variables between the last ones, which will be optimized independently. Let us call 𝜈_𝑗𝑖(𝑘) Thus, let us define the sets

𝜈_𝑗𝑖(𝑘),all𝑖∈𝐼. such that for the 𝑁_𝑖=𝑁⁻_𝑖𝑈𝑁⁺_𝑖,𝑁⁻_𝑖=1,2,…,𝑛_𝑖,𝑁⁺_𝑖=𝑛_𝑖+1,𝑛_𝑖+2,…,𝑛_𝑖+𝑎_𝑖,all𝑖∈𝐼,(2.15) current sampling instant, 𝑘th and 𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)≠0,if𝑥_𝑖(𝑘)≠0,if𝑗_𝑖∈𝑁⁻_𝑖,𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦_𝑖(𝑘)≠0,if𝑗_𝑖∈𝑁⁺_𝑖,(2.16)

Then from (2.14), one has the auxiliary inputs 𝑘=0,1,2,…. with the real inputs for the nonsaturated case being𝜈_𝑗𝑖(𝑘)≡⎧⎪⎪⎪⎨⎪⎪⎪⎩𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)^−𝑖𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘)𝜈_𝑗𝑖(𝑘);𝑗∈𝑁⁻_𝑖;𝑖∈𝐼,𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)⁻¹×𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)𝜈_𝑗𝑖(𝑘),𝑗∈𝑁⁺_𝑖;𝑖∈𝐼;𝑘=0,1,2,3,…,(2.17)

Remark 2.4. Note that the minimum energy control (i.e., 𝑢_𝑖(𝑘)≡⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)⁻¹𝜈_𝑗𝑖(𝑘);if𝑗_𝑖∈𝑁_𝑖,ifthereexistssome𝑗_𝑖∈𝑁⁻_𝑖suchthat𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)≠0,𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦_𝑖(𝑘)⁻¹𝜈_𝑗𝑖(𝑘);if𝑗_𝑖∈𝑁_𝑖,ifthereexistssome𝑗_𝑖∈𝑁⁺_𝑖suchthat𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)≠0,0;if𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)=0forall𝑗_𝑖∈𝑁⁻_𝑖andalso𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)=0forall𝑗_𝑖∈𝑁⁺_𝑖,all𝑖∈𝐼;𝑘=0,1,2,3,….(2.18) some 𝑢_𝑖(⋅)=0,) has been chosen when the system becomes insensitive to the corresponding control component. In fact, any arbitrary control could be applied.

Remark 2.5. Since there are 𝑖∈𝐼 equivalent inputs, deals with an over determined problem when the (𝑛+𝑎) inputs are solved from them. Their solution cannot generally be satisfied exactly, and there exist many possible ways of defining the “best” approximate solution if there are 𝑝<𝑛+𝑎 equivalent nonzero inputs. Also, note that the decomposition of the 𝑝<𝑠≤𝑛+𝑎-matrices (here, into a sum of dyads, as shown in (2.10)) is not unique. Then, the least-squares technique (which is known to give the best approximation) can be considered by taking into account the overall control term in (2.13) instead of the control itself. Namely, 𝐵_𝑖 where 𝑢_𝑖(𝑘)≡⎧⎪⎨⎪⎩𝜈^𝑇_𝑖(𝑘)𝐵_𝑖𝑥(𝑘)𝑥^𝑇(𝑘)𝐵^𝑇_𝑖𝐵_𝑖𝑥(𝑘);if𝐵_𝑖𝑥(𝑘)≠0,0,otherwise,all𝑖∈𝐼;𝑘=0,1,2,3,…,(2.19)
Note that in general, the (𝜈_𝑖(𝑘)≡[0,…,0^{∑𝑖−1𝑗=1(𝑛𝑗+𝑎𝑗)zeros},𝜈_𝑙𝑖(𝑘),…,𝜈_{𝑛𝑖,𝑖}(𝑘),𝜈_{𝑛𝑖+1,𝑖}(𝑘),…,𝜈_{𝑛𝑖+𝑎𝑖,𝑖}(𝑘),0,…,^{∑𝑝𝑗=𝑖+1(𝑛𝑗+𝑎𝑗)zeros}0]^𝑇∈𝑅^𝑛+𝑎,all𝑖∈𝐼;𝑘=0,1,2,3,….(2.20)) vector has only nonzero components. Equation (2.19) can be rewritten as 𝑘𝑢_𝑖(𝑘)≡[^{𝑛𝑖+𝑎𝑖}_𝑗=1{𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗(𝑖)+𝑝^𝑇_{𝑗−𝑛𝑖}0(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)𝛿_𝑗(𝑖)}²]⁻¹×[^{𝑛𝑖+𝑎𝑖}_𝑗=1𝜈_𝑗𝑖(𝑘){𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗(𝑖)+𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)𝛿_𝑗(𝑖)}],(2.21)

2.4. Two Alternative Interpretation Schemes of the Equivalent Linear System

Equation (2.18) leads to the following equivalences to the linear system (2.13): where𝛿(𝑖)≡1−𝛿_𝑗(𝑖),with𝛿_𝑗(𝑖)=1if1≤𝑗≤𝑛and0,if𝑛_𝑖+1≤𝑗≤𝑛_𝑖+𝑎_𝑖all𝑖∈𝐼;𝑘=0,1,2,…. where the 𝑥(𝑘+1)=𝐴𝑥(𝑘)+^𝑝_𝑖=1𝑐_{𝑗+𝑙𝑖−1}(𝑘)𝜈_𝑗𝑖(𝑘)(2.22)=𝐴𝑥(𝑘)+^𝑝_𝑖=1𝑒_{𝑗+𝑙𝑖−1}𝜈_𝑗𝑖(𝑘)+𝜉(𝑘),(2.23)-vectors (𝑛+𝑎) and 𝑐_(⋅)(⋅) are defined in Appendix A.

Equations (2.22)-(2.23) provide two alternative interpretations of (2.13). Representation (2.22) involves an extended state-dependent control vector, while (2.23) includes a known deterministic input-dependent vector. These equivalent representations will lead to two different suboptimal control designs in the next section.

2.5. Controllability of the Bilinear System with Interconnections

The controllability outside the origin of the bilinear system (2.3)–(2.5) can be studied by splitting the problem as follows:

(i)controllability of the feed-forward linear system,(ii)capability of the 𝜉(⋅)-sequences of actually acting upon the equivalent 𝑢_(⋅)(⋅)-sequences.The following result on controllability conditions is useful as a previous requirement for optimization.

Theorem 2.6. (controllability of the bilinear system involving interconnections). Controllability of the bilinear system (2.3)–(2.5) is ensured on an interval 𝜈_(⋅)(⋅), [𝑘^,𝑘^]∩𝑅⁺, if the following set of conditions holds.
(i)The controllability grammian of the overall system is positive definite in 𝑘^=𝑘^+𝑚𝑎𝑥_{1≤𝑖≤𝑝}(𝑛_𝑖+𝑎_𝑖)(ii)[𝑘^,𝑘^]. with max[rank(∑^𝑝_𝑙=1𝑀_𝑠1𝐶₁)]≥max[rank(𝐿_𝑖𝑗𝐶_𝑗)](iii)All the pairs of block diagram, 1≤𝑠≤𝑝,1≤𝑖≤𝑝and1≤𝑗≤𝑝., are nonsingular.(iv)[(𝐴^𝑇_𝑖,𝑀^𝑇_𝑖),𝐵_𝑖],𝑖∈𝐼, (𝑏^∗𝑇_𝑗𝑖(𝑖)𝐶^∗_𝑖(𝑘)≠0 diagram 𝑏^∗𝑇_𝑗𝑖(𝑖), (𝐴⁻¹_𝑖,𝑀⁻¹_𝑖)𝐶^∗_𝑖(𝑘)≠0, all𝑖∈𝐼, where all𝑘^≤𝑘≤𝑘^ are, respectively, the subvectors of the matrices 𝐶^∗_(⋅)and𝑏^∗_(⋅) in (2.3)-(2.4) associated with each subsystem 𝐶_(⋅)and𝐵_(⋅).

Sketch of the Proof
Since the pair 𝑆_(⋅) is controllable for all arbitrary matrix [𝐴+𝐵_(⋅)𝐺_(⋅),𝐵_(⋅)] if the pair 𝐺_(⋅) is controllable [29, 38], then the controllability of the pair [𝐴,𝐵_(⋅)] implies the controllability of the pair [𝐴^𝑑,𝐵_(⋅)] if there exists a [𝐴^𝑑+𝐴^∘,𝐵_(⋅)],-matrix such that 𝐺_(⋅) with 𝐴^∘=𝐵_(⋅)𝐺_(⋅), 𝐴^∘=𝐴−𝐴^𝑑  block diagram 𝐴^𝑑=. For the equality above relating ((𝐴₁,𝑀₁),(𝐴₂,𝑀₂),…,(𝐴_𝑝,𝑀_𝑝)) and 𝐴^∘,𝐵_(⋅), to hold, the Rouché-Fröbenius theorem implies that condition (ii) must be fulfilled according to (2.9) and (2.11). Also, since controllability grammian of 𝐺_(⋅) on [𝐴^𝑑,𝐵_(⋅)] block diagram controllability grammian [[block diagram [𝑘^,𝑘^]≡],[block diagram (𝐴₁,𝑀₁),𝐵₁], (𝐴₂,𝑀₂),𝐵₂, [block diagram …, which is positive if proposition (i) holds. This proves sufficiency of (i)-(ii) versus controllability of the equivalent feed-forward system. Propositions (iii)-(iv) ensure from Theorem 2.1 that controllability of the equivalent linear system implies controllability of the system (2.3)-(2.4) by allowing the determination of the scalars (𝐴_𝑝,𝑀_𝑝),𝐵_𝑝]] from 𝑢_𝑖(⋅), 𝜈_𝑗𝑖(⋅), in order to drive each subsystem 𝑖∈𝐼 from an arbitrary initial point to a predefined final state.

Remark 2.7. (2.5.1) If condition (iv) is such that the matrices 𝑆_𝑖 and 𝐴⁻¹_𝑖, some 𝑀⁻¹_𝑖, the same extensions referenced in Theorem 2.1(iii) are applicable (namely, see Remark 2.2).
(2.5.2) Theorem 2.6(iv) ensures the existence of an admissible control sequence which leads to the system outside the hyperplane of 𝑖∈𝐼 of insensitivity of the subsystem 𝑏^∗𝑇_{𝑙𝑖+𝑗𝑖−1}(𝑖)𝑥_𝑖(𝑘)≠0 to the 𝑆_𝑖-control [17, 28]. A weak additive signal can be added to the right-hand side of (2.3) at the 𝑢_𝑖(⋅) sample in case of the failure of condition (iv) in Theorem 2.6 at the next step.

3. Optimization Techniques

In this section, a regulation criterion which involves quadratic weighting terms in the state variables is given [1, 39]. The loss function used is 𝑘th where “𝐽_𝑁=12^𝑁_𝑘=0𝑥^𝑇(𝑘)𝑄(𝑘)𝑥(𝑘);𝑄(⋅)=𝑄^𝑇(⋅)≥0;𝑁<∞,(3.1)” denotes positive semidefiniteness. The reason of choosing a finite-time horizon N is then seen from a practical context viewpoint.

In the sequel, the various inverse matrices which take place are assumed to exist. Because of their structures, this hypothesis is not strong.

3.1. Centralized Control

The overall system comprises subsystems, which are interconnected, and the implementation design must deal with the interactions that exist. Centralized control techniques are used for optimization. In the centralized approach, contrarily to the decentralized case [19, 27, 40], a central coordinator exists to take into account the interactions; namely, the coordinator is supposed to be available to the local controller to supply information to each control station from the remaining ones.

The philosophy involved is to define transformed state and input variables, which allows the redefinition of the loss function including quadratic terms of the redefined state and input variables. The idea of using transformations of variables was first pointed out in [1] for the linear continuous case. Its major interest appears in the case of constrained input sequences because, despite this fact, a recursive Riccati expression can be found leading to the optimal feedback solution. Thus, three transformations of variables are made on the equivalent feed-forward linear system of (2.22) (subsequently (2.23) will be invoked) to redefine the state vector equation of standard type (see (B.2) in Appendix B); that is, ≥0 On the other hand, such transformations are applied to the loss function (3.1) to obtain a standard performance index on a given planning horizon 𝑓(𝑘+1)=𝐴⁽³⁾(𝑘)𝑓(𝑘)+𝐶⁽³⁾(𝑘)𝜈⁽³⁾(𝑘);𝑘=0,1,2,….(3.2)

The necessary manipulation on (2.22)–(3.1) to derive (3.2)-(3.3) is detailed in Appendix B.

Now, since (3.3) is a quadratic criterion of standard type (i.e., it involves quadratic terms in both the input and the state variables), a Riccati-matrix sequence may be found, which leads to the optimal solution. The main optimization results are now summarized.

Theorem 3.1. (optimization of the system equivalent inputs). Assume that Theorem 2.6 holds. Then, the optimal equivalent input sequence for the equivalent feed-forward linear system with respect to the loss function (3.1) is given by the following expression: 𝐽_𝑁=12^𝑁_𝑘=0𝑓^𝑇(𝑘)𝑄⁽³⁾(𝑘)𝑓(𝑘)+𝜈^(3)𝑇(𝑘)𝑅⁽³⁾(𝑘)𝜈⁽³⁾(𝑘);𝑁<∞.(3.3) where the optimal equivalent redefined input sequence is 𝜈^∗_𝑗𝑖(𝑘)=𝜈^∗(3)_𝑗𝑖(𝑘+1)−𝐶^𝑇_𝑖(𝑘)𝑄(𝑘+1)𝐶_𝑖(𝑘)⁻¹𝐶^𝑇_𝑖(𝑘)𝑄(𝑘+1)𝐴𝑥(𝑘)−𝑒^𝑇_𝑗𝑖𝑅⁽³⁾⁻¹(𝑘+1)𝑊^(2)𝑇(𝑘+1)𝐴𝑥(𝑘);𝑖∈𝐼;𝑘=0,1,2,…,𝑁,(3.4) with 𝜈^(3)∗_𝑗𝑖(𝑘+1)=−𝑒^𝑇_𝑗𝑖𝑅⁽³⁾(𝑘+1)+𝐶^(3)𝑇(𝑘+1)𝑃(𝑘+2)𝐶⁽³⁾(𝑘+1)⁻¹×𝐶^(3)𝑇(𝑘+1)𝑃(𝑘+2)𝐴⁽³⁾(𝑘+1)𝐴𝑥(𝑘);𝑘=0,1,2,…,𝑁(3.5) being the recursive Riccati matrices; that is, 𝑃(⋅)

Outline of the Proof
It is obvious by direct substitution while dealing with the aforementioned optimization strategies [1, 29, 38] applied to the feed-forward linear equivalent system of (2.22).

Corollary 3.2 (generation of the optimal input sequences). Under the same hypotheses as in Theorem 3.1, consider the sequences of scalars 𝑃(𝑘)=𝑄⁽³⁾(𝑘)+𝐴⁽³⁾(𝑘)𝑃(𝑘+1)×{𝐼−𝐶⁽³⁾(𝑘)𝑅⁽³⁾(𝑘)+𝐶⁽³⁾(𝑘)𝑃(𝑘+1)𝐶⁽³⁾(𝑘)⁻¹𝐶⁽³⁾(𝑘)𝑃(𝑘+1)}𝐴⁽³⁾(𝑘);𝑃(𝑁)=𝑄⁽³⁾(𝑁);𝑘=0,1,…,𝑁.(3.6) with 𝑢_𝑖(𝑘)≡⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)⁻¹𝜈^∗_𝑗𝑖(𝑘);all𝑗∈𝑁_𝑖,𝑖𝑓𝑗_𝑖∈𝑁⁻_𝑖,𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)⁻¹𝜈^∗_𝑗𝑖(𝑘);all𝑗∈𝑁_𝑖,𝑖𝑓𝑗_𝑖∈𝑁⁺_𝑖,0;𝑖𝑓𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)=0,𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)=0∀𝑗∈𝑁_𝑖;𝑘=0,1,2,…,𝑁(3.7) being the 𝜈^∗_𝑗𝑖(𝑘),𝑘=0,1,2,… component of the optimal equivalent inputs obtained from Theorem 3.1.
Then, the following propositions hold.
(a) The optimal input sequence with respect to the loss function (3.1) is 𝑗_𝑖th if the interval 𝑢^∗_𝑖(𝑘)=𝑢_𝑖(𝑘),𝑘=0,1,2,… is allowed as input definition domain, (−∞,∞)⊂ℝ. (b)all𝑖∈𝐼  If|𝑢_𝑖(𝑘)|≤𝑀<∞, then the optimal input sequence with respect to the loss function (3.1) is all𝑖∈𝐼,  𝑢^∗_𝑖=sat[𝑢_𝑖(𝑘)], where the saturation function is defined by all𝑖∈𝐼,𝑘=0,1,2,…;

Proof. It follows from Theorem 3.1 and the fact that, from (3.7), propositions (a)-(b) in Corollary 3.2 imply that sat𝑔(𝑥)≡⎧⎪⎨⎪⎩𝑀sgn𝑔(𝑥);𝑖𝑓|||𝑔(𝑥)|||>𝑀,𝑔(𝑥);𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.(3.8) belongs to the boundary of its definition domain (and vice versa), and the fact that the optimal Hamiltonian function associated with the loss function (3.1) is a strictly convex function of the equivalent input sequence 𝑢_(⋅)(⋅).

Remark 3.3. (3.1.1) Note that the equivalent and the true input sequences are obtained in a centralized way because of the coupling between the different 𝜈_(⋅)(⋅) through the states of the remaining 𝑆_𝑖
(3.1.2) Also, note that (2.7) can be rewritten as𝑆_𝑗,𝑖,𝑗∈𝐼. where 𝑥(𝑘+1)=𝐴𝑥(𝑘)+^𝑝_𝑗=1𝐻_𝑗(𝑘)𝑢_𝑗(𝑘);𝑘=0,1,2,…,(3.9) By applying optimization techniques to (3.9), one could obtain, at first, an optimal 𝐻_𝑖(𝑘)≡𝐵_𝑖𝑥(𝑘);all𝑖∈𝐼;𝑘=0,1,2,….(3.10)-sequence. But this would be nonsense because the Riccati-matrix sequence would be dependent on the future measurements of the state vector.
(3.1.3) Due to the structure of the 𝑢_(⋅)(⋅)-vectors which are state-dependent (see (A.1)) and so unknown in advance, both the optimal equivalent and the true inputs cannot be determined through Riccati-matrix sequence from the results in Theorem 3.1 and Corollary 3.2.
Suboptimal schemes are now presented which are alternative to the optimal one of Theorem 3.1.

3.2. Centralized Suboptimal Schemes

3.2.1. Centralized Suboptimal Modified Scheme 1

It consists of the following variations on the optimal scheme from Theorem 3.1 and Corollary 3.2:

(i)to use finite (sufficient small for the problem coherence) time-sliding optimization horizons. Namely, the time horizon for optimization is one-step advanced when the input is determined at each step. Only the first input associated with each optimization horizon is in fact applied;(ii)to approximate the 𝑐_(⋅)-vectors by their values by using each current value at the first sample of each optimization horizon; that is, 𝐶_(⋅)(⋅), 𝐶(𝑡)=𝐶(𝑘)

3.2.2. Centralized Suboptimal Modified Scheme 2

It is based on the application of (2.23) and (A.2) in Appendix A, by implementing the optimization scheme while neglecting the dependence of all𝑘≤𝑡≤𝑘+𝑁;𝑘=0,1,2,…. on the 𝜉(⋅)-optimal equivalent input sequence. This leads to a suboptimal Riccati-matrix sequence, which does not depend on the state vector, and thus, implementable. The presence of the deterministic disturbance 𝜈⁽³⁾_(⋅)(⋅) obliges to use an additional vector (denoted by 𝜉(⋅)) in the optimization procedure in Theorem 3.1 in order to maintain a Riccati-type solution. The costate associated with the Hamiltonian of the loss function (3.1) verifies the so-called “modified Riccati transformation” [26]; namely, 𝜌(⋅) Thus, by taking into account that now 𝜂(𝑡)(costate)≡𝑃(𝑡)𝑓(𝑡)−𝜌(𝑡). one finds out the following result instead of Theorem 3.1. See Appendix B for particular mathematical details.

Theorem 3.4. (optimization of the equivalent inputs for the system representation including a deterministic disturbance). Under the same assumptions as in Theorem 3.1, it follows that the optimal equivalent input sequence for the equivalent feed-forward linear representation involving a deterministic disturbance with respect to the loss function (3.1) is given by 𝑓(𝑘)≡𝐴𝑥(𝑘−1)+𝜉(𝑘−1), where the optimal equivalent redefined input sequence is 𝜈^∗_𝑗𝑖(𝑘)=𝜈^(3)∗_𝑗𝑖(𝑘+1)−𝑒^𝑇_𝑗𝑖𝑄(𝑘+1)𝑒_𝑗𝑖⁻¹𝑒^𝑇_𝑗𝑖𝑄(𝑘+1)−𝑒^𝑇_𝑗𝑖𝑅⁽³⁾⁻¹(𝑘+1)𝑊^(2)𝑇(𝑘+1)𝐴𝑥(𝑘)+𝜉(𝑘);𝑎𝑙𝑙𝑖∈𝐼;𝑘=0,1,2,…,(3.11) The recursive Riccati-matrix sequence is now given by 𝜈^∗(3)_𝑗𝑖(𝑘)=𝑒^𝑇_𝑗𝑖𝑅⁽³⁾(𝑘)+𝐶^(3)𝑇(𝑘)𝑃(𝑘+1)𝐶⁽³⁾(𝑘)⁻¹𝐶^(3)𝑇(𝑘)×𝜌(𝑘+1)−𝑃(𝑘+1)𝐴⁽³⁾(𝑘)𝐴𝑥(𝑘−1)+𝜉(𝑘−1)+𝜉(𝑘);𝑎𝑙𝑙𝑖∈𝐼;𝑘=1,2,3,….(3.12)

Proof. It follows immediately from applying the three transformations of variables of Appendix B to (2.23) and (A.2), and from the fact that 𝑃(𝑘)=𝑄⁽³⁾(𝑘)+𝐴⁽³⁾(𝑘)×{𝑃(𝑘+1)−𝑃(𝑘+1)𝐶⁽³⁾(𝑘)𝑅⁽³⁾(𝑘)+𝐶^(3)𝑇(𝑘)𝑃(𝑘+1)𝐶⁽³⁾(𝑘)⁻¹𝐶⁽³⁾(𝑘)𝑃(𝑘+1)}𝐴⁽³⁾(𝑘);𝑃(𝑁)=𝑄⁽³⁾(𝑁);𝑘=0,1,2,…,𝑁,(3.13)𝜌(𝑘+1)={𝐼+𝑄⁽³⁾(𝑘)−𝑃(𝑘)𝐴⁽³⁾⁻¹(𝑘)𝑅⁽³⁾⁻¹(𝑘)𝐶⁽³⁾(𝑘)𝐶^(3)𝑇(𝑘)𝐴^(3)−𝑇(𝑘)}−𝐴^(3)−𝑇(𝑘){𝜌(𝑘)+𝑃(𝑘)−𝑄⁽³⁾(𝑘)𝐴⁽³⁾⁻¹(𝑘)𝜉(𝑘)};𝜌(𝑁)=0;𝑘=0,1,2,…,𝑁.(3.14) or equivalently 𝑓(𝑘)≡𝐴𝑥(𝑘−1)+𝜉(𝑘−1), and the fact that the 𝑓(𝑘+1)=𝐴⁽³⁾(𝑘)𝑓(𝑘)+𝐶⁽³⁾(𝑘)𝜈⁽³⁾(𝑘)+𝜉(𝑘);𝑘=0,1,2,…,(3.15) vectors are now 𝐶⁽³⁾(⋅)

Remark 3.5. (3.2.1) Note that Corollary 3.2 also applies here to obtain the optimal system input sequence under constraints (see (3.7) and (3.8)).
(3.2.2) The 𝐶⁽³⁾(𝑘)=𝐶⁽³⁾≡𝐴𝑒_{𝑗1+𝑙1−1},𝐴𝑒_{𝑗2+𝑙2−1},…,𝐴𝑒_{𝑗𝑝+𝑙𝑝−1}with𝑘=0,1,2,3,….(3.16)-vector (related to the deterministic disturbance in the modified Riccati transformation is dependent on the 𝜌(⋅)-sequence. To solve such a dependence, the 𝜈⁽³⁾(⋅)-vector of (3.14) can be decomposed into 𝜌(𝑘+1) being 𝜌(𝑘+1)=𝑚(𝑘)𝜌(𝑘)+𝑀(𝑘)𝜉(𝑘);𝑘=0,1,2,…,(3.17), 𝑀(⋅) matrices defined by𝑀(⋅)∈𝑅^{(𝑛+𝑎)𝑥(𝑛+𝑎)} Thus, the following result yields.

Corollary 3.6 (useful implementability result for obtaining the optimal equivalent input). If the same assumptions from Theorem 3.4 are fulfilled, then the optimal equivalent input for the system interpretation including a deterministic disturbance with respect to the loss function (3.1) can be rewritten as follows (see (3.11)): 𝑀(𝑘)=𝐼+𝑄⁽³⁾(𝑘)−𝑃(𝑘)𝐴⁽³⁾⁻¹(𝑘)𝑅⁽³⁾⁻¹(𝑘)𝐶⁽³⁾(𝑘)𝐶^(3)𝑇(𝑘)𝐴^(3)−𝑇(𝑘)⁻¹;𝑘=0,1,2,3,…,𝑀(𝑘)=𝐼+𝑄⁽³⁾(𝑘)−𝑃(𝑘)𝐴⁽³⁾⁻¹(𝑘)𝑅⁽³⁾⁻¹(𝑘)𝐶⁽³⁾(𝑘)𝐶^(3)𝑇(𝑘)𝐴^(3)−𝑇(𝑘)⁻¹×𝐴^(3)−𝑇(𝑘)𝑃(𝑘)−𝑄⁽³⁾(𝑘)𝐴⁽³⁾⁻¹(𝑘);𝑘=0,1,2,3,….(3.18)where the matrix 𝜈^∗_𝑗𝑖(𝑘)={1+𝑒^𝑇_𝑗𝑖{𝑅⁽³⁾(𝑘+1)+𝐶^(3)𝑇(𝑘+1)𝑃(𝑘+2)𝐶⁽³⁾(𝑘+1)⁻¹𝐶^(3)𝑇(𝑘+1)𝑃(𝑘+2)𝐴⁽³⁾(𝑘+1)+𝑅⁽³⁾⁻¹(𝑘+1)𝑊⁽²⁾(𝑘+1)𝐹(𝑘+1)+𝑒^𝑇_𝑗𝑖𝑄(𝑘+1)𝑒_𝑗𝑖⁻¹}𝑒_𝑗𝑖𝑄(𝑘+1)𝐹(𝑘+1)}⁻¹×{𝑒^𝑇_𝑗𝑖𝑅⁽³⁾(𝑘+1)+𝐶^(3)𝑇(𝑘+1)𝑃(𝑘+2)𝐶⁽³⁾(𝑘+1)⁻¹𝐶^(3)𝑇(𝑘+1)×𝑀(𝑘+1)−𝑃(𝑘+2)𝜉(𝑘+1)+𝑀(𝑘+1)𝜌(𝑘+1)−𝑃(𝑘+2)𝐴⁽³⁾(𝑘+1)𝐴𝑥(𝑘)−𝑒^𝑇_𝑗𝑖𝑄(𝑘+1)𝑒_𝑗𝑖⁻¹𝑒^𝑇_𝑗𝑖𝑄(𝑘+1)𝐴𝑥(𝑘)−𝑒^𝑇_𝑗𝑖𝑅⁽³⁾⁻¹(𝑘+1)𝑊^(2)𝑇(𝑘)𝐴𝑥(𝑘)};all𝑖∈𝐼;𝑘=0,1,2,3,…,(3.19) is partitioned as 𝐹(⋅)∈𝑅^{(𝑛+𝑎)×𝑝} So that each vector 𝐹(𝑘)≡𝐹₁(𝑘),𝐹₂(𝑘),…,𝐹_𝑝(𝑘);𝑘=0,1,2,3,….(3.20)𝐹_𝑖(𝑘)∈𝑅^𝑛+𝑎

Proof (outline). Equation (3.19) can be obtained after direct calculations if (3.17) and (3.18) are substituted into (3.11) while taking Remark 3.5 into account.

3.3. Considerations about Implementation

The following remarks must be pointed out.

(1)The system implementation scheme, which has been reported in the paper, has appeared to be suboptimal twice. First of all, the optimization scheme is dependent on the state/input vectors (some estimates must be made) according to the two alternative interpretations of the feed-forward linear system (2.13). This fact is due to the nature of the decomposition methods which have been applied and also due to the over determination problem one must deal with when the system inputs are generated from the associate equivalent ones. The performance degradation (i.e., the optimality losses being inherent to the applied suboptimization procedures) could be studied through direct calculations. However, the hypotheses that have been taken into account allow the system implementation.(2)In summary, the steps that the designer ought to follow in the implementation environment are as follows:(a)to set up a finite-time sliding optimization horizon, 𝐹_𝑖(𝑘)≡^𝑝_𝑖=1^{𝑛𝑖+𝑎𝑖}_{𝑗=1,𝑗≠𝑖}{𝑏^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖}⁻¹×{𝑏^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖};all𝑖∈𝐼,𝑘=0,1,2,….(3.21), [𝑘,𝑘+𝑁], integers,(b)to test the system controllability according to Theorem 2.6 and/or Corollary3.2,(c)to estimate the state vector (suboptimal modified scheme 1 of Section 3.2.1) and the equivalent deterministic disturbance (suboptimal modified scheme 2 of Section 3.2.2) on the optimization horizon 1≤𝑁<∞ To do that, the estimated system is chosen from (3.2) as [𝑘,𝑘+𝑁]. and the equivalent input estimates become 𝑓̂‌_𝑘+1=𝐴⁽³⁾_1,𝑘+1𝑓̂‌_𝑘+𝐶⁽³⁾_1,𝑘𝜈⁽³⁾_𝑘−1,𝑓̂‌_𝑘+𝑗=𝐴⁽³⁾_1,𝑘+𝑗𝑓̂‌_{𝑘+𝑗−1}+𝐶⁽³⁾_{𝑘+𝑗−1}𝜈⁽³⁾_{𝑘+𝑗−1};𝑗∈2,𝑁,(3.22)(d)to compute the Riccati matrix from (3.6) or (3.13), respectively, and the system matrices from (B.6) and (B.8) of Appendix B, or (B.8) and (3.16), respectively.(e)to generate the suboptimal input sequence 𝜈̂_𝑘⁽³⁾𝑘=𝜈̂_𝑘−1⁽³⁾𝑘−1=𝜈⁽³⁾_𝑘−1,(3.23)𝜈̂_𝑘+𝑗⁽³⁾𝑘=𝜈̂_𝑘+𝑗⁽³⁾𝑘−1;allinteger𝑘≠𝑗,𝑗∈1,𝑁,(3.24)𝜈̂_𝑘+𝑗⁽³⁾𝑘=0;if𝑘=𝑗,(3.25) from (3.4) or (3.11), respectively.(f)to obtain the suboptimal input sequence 𝜈^∗_𝑗𝑖(𝑘),alli∈𝐼,𝑘=0,1,2,… by applying Corollary 3.2. One must have only input controls as inputs to the system; namely, p inputs.

Remark 3.7. (3.3.1) Note that in (3.22)–(3.25), the following system representation is assumed to be 𝑢^∗_𝑖(𝑘),alli∈𝐼,𝑘=0,1,2,3,…,𝑁, instead of (3.2), and the subscript “𝑓_𝑘+1=𝐴⁽³⁾_𝑘𝑓_𝑘+𝐶⁽³⁾_𝑘𝜈⁽³⁾_𝑘;𝑘=0,1,2,…(3.26)” in the matrices is related to the suboptimal modified scheme 1. Similar expressions can be derived for the alternative interpretation scheme (3.15) of Section 3.2.2.
(3.3.2) Also, note that such “fictitious” states must be estimated in order to later reupdate the input sequence. In fact, only the 1 equivalent input is applied in each optimization horizon.
(3.3.3) The usual asymptotic stability tests via Lyapunov's theorem (which is applied to linear time-invariant optimal systems with respect to quadratic criteria) are not applicable here because the system under study is neither invariant nor linear. Besides, standard stability proofs for optimal regulators cannot be applied to the optimal scheme, which has been obtained because of its nonimplementability. The suboptimal modified schemes can lead to instability (or at least stability can not be proved) when 𝜈⁽³⁾(𝑘−1) . But this fact does not affect the system stability in the current context of this paper since a finite-time sliding for the optimization horizon is used.
In particular, these circumstances make the implementation of a decentralized design quite difficult because some links of the equivalent feed-forward system (2.13) are cut off. However, stability can be ensured if 𝑁→∞ is a Hurwitz's matrix and 𝐴 is a bounded sequence, if a saturation type rule is applied when computing the optimal redefined equivalent inputs 𝑢(⋅) in the various schemes.

4. Numerical Simulation

4.1. Example 1

This first example implemented deals with the following ninth-order discrete model: 𝜈⁽³⁾_(⋅)(⋅)𝑆₁∶𝑥₁(𝑘+1)≡𝐴₁(𝑇)+𝑢₁(𝑘)𝐵₁(𝑇)𝑥₁(𝑘)+𝐷₁(𝑇)𝑧₁(𝑘);𝑢₁∈𝑅,𝑦₁(𝑘)≡𝑥₁(𝑘);𝑥₁(⋅)∈𝑅^2𝑥1,𝐻₁∶𝑧₁(𝑘+1)≡𝑀₁(𝑇)+𝑢₁(𝑘)𝑃₁(𝑇)𝑧₁(𝑘)+𝐿₁₁(𝑇)𝑥₁(𝑘);𝑧₁(⋅)∈𝑅^2𝑥1,𝐻₂∶𝑧₂(𝑘+1)≡𝑀₂(𝑇)𝑧₂(𝑘)+𝐿₂₁(𝑇)𝑥₁(𝑘);𝑧₂(⋅)∈𝑅^3𝑥1,𝐻₃∶𝑧₃(𝑘+1)≡𝑀₃(𝑇)𝑧₃(𝑘)+𝐿₃₁(𝑇)𝑥₁(𝑘);𝑧₃(⋅)∈𝑅^2𝑥1,(4.1) Table 1 displays the entries of the constant discrete matrices related to the composite homogeneous bilinear structures all𝑘≥0,𝑘∈𝑧., 𝑆₁, 𝐻_𝑖, in (4.1). The simulations have been performed with the centralized suboptimal modified scheme 1, including free control type (namely, without saturation) only according to Corollary 3.2.

Table 1 

Matrix entries of the simulated system structure.

The steady state of the signal trajectory under no control leads to the following parameters: 𝐿₂₁≡⎡⎢⎣00(−6.9022+44.5212𝑇)𝑇0⎤⎥⎦, 𝑦^∗_1,st=21.8, 𝑦^∗_2,st=238.63,  𝑢^∗_st=400 samples. The system performance is studied as a function of the optimization horizon, 𝑁^∗_st=57 Also, the plant settings are: Initial conditions, 𝑁., 𝑥₁(0)=[20.0,250.0]^𝑇, 𝑧₁(0)=[0.0,0.0]^𝑇, and 𝑧₂(0)=[0.0,0.0,0.0]^𝑇; sampling period 𝑧₃(0)=[0.0,0.0]^𝑇 seconds and working horizon 𝑇=0.03 samples.

The obtained results are shown in Table 2 where the deviations (percent) are related to the normal operation signals (i.e., WH=150, 𝑦^∗_1,st and 𝑦^∗_2,st) and the CCT is the average time the computer needs to perform a computation iteration of the optimization horizon.

Table 2 

Typical results of the discretely controlled plant without saturation as a function of the optimization horizon.

Exhaustive simulations about variations of the optimization horizon have shown that the discrete control improves the system response from a better transient characteristic point of view. However, there exist deviations of the steady-state signal tracking related to the no-control operation mode signals. Also, the experiments become time-consuming and memory-storage expensive as the optimization horizon does increase. For instance, some of the simulated examples have been plotted in Figure 2, consisting of the two components of the plant output and the waveform of the control input to the system (see Table 2 for the definitions of the involved signals related to the optimization horizon sizes).

(a) Relation between the transient characteristics of the first output component and the optimization horizon

(b) Second output component transient characteristics as a function of the optimization horizon

(c) Input control performance as affected by the optimization horizon

Figure 2 

4.2. Example 2

A discrete bilinear system which corresponds to the discretization of a bilinear continuous one for small sampling period and zero-order hold (ZOH) is given as defined in Section 2.2, by − where 𝑥(𝑘+1)=𝐴(𝑇)+𝑢(𝑘)𝑁(𝑇)𝑥(𝑘)+𝐸(𝑇)𝑢(𝑘)+𝑤(𝑘),𝑦(𝑘)=𝐶𝑥(𝑘),(4.2), 𝑥(𝑘)≡[𝑥^𝑇₁(𝑘),𝑥^𝑇₂(𝑘),𝑥^𝑇₃(𝑘)]^𝑇 and there exists a deterministic disturbance 𝑦(𝑘)≡[𝑦^𝑇₁(𝑘),0,0]^𝑇 so that 𝑤(𝑘)≡[𝑤^𝑇₁(𝑘),0,0]^𝑇 being 𝑤₁(𝑘)≡𝑆₁(𝑇)𝑢²(𝑘)+[𝑆₂(𝑇)+𝑢(𝑘)𝑆₃(𝑘)]𝑞(𝑘), the sampling period, 𝑇 denotes transpose and the sampling instant is defined as (⋅)^𝑇

In (4.2), the coefficients of the block matrices are composed of the following structures: 𝑡_𝑘+1=𝑡_𝑘+𝑇. In particular, (4.2) and (4.3) correspond to the linearized discretization of an electromechanical system consisting of two amplidynes acting on the DC-motor shown in Figure 3, with a regulation mechanism incorporated as introduced in Section 3. This regulation scheme is not restricted to considering either the relevant control current or voltage to be constant, which emphasizes the interest of the bilinear modelling.

Figure 3 

Electromechanical system with interconnected subsystems.

Figure 4 plots a typical tracking response output of the system as a function of the sampling period, for a fixed optimization horizon 𝐴(𝑇)=⎡⎢⎢⎢⎣𝐴₁(𝑇)0𝑀₁(𝑇)𝑃₂(𝑇)𝐴₂(𝑇)𝑅₂(𝑇)𝑃₃(𝑇)𝑅₃(𝑇)𝐴₃(𝑇)⎤⎥⎥⎥⎦;𝑁(𝑇)=⎡⎢⎢⎢⎣𝑁₁(𝑇)0𝑆₂(𝑇)000000⎤⎥⎥⎥⎦;𝐸(𝑇)=⎡⎢⎢⎢⎣𝐸₁(𝑇)𝐸₂(𝑇)𝐸₃(𝑇)⎤⎥⎥⎥⎦;𝐶=𝐼00.(4.3) samples and 𝑁=50 being unity in the quadratic regulation performance criterion (see (3.1)).

Figure 4 

Motor angular speed output as function of the sampling period with respect to a step reference signal of 277 rad/seg.

The control action implemented can be used as an alternative to the traditional strategies of linear control on electrical DC-machines. In general, it has been noticed that the discrete control implies better transients related to the linear approach because those transient responses are faster, that is, present lower settling times with greater shooting parameters. Besides, it has also been observed that the discrete control performance improves, in general, as the optimization horizon is increased for different experiments, with respect to that of the nominal plant when no control action is reached.

5. Conclusions

A multivariable invariant discrete-time bilinear system being composed of interconnected subsystems has been studied. An equivalent feed-forward linear system with equivalent inputs, which are derived from products state-input, has been given. Then, the system has been suboptimized with respect to a quadratic finite-time optimization horizon in order to drive each subsystem from any arbitrary initial point to a predefined final state.

The suboptimization has been made by neglecting either the time-dependence of the control vectors on the state vector (modified suboptimal scheme 1) or the dependence of a deterministic disturbance vector on the equivalent input sequence (modified suboptimal scheme 2). These approaches have effect only on the system implementability rather than on its stability. Besides, stability proofs lead to drawbacks when the optimization horizon is infinite because of the suboptimal real implementation. This is also translated into drawbacks when implementing decentralized expected results. The proposed suboptimal schemes have been proven by means of realistic examples.

Appendices

A. Vector Variation for the Two Alternative Interpretation Schemes of the Feed-Forward Linear Systems

Section 2.4 deals with two alternative representations of the feed-forward linear system (2.13). This appendix is devoted to the derivation of the vectors 𝑄 and 𝐶_(⋅)(⋅) involved in this section.

Substituting (2.18) into (2.13) and grouping terms, one obtains

𝜉(⋅) or alternatively 𝑐_{𝑗+𝑙𝑖−1}(𝑘)≡𝑒_{𝑗𝑖+𝑙𝑖−1}+∑^{𝑛𝑖+𝑎𝑖}_{𝑗=1/𝑗≠𝑗𝑖}{{𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+∑^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖}⁻¹×{𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+∑^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖}𝑒_{𝑗+𝑙𝑖−1}},(A.1) where 𝜉(𝑘)≡⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩^𝑝_𝑖=1^{𝑛𝑖+𝑎𝑖}_{𝑗=1/𝑗≠𝑗𝑖}{𝑏^𝑇_𝑗𝑖(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗𝑖−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖}⁻¹×{𝑏^𝑇_𝑗(𝑖)𝑥_𝑖(𝑘)𝛿_𝑗𝑖+[𝑝^𝑇_{𝑗−𝑛𝑖}(𝑖)𝑧_𝑖(𝑘)+^𝑝_𝑙=1𝑚^𝑇_{𝑗−𝑛𝑖}(𝑖,𝑙)𝑦₁(𝑘)]𝛿_𝑗𝑖}𝜈_𝑗𝑖(𝑘)𝑒_{𝑗+𝑙𝑖−1}};ifatleastanadmissible𝑗_𝑖∈𝑁_𝑖existsatthecurrentsamplinginstant,0;otherwise,(A.2)

B. Derivation of the Optimization Equations

In the sequel, in order not to repeat tedious notation, the 𝛿_𝑗𝑖≡1−𝛿_𝑗𝑖,with𝛿_𝑗𝑖=1,if𝑗_𝑖∈𝑁_𝑖,and𝛿_𝑗𝑖=0,if𝑗_𝑖∈𝑁⁺_𝑖,all𝑖∈𝐼,𝑘=0,1,2,….-index (related to the system inputs and associate equations of the feed-forward linear system 2.13)) will be denoted by the subscript 𝑗_𝑖 Also, three modified auxiliary inputs 𝑖,𝑖∈𝐼., 𝜈⁽¹⁾(𝑘) and 𝜈⁽²⁾(𝑘), are calculated from 𝜈⁽³⁾(𝑘) and introduced due to the bilinear terms.

B.1. First Transformation

Let us define new variables as 𝜈(𝑘), with 𝐶⁽¹⁾_𝑗(𝑘)≡𝐴𝐶_𝑗(𝑘−1);𝜈⁽¹⁾_𝑗(𝑘)≡𝜈_𝑗(𝑘−1),𝑓(𝑘)≡𝐴𝑥(𝑘−1)=𝑥(𝑘)−^𝑝_𝑖=1𝐶_𝑖(𝑘−1)𝜈_𝑖(𝑘−1)(B.1) and 𝑗∈𝐼.

Substituting (B.1) into (2.22), one obtains 𝑘≥1 Also, the loss function (3.1) can be equivalently rewritten as 𝑓(𝑘+1)=𝐴𝑓(𝑘)+^𝑝_𝑖=1𝐶⁽¹⁾_𝑖(𝑘)𝜈⁽¹⁾_𝑖(𝑘);𝑘=0,1,2,3,….(B.2) where 𝐽_𝑁=12{^𝑁_𝑘=0𝑓^𝑇(𝑘)𝑄(𝑘)−^𝑝_𝑖=1𝑟⁻¹_𝑖(𝑘)ℎ_𝑖(𝑘)ℎ^𝑇_𝑖(𝑘)𝑓(𝑘)+^𝑝_𝑖=1{𝑟_𝑖(𝑘)𝜈⁽¹⁾_𝑖(𝑘)+𝑟⁻¹_𝑖(𝑘)ℎ^𝑇_𝑖(𝑘)𝑓(𝑘)²+2^𝑝_{𝑗=1/𝑗>𝑖}𝐶^(1)𝑇_𝑖(𝑘)𝐴^−𝑇ℎ_𝑗(𝑘)𝜈⁽¹⁾_𝑖(𝑘)𝜈⁽¹⁾_𝑗(𝑘)}};𝑁<∞,(B.3)𝑟_𝑗(𝑘)≡𝐶^(1)𝑇_𝑗(𝑘)𝐴^−𝑇𝑄(𝑘)𝐴⁻¹𝐶⁽¹⁾_𝑗(𝑘)=𝐶^𝑇_𝑗(𝑘−1)𝑄(𝑘)𝐶_𝑗(𝑘−1),ℎ_𝑗(𝑘)≡𝑄(𝑘)𝐴⁻¹𝐶⁽¹⁾_𝑗(𝑘)=𝑄(𝑘)𝐶_𝑗(𝑘−1),(B.4).

B.2. Second Transformation

The following new variables are defined as follows: all𝑗∈𝐼;𝑘=1,2,3,…𝐶⁽²⁾_𝑗(𝑘)≡𝐶⁽¹⁾_𝑗(𝑘),𝜈⁽²⁾_𝑗(𝑘)≡𝜈⁽¹⁾_𝑗(𝑘)+𝑟⁻¹_𝑗(𝑘)ℎ^𝑇_𝑗(𝑘)𝑓(𝑘)=𝜈_𝑗(𝑘−1)+𝑥^𝑇(𝑘−1)𝐴^𝑇𝑄(𝑘)𝐶_𝑗(𝑘−1)𝐶^𝑇_𝑗(𝑘−1)𝑄(𝑘)𝐶_𝑗(𝑘−1)⁻¹,𝐴⁽²⁾(𝑘)≡𝐴−^𝑝_𝑖=1𝑟⁻¹_𝑖(𝑘)𝐶⁽¹⁾_𝑖(𝑘)ℎ^𝑇_𝑖(𝑘)=𝐴−^𝑝_𝑖=1𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘)𝐶_𝑖(𝑘−1)⁻¹𝐴𝐶_𝑖(𝑘−1)𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘),(B.5).

From (B.5), one has equivalently to (B.3);all𝑗∈𝐼;𝑘=1,2,3,… where 𝑓(𝑘+1)=𝐴⁽²⁾(𝑘)𝑓(𝑘)+𝐶⁽²⁾(𝑘)𝜈⁽²⁾(𝑘);𝑘=0,1,2,3,…,(B.6) Also, taking into account the mentioned (B.5), (B.3) becomes as follows: 𝐶⁽²⁾(⋅)≡[𝐶⁽²⁾₁(⋅),𝐶⁽²⁾₂(⋅),…,𝐶⁽²⁾_𝑝(⋅)]∈𝑅^{(𝑛+𝑎)𝑥(𝑛+𝑎)},𝜈⁽²⁾(⋅)≡[𝜈⁽²⁾₁(⋅),𝜈⁽²⁾₂(⋅),…,𝜈⁽²⁾_𝑝(⋅)]∈𝑅^𝑝.(B.7) with 𝐽_𝑁=12^𝑁_𝑘∑₌₀{𝑓^𝑇(𝑘)𝑄⁽²⁾(𝑘)𝑓(𝑘)+∑^𝑝_𝑖=1𝑟_𝑖(𝑘)𝜈⁽²⁾²_𝑖(𝑘)+2∑^𝑝_{𝑗=1/𝑗>𝑖}𝑝_𝑖𝑗(𝑘)𝜈²_𝑖(𝑘)𝜈²_𝑗(𝑘)−2∑^𝑝∑_𝑖=1^𝑝_{𝑗=1/𝑗>𝑖}𝑓^𝑇(𝑘)𝑝_𝑖𝑗(𝑘)𝑡^𝑇_𝑗(𝑘)𝜈²_𝑖(𝑘)+𝑡^𝑇_𝑖(𝑘)𝜈²_𝑗(𝑘)};𝑁<∞(B.8)𝑝_𝑖𝑗(𝑘)=𝑝_𝑗𝑖(𝑘)≡𝐶^(2)𝑇_𝑖𝐴^−𝑇ℎ_𝑗(𝑘)=𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘)𝐶_𝑗(𝑘−1);𝑗≠𝑖,𝑡_𝑖(𝑘)≡𝑟⁻¹_𝑖(𝑘)ℎ^𝑇_𝑖(𝑘)=𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘)𝐶_𝑖(𝑘−1)⁻¹𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘),(B.9)

Equation (B.8) can be compactly rewritten as forall𝑖,𝑗∈𝐼;𝑘=0,1,2,…. where the matrices involved are 𝐽_𝑛=12^𝑁_𝑘=0𝑓^𝑇(𝑘)𝑄⁽²⁾(𝑘)𝑓(𝑘)+𝜈^(2)𝑇(𝑘)𝑅⁽²⁾(𝑘)𝜈⁽²⁾(𝑘)+2𝑓^𝑇(𝑘)𝑊⁽²⁾(𝑘)𝜈⁽²⁾(𝑘);𝑁<∞,(B.10) with 𝑄⁽²⁾(𝑘)≡𝑄(𝑘)−^𝑝_𝑖=1{𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘)𝐶_𝑖(𝑘−1)⁻¹𝑄(𝑘)𝐶_𝑖(𝑘−1)𝐶^𝑇_𝑖(𝑘−1)𝑄(𝑘)−2^𝑝_{𝑗=1/𝑗>𝑖}{𝐶^𝑇_𝑗(𝑘−1)𝑄(𝑘)𝐶_𝑖(𝑘−1)⁻¹𝐶^𝑇_𝑗(𝑘−1)𝑄(𝑘)𝐶_𝑗(𝑘−1)⁻¹×𝐶_𝑖(𝑘−1)𝑄(𝑘)𝐶_𝑗(𝑘−1)𝑄(𝑘)𝐶_𝑖(𝑘−1)𝐶^𝑇_𝑗(𝑘−1)𝑄(𝑘)}}(B.11)𝑅⁽²⁾(𝑘)=𝑅⁽²⁾_𝑖𝑗(𝑘);𝑅⁽²⁾_𝑖𝑗(𝑘)=𝑅⁽²⁾_𝑗𝑖(𝑘)≡⎧⎪⎨⎪⎩𝑟_𝑖(𝑘),if𝑗=𝑖,𝑝_𝑖𝑗(𝑘),if𝑗≠𝑖,(B.12)𝑊⁽²⁾(𝑘)≡𝑊₁(𝑘)𝑊₂(𝑘)(B.13)𝑊₁(𝑘)≡𝑡^𝑇₁(𝑘),𝑡^𝑇₂(𝑘),…,𝑡^𝑇_𝑝(𝑘),𝑊₂(𝑘)≡𝑊_𝑖𝑗(𝑘);𝑊_𝑖𝑗(𝑘)=𝑊_𝑗𝑖(𝑘)≡⎧⎪⎨⎪⎩0if𝑗=𝑖−𝑝_𝑖𝑗(𝑘),if𝑗≠𝑖(B.14), and 𝑟_𝑖(⋅),𝑝_𝑖𝑗(⋅), in (B.12) to (B.14); 𝑡_𝑖(⋅);𝑗≠𝑖, are given, respectively, in (B.4) and (B.9).

Note that the factorization of matrix 𝑖,𝑗∈𝐼in (B.13) is possible from (B.8).

B.3. Third Transformation

Although it is not necessary, this transformation of variables provides immediately a loss function (redefinition of (3.1)), which includes weighting terms associated with both the state vector and the transformed input. Let us define them as follows: 𝑊⁽²⁾(⋅) Substitution of (B.15) and (B.16) into (B.6) yields directly (3.2). Also, substituting (B.15) into (B.10), one obtains (3.3) if the weighting matrices are defined by 𝐶⁽³⁾(𝑘)≡𝐶⁽²⁾(𝑘),𝜈⁽³⁾(𝑘)≡𝜈⁽²⁾(𝑘)+𝑅⁽²⁾⁻¹(𝑘)𝑊^(2)𝑇(𝑘)𝑓(𝑘),(B.15)𝐴⁽³⁾(𝑘)≡𝐴⁽²⁾(𝑘)−𝐶⁽²⁾(𝑘)𝑅⁽²⁾⁻¹(𝑘)𝑊^(2)𝑇(𝑘)∀𝑘=0,1,2,3,….(B.16) with 𝑄⁽³⁾(𝑘)≡𝑄⁽²⁾(𝑘)−𝑊⁽²⁾(𝑘)𝑅⁽²⁾⁻¹(𝑘)𝑊^(2)𝑇(𝑘),𝑅⁽³⁾(𝑘)≡𝑅⁽²⁾(𝑘)(B.17) and 𝑄⁽³⁾(⋅)≥0

Acknowledgments

The authors are very grateful to the Spanish MEC for supporting this work through the research Projects DPI2006-00714 and DPI2006-01677. They are also grateful to the UPV/EHU and the Basque Government MEC for its partial support through Projects EHU 06/88 and S-PE06UN10, respectively. They would like also to thank the comments and suggestions of the anonymous reviewers who have helped to improve the previous version of this manuscript.