# Suboptimal Regulation of a Class of

Authors: M. de la Sen; Aitor J. Garrido; J. C. Soto; O. Barambones; I. Garrido
Journal: Mathematical Problems in Engineering
 (2008)

---

## 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.


---

## Body

## 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 Section2,
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: (2.1)x(k+1)=[A+u(k)B]x(k)+Cu(k);k=0,1,2,…,x(⋅)∈Rn;u(⋅)∈R.
Equation (2.1) can
equivalently be described, if rank (B,C)=1        and then rank  (B)=1   if  C  and  h  are both nonzero, with B=ChT (C,h being unique  n-vectors within a
scalar factor), by the linear feed-forward system with multiplicative input:
(2.2)x(k+1)=Ax(k)+Cν(k);ν(k)≡{u(k)hTx(k);ifC=0,u(k)[hTx(k)+1];ifC≠0
with  k=0,1,2,…;   and where  x(k) and u(k) are the state and the control at time  k, respectively;
and A and B are real constant coefficient matrices of compatible
dimensions.Note that for
arbitrary ranks ofB or (B,C) in (2.1), the decomposition B=∑i=1peibiT can be used, with  ei and  biT           being, respectively, the  ith  unity
vector (i.e., the ith component is unity and the remaining ones are
zero) and the  ith  row vector of B. 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 (A,C) in (2.2) is controllable.(ii)
If(AT,h)  is
not controllable, then the homogeneous bilinear system is not controllable.(iii)
If(A,C) and (AT,h) are both controllable pairs and if hTC, hTA−1C≠0, then the homogenous bilinear system is
controllable in Rn outside the origin.(iv)
The inhomogeneous bilinear system
is not controllable inRn if
(a)
rank [controllability matrix of(AT,h)] = rank [controllability matrix of (AT,h);1]; (1≡n-vector
with every element being unity),(b)
∑i=1nai=1, with  a(⋅)  being the elements of the last row of the A-matrix in the
canonical state-space representation of system (2.2).(v)
The inhomogeneous bilinear system is controllable inRn, if (A,C) is a controllable pair, hTC, hTA−1C≠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 A−1 does not
exist. Such a paper studies the controllability conditions for three cases of
violation of the conditions given in [37] (namely, (A,C)  and      (AT,h)  are controllable pairs and  k0=g.c.d{i/hTAi−1C≠0,0<i≤n2}=1) (g.c.d  denotes the greatest common divisor). Such cases
include the  following:
(a)
rank⌈(A,C)=n,       rank   Ω(hT,A)<n;(b)
rank⌈(A,C)<n,  rank   Ω(hT,A)=n; and(c)
rank⌈(A;C)<n,  rank   Ω(hT,A)<n, with  ⌈(A,C)≡[C,AC,…,An−1C] and  Ω(hT,A)≡[h,ATh,…,A(n−1)Th]T.


### 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  (Si)  of homogeneous type is
considered
(2.3)Si:xi(k+1)=[Ai+ui(k)Bi]xi(k)+Dizi(k),(2.4)yi(k)=Cixi(k);i∈I;I={1,2,…,p};k=0,1,2,…
with the interconnection bilinear
subsystems (Hi), which include linear and bilinear coupling terms,
given by
(2.5)Hi:zi(k+1)=[Mi+ui(k)Pi]zi(k)+∑j=1p[Lij+ui(k)Mij]yj(k);i∈I;k=0,1,2,…,
where  xi(k)∈Rni,zi(k)∈Rai,yi(k)∈Rri,ui(k)∈R  are, respectively, the state vector of the  S̲i-subsystems,
the state vectors of the interaction Hi-subsystems,
the corresponding output and input vectors at time  k, 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(2.6)x¯(k)≡[x¯1T(k),x¯2T(k),…,x¯pT(k)]T∈Rn+a,wherex¯i(k)≡[xiT(k),ziT(k)]T∈Rni+ai,
all  i∈I  with n≡∑j=1pnj,a≡∑j=1paj, and
the extended output vector is defined as  y¯(k)≡[y1T(k),y2T(k),…,ypT(k)]T∈Rr,
with r≡∑j=1prj.
Thus, the following multivariable structure, equivalent to (2.3)–(2.5), is deduced:
(2.7)x¯(k+1)=A¯x¯(k)+∑i=1pui(k)Bi¯x¯(k);k=0,1,2,3,…,(2.8)y¯(k)=C¯x¯(k);k=0,1,2,3,…,
where the matrices of parameters
have the following structures:
(2.9)A¯=[A¯1T,A¯2T,…,A¯pT]T;Ai¯=[A¯i1,A¯i2,…,A¯ii,…,A¯ip]A¯ii=[AiDiLiiCiMi]A¯ij=[00LijCj0](forj≠i)B¯i=[0i−1︸,…,0,B˜iT,0p−1︸,…,0]T;B˜i=[B˜i1,B˜i2,…,B˜ip]B˜ii=[Bi0MiiCiPi]B˜ij=[00MijCj0](forj≠i)C¯=[C¯1T,C¯2T,…,C¯pT];C¯i=[0i−1︸,…,0,Ci,0,0p−1︸,…,0] for all j,i∈I.Remark 2.3.
In general,
the number of subsystemsS(⋅) and the number of interaction
subsystems H(⋅)  are not equal. Thus, assume that  p1
and
p2(p1≠p2) are, respectively, the number of subsystems St,t∈I1={1,2,3,…,p1}, and the number of subsystems  Hs,s∈I2={1,2,…,p2}.      Then, an extended state vector    ξ(⋅)  can be obtained similarly,
if  p=max(p1p2)  and     I=I1UI2.
Taking
into account the dimensions and the structural zeros of theB¯i-matrices, they can  be decomposed into a sum of
dyads as
(2.10)B¯i=∑j=1ni+aiej+li−1b¯j+li−1T(i);for alli∈I,
where ej is the
unity jth       vector, b¯jT(i) is the jth  vector of the B¯i-matrix,    i∈I, and the
lower subscripts 1j,allj∈I, are related to the first nonstructural
zero rows of the B¯i-matrix. If they have additional (nonstructural)
rows being zero,  the corresponding  e(⋅)-vectors can be neglected in (2.10).
From (2.7), it is clear that the rows of the Bi-matrices have the
following structures:
(2.11)b¯jT(i)≡{[0,…,0︷(i−1)block,ej−li+1TBi,0,…,0︷(p−i)blocks];if1i≤j<1i+ni,[ej−li−ni+1TMi1C1,0,…,ej−li.ni+1TMi(i−1)Ci−1,0︷(i−1)blocks,ej−li−ni+1TPi,ej−li−ni+1TMi(i+1)Ci+1,0,…,ej−li−ni+1TMipCi,0︷(p−1)blocks];ifli+nii≤j<1i+ni+ai,for alli∈I. Then one deduces from (2.11) that
(2.12)b¯j+1i−1T(i)x¯≡{bjT(i)xi(k);if1≤j≤ni,pj−niT(i)zi(k)+∑l=1pmj−ni+1Tyl(k);ifni+1≤j≤ni+ai,0;ifni+ai+1≤j≤n+a,withk=0,1,2,3,…,    when bjT(i), pjT(i),    mjT(i,1)  denote, respectively, the jth vector row of the matrices  Bi,Pi,andMi,forallI,1∈I.   From (2.7)
through (2.12) one obtains
(2.13)x¯(k+1)=A¯x¯(k)+∑i=1p∑j=1ni+aeej+li−1νji(k);k=0,1,2,…,
where the
sequences of scalars νji(k) are defined by
(2.14)νji(k)≡{ui(k)bjT(i)xi(k);if1≤j≤ni,ui(k)[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)];ifni+1≤j≤ni+ai,foralli∈I,k=0,1,2,…, are considered as equivalent
(primary) inputs for the equivalent linear systems (2.13) to be determined prior
to the generation of the u(⋅)(⋅)-sequence in the
design implementation.


### 2.3. Obtaining theui(k)-Sequences from the γji(k)-Sequences

Equations
(2.13) and (2.14) constitute a feed-forward linear equivalent system with
multiplicative inputs. But the  pui(k)-inputs cannot, in
general, be determined from the (n+a)>p equivalent νji(k)-inputs. Then, it is
necessary to take a set of p variables between the last ones, which will
be optimized independently. Let us call νji(k),alli∈I.  Thus, let us define the sets(2.15)Ni=Ni−UNi+,Ni−={1,2,…,ni},Ni+={ni+1,ni+2,…,ni+ai},alli∈I, such that for the  kth
current sampling instant,
(2.16)bjiT(i)xi(k)≠0,ifxi(k)≠0,ifji∈Ni−,pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)yi(k)≠0,ifji∈Ni+,
and k=0,1,2,….Then from (2.14), one has the
auxiliary inputs
(2.17)νji(k)≡{[bjiT(i)xi(k)]−ibjT(i)xi(k)νji(k);j∈Ni−;i∈I,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]−1×[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)νji(k)],j∈Ni+;i∈I;k=0,1,2,3,…,
with the real
inputs for the nonsaturated case being(2.18)ui(k)≡{[bjiT(i)xi(k)]−1νji(k);ifji∈Ni,if there exists someji∈Ni−such thatbjiT(i)xi(k)≠0,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)yi(k)]−1νji(k);ifji∈Ni,if there exists someji∈Ni+such thatpji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)≠0,0;ifbjiT(i)xi(k)=0for allji∈Ni−and alsopji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)=0for allji∈Ni+,alli∈I;k=0,1,2,3,….Remark 2.4.
Note
that the minimum energy control (i.e.,ui(⋅)=0,   some i∈I) 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(n+a) equivalent inputs, deals with an over determined problem
when the p<n+a 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 p<s≤n+a equivalent nonzero inputs. Also, note that the decomposition of the B¯i-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,
(2.19)ui(k)≡{ν¯iT(k)B¯ix¯(k)x¯T(k)B¯iTB¯ix¯(k);ifB¯ix¯(k)≠0,0,otherwise,alli∈I;k=0,1,2,3,…,
where
(2.20)ν¯i(k)≡[0,…,0,︷∑j=1i−1(nj+aj)zerosνli(k),…,νni,i(k),νni+1,i(k),…,νni+ai,i(k),0,…,0︷∑j=i+1p(nj+aj)zeros]T∈Rn+a,alli∈I;k=0,1,2,3,….
Note that in
general, the (k) vector has only nonzero components. Equation (2.19) can be
rewritten as
(2.21)ui(k)≡[∑j=1ni+ai{bjT(i)xi(k)δj(i)+[pj−niT0(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯j(i)}2]−1×[∑j=1ni+aiνji(k){bjT(i)xi(k)δj(i)+[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯j(i)}],whereδ¯(i)≡1−δj(i),withδj(i)=1if1≤j≤nand0,ifni+1≤j≤ni+aialli∈I;k=0,1,2,….


### 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):
(2.22)x¯(k+1)=A¯x¯(k)+∑i=1pcj+li−1(k)νji(k)(2.23)=A¯x¯(k)+∑i=1pej+li−1νji(k)+ξ(k),
where the (n+a)-vectors  c(⋅)(⋅) 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 theu(⋅)(⋅)-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 [k′,k′′]∩R+,   k′′=k′+max1≤i≤p(ni+ai),
if the following set of conditions holds.
(i)
The controllability
grammian of the overall system is positive definite in[k′,k′′].(ii)
max[rank(∑l=1pMs1C1)]≥max[rank(LijCj)] with   1≤s≤p,1≤i≤pand1≤j≤p.(iii)
All the pairs  of
block diagram,[(AiT,MiT),B˜i],i∈I,     are nonsingular.(iv)
(bji*T(i)Ci*(k)≠0, bji*T(i)   diagram (Ai−1,Mi−1)Ci*(k)≠0,  alli∈I, allk′≤k≤k′′,   where  C(⋅)*andb(⋅)*  are, respectively, the subvectors
of the matrices  C(⋅)andB(⋅) in (2.3)-(2.4) associated with each subsystem  S(⋅).Sketch of the Proof
Since
the pair[A+B(⋅)G(⋅),B(⋅)]    is controllable for all arbitrary matrix    G(⋅) if the pair [A,B(⋅)]   is controllable [29, 38], then the
controllability of the pair        [Ad,B¯(⋅)] implies the
controllability of the pair [Ad+A°,B¯(⋅)], if there exists a G(⋅)-matrix
such that A°=B¯(⋅)G(⋅)            with      A°=A¯−Ad,    Ad=   block diagram ((A1,M1),(A2,M2),…,(Ap,Mp)). For the equality above relating A°,B¯(⋅),        and               G(⋅)       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 [Ad,B¯(⋅)]
on   [k′,k′′]≡   block diagram controllability grammian [[block
diagram    (A1,M1),B˜1],[block diagram (A2,M2),B˜2], …, [block diagram (Ap,Mp),B˜p]], 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 ui(⋅) from  νji(⋅),   i∈I,   in order to drive each subsystem      Si  from an arbitrary initial point to a predefined final state.Remark 2.7.
(2.5.1) If condition
(iv) is such that the matricesAi−1 and Mi−1,
some i∈I, the same extensions referenced in
Theorem 2.1(iii) are applicable (namely, see Remark 2.2).
(2.5.2) Theorem2.6(iv)
ensures the existence of an admissible control sequence which leads to the system outside the hyperplane of   bli+ji−1*T(i)xi(k)≠0            of insensitivity of the subsystem     Si  to the ui(⋅)-control
[17, 28]. A weak
additive signal can be added to the right-hand side of (2.3) at the kth sample in case of the failure of condition (iv) in Theorem 2.6 at
the next step.


## 2.1. Controllability of a Class of Bilinear Systems

In [17],  the following
typical class of bilinear systems was considered: (2.1)x(k+1)=[A+u(k)B]x(k)+Cu(k);k=0,1,2,…,x(⋅)∈Rn;u(⋅)∈R.
Equation (2.1) can
equivalently be described, if rank (B,C)=1        and then rank  (B)=1   if  C  and  h  are both nonzero, with B=ChT (C,h being unique  n-vectors within a
scalar factor), by the linear feed-forward system with multiplicative input:
(2.2)x(k+1)=Ax(k)+Cν(k);ν(k)≡{u(k)hTx(k);ifC=0,u(k)[hTx(k)+1];ifC≠0
with  k=0,1,2,…;   and where  x(k) and u(k) are the state and the control at time  k, respectively;
and A and B are real constant coefficient matrices of compatible
dimensions.Note that for
arbitrary ranks ofB or (B,C) in (2.1), the decomposition B=∑i=1peibiT can be used, with  ei and  biT           being, respectively, the  ith  unity
vector (i.e., the ith component is unity and the remaining ones are
zero) and the  ith  row vector of B. 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 (A,C) in (2.2) is controllable.(ii)
If(AT,h)  is
not controllable, then the homogeneous bilinear system is not controllable.(iii)
If(A,C) and (AT,h) are both controllable pairs and if hTC, hTA−1C≠0, then the homogenous bilinear system is
controllable in Rn outside the origin.(iv)
The inhomogeneous bilinear system
is not controllable inRn if
(a)
rank [controllability matrix of(AT,h)] = rank [controllability matrix of (AT,h);1]; (1≡n-vector
with every element being unity),(b)
∑i=1nai=1, with  a(⋅)  being the elements of the last row of the A-matrix in the
canonical state-space representation of system (2.2).(v)
The inhomogeneous bilinear system is controllable inRn, if (A,C) is a controllable pair, hTC, hTA−1C≠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 A−1 does not
exist. Such a paper studies the controllability conditions for three cases of
violation of the conditions given in [37] (namely, (A,C)  and      (AT,h)  are controllable pairs and  k0=g.c.d{i/hTAi−1C≠0,0<i≤n2}=1) (g.c.d  denotes the greatest common divisor). Such cases
include the  following:
(a)
rank⌈(A,C)=n,       rank   Ω(hT,A)<n;(b)
rank⌈(A,C)<n,  rank   Ω(hT,A)=n; and(c)
rank⌈(A;C)<n,  rank   Ω(hT,A)<n, with  ⌈(A,C)≡[C,AC,…,An−1C] and  Ω(hT,A)≡[h,ATh,…,A(n−1)Th]T.


## 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  (Si)  of homogeneous type is
considered
(2.3)Si:xi(k+1)=[Ai+ui(k)Bi]xi(k)+Dizi(k),(2.4)yi(k)=Cixi(k);i∈I;I={1,2,…,p};k=0,1,2,…
with the interconnection bilinear
subsystems (Hi), which include linear and bilinear coupling terms,
given by
(2.5)Hi:zi(k+1)=[Mi+ui(k)Pi]zi(k)+∑j=1p[Lij+ui(k)Mij]yj(k);i∈I;k=0,1,2,…,
where  xi(k)∈Rni,zi(k)∈Rai,yi(k)∈Rri,ui(k)∈R  are, respectively, the state vector of the  S̲i-subsystems,
the state vectors of the interaction Hi-subsystems,
the corresponding output and input vectors at time  k, 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(2.6)x¯(k)≡[x¯1T(k),x¯2T(k),…,x¯pT(k)]T∈Rn+a,wherex¯i(k)≡[xiT(k),ziT(k)]T∈Rni+ai,
all  i∈I  with n≡∑j=1pnj,a≡∑j=1paj, and
the extended output vector is defined as  y¯(k)≡[y1T(k),y2T(k),…,ypT(k)]T∈Rr,
with r≡∑j=1prj.
Thus, the following multivariable structure, equivalent to (2.3)–(2.5), is deduced:
(2.7)x¯(k+1)=A¯x¯(k)+∑i=1pui(k)Bi¯x¯(k);k=0,1,2,3,…,(2.8)y¯(k)=C¯x¯(k);k=0,1,2,3,…,
where the matrices of parameters
have the following structures:
(2.9)A¯=[A¯1T,A¯2T,…,A¯pT]T;Ai¯=[A¯i1,A¯i2,…,A¯ii,…,A¯ip]A¯ii=[AiDiLiiCiMi]A¯ij=[00LijCj0](forj≠i)B¯i=[0i−1︸,…,0,B˜iT,0p−1︸,…,0]T;B˜i=[B˜i1,B˜i2,…,B˜ip]B˜ii=[Bi0MiiCiPi]B˜ij=[00MijCj0](forj≠i)C¯=[C¯1T,C¯2T,…,C¯pT];C¯i=[0i−1︸,…,0,Ci,0,0p−1︸,…,0] for all j,i∈I.Remark 2.3.
In general,
the number of subsystemsS(⋅) and the number of interaction
subsystems H(⋅)  are not equal. Thus, assume that  p1
and
p2(p1≠p2) are, respectively, the number of subsystems St,t∈I1={1,2,3,…,p1}, and the number of subsystems  Hs,s∈I2={1,2,…,p2}.      Then, an extended state vector    ξ(⋅)  can be obtained similarly,
if  p=max(p1p2)  and     I=I1UI2.
Taking
into account the dimensions and the structural zeros of theB¯i-matrices, they can  be decomposed into a sum of
dyads as
(2.10)B¯i=∑j=1ni+aiej+li−1b¯j+li−1T(i);for alli∈I,
where ej is the
unity jth       vector, b¯jT(i) is the jth  vector of the B¯i-matrix,    i∈I, and the
lower subscripts 1j,allj∈I, are related to the first nonstructural
zero rows of the B¯i-matrix. If they have additional (nonstructural)
rows being zero,  the corresponding  e(⋅)-vectors can be neglected in (2.10).
From (2.7), it is clear that the rows of the Bi-matrices have the
following structures:
(2.11)b¯jT(i)≡{[0,…,0︷(i−1)block,ej−li+1TBi,0,…,0︷(p−i)blocks];if1i≤j<1i+ni,[ej−li−ni+1TMi1C1,0,…,ej−li.ni+1TMi(i−1)Ci−1,0︷(i−1)blocks,ej−li−ni+1TPi,ej−li−ni+1TMi(i+1)Ci+1,0,…,ej−li−ni+1TMipCi,0︷(p−1)blocks];ifli+nii≤j<1i+ni+ai,for alli∈I. Then one deduces from (2.11) that
(2.12)b¯j+1i−1T(i)x¯≡{bjT(i)xi(k);if1≤j≤ni,pj−niT(i)zi(k)+∑l=1pmj−ni+1Tyl(k);ifni+1≤j≤ni+ai,0;ifni+ai+1≤j≤n+a,withk=0,1,2,3,…,    when bjT(i), pjT(i),    mjT(i,1)  denote, respectively, the jth vector row of the matrices  Bi,Pi,andMi,forallI,1∈I.   From (2.7)
through (2.12) one obtains
(2.13)x¯(k+1)=A¯x¯(k)+∑i=1p∑j=1ni+aeej+li−1νji(k);k=0,1,2,…,
where the
sequences of scalars νji(k) are defined by
(2.14)νji(k)≡{ui(k)bjT(i)xi(k);if1≤j≤ni,ui(k)[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)];ifni+1≤j≤ni+ai,foralli∈I,k=0,1,2,…, are considered as equivalent
(primary) inputs for the equivalent linear systems (2.13) to be determined prior
to the generation of the u(⋅)(⋅)-sequence in the
design implementation.


## 2.3. Obtaining theui(k)-Sequences from the γji(k)-Sequences

Equations
(2.13) and (2.14) constitute a feed-forward linear equivalent system with
multiplicative inputs. But the  pui(k)-inputs cannot, in
general, be determined from the (n+a)>p equivalent νji(k)-inputs. Then, it is
necessary to take a set of p variables between the last ones, which will
be optimized independently. Let us call νji(k),alli∈I.  Thus, let us define the sets(2.15)Ni=Ni−UNi+,Ni−={1,2,…,ni},Ni+={ni+1,ni+2,…,ni+ai},alli∈I, such that for the  kth
current sampling instant,
(2.16)bjiT(i)xi(k)≠0,ifxi(k)≠0,ifji∈Ni−,pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)yi(k)≠0,ifji∈Ni+,
and k=0,1,2,….Then from (2.14), one has the
auxiliary inputs
(2.17)νji(k)≡{[bjiT(i)xi(k)]−ibjT(i)xi(k)νji(k);j∈Ni−;i∈I,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]−1×[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)νji(k)],j∈Ni+;i∈I;k=0,1,2,3,…,
with the real
inputs for the nonsaturated case being(2.18)ui(k)≡{[bjiT(i)xi(k)]−1νji(k);ifji∈Ni,if there exists someji∈Ni−such thatbjiT(i)xi(k)≠0,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)yi(k)]−1νji(k);ifji∈Ni,if there exists someji∈Ni+such thatpji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)≠0,0;ifbjiT(i)xi(k)=0for allji∈Ni−and alsopji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)=0for allji∈Ni+,alli∈I;k=0,1,2,3,….Remark 2.4.
Note
that the minimum energy control (i.e.,ui(⋅)=0,   some i∈I) 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(n+a) equivalent inputs, deals with an over determined problem
when the p<n+a 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 p<s≤n+a equivalent nonzero inputs. Also, note that the decomposition of the B¯i-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,
(2.19)ui(k)≡{ν¯iT(k)B¯ix¯(k)x¯T(k)B¯iTB¯ix¯(k);ifB¯ix¯(k)≠0,0,otherwise,alli∈I;k=0,1,2,3,…,
where
(2.20)ν¯i(k)≡[0,…,0,︷∑j=1i−1(nj+aj)zerosνli(k),…,νni,i(k),νni+1,i(k),…,νni+ai,i(k),0,…,0︷∑j=i+1p(nj+aj)zeros]T∈Rn+a,alli∈I;k=0,1,2,3,….
Note that in
general, the (k) vector has only nonzero components. Equation (2.19) can be
rewritten as
(2.21)ui(k)≡[∑j=1ni+ai{bjT(i)xi(k)δj(i)+[pj−niT0(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯j(i)}2]−1×[∑j=1ni+aiνji(k){bjT(i)xi(k)δj(i)+[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯j(i)}],whereδ¯(i)≡1−δj(i),withδj(i)=1if1≤j≤nand0,ifni+1≤j≤ni+aialli∈I;k=0,1,2,….


## 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):
(2.22)x¯(k+1)=A¯x¯(k)+∑i=1pcj+li−1(k)νji(k)(2.23)=A¯x¯(k)+∑i=1pej+li−1νji(k)+ξ(k),
where the (n+a)-vectors  c(⋅)(⋅) 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 theu(⋅)(⋅)-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 [k′,k′′]∩R+,   k′′=k′+max1≤i≤p(ni+ai),
if the following set of conditions holds.
(i)
The controllability
grammian of the overall system is positive definite in[k′,k′′].(ii)
max[rank(∑l=1pMs1C1)]≥max[rank(LijCj)] with   1≤s≤p,1≤i≤pand1≤j≤p.(iii)
All the pairs  of
block diagram,[(AiT,MiT),B˜i],i∈I,     are nonsingular.(iv)
(bji*T(i)Ci*(k)≠0, bji*T(i)   diagram (Ai−1,Mi−1)Ci*(k)≠0,  alli∈I, allk′≤k≤k′′,   where  C(⋅)*andb(⋅)*  are, respectively, the subvectors
of the matrices  C(⋅)andB(⋅) in (2.3)-(2.4) associated with each subsystem  S(⋅).Sketch of the Proof
Since
the pair[A+B(⋅)G(⋅),B(⋅)]    is controllable for all arbitrary matrix    G(⋅) if the pair [A,B(⋅)]   is controllable [29, 38], then the
controllability of the pair        [Ad,B¯(⋅)] implies the
controllability of the pair [Ad+A°,B¯(⋅)], if there exists a G(⋅)-matrix
such that A°=B¯(⋅)G(⋅)            with      A°=A¯−Ad,    Ad=   block diagram ((A1,M1),(A2,M2),…,(Ap,Mp)). For the equality above relating A°,B¯(⋅),        and               G(⋅)       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 [Ad,B¯(⋅)]
on   [k′,k′′]≡   block diagram controllability grammian [[block
diagram    (A1,M1),B˜1],[block diagram (A2,M2),B˜2], …, [block diagram (Ap,Mp),B˜p]], 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 ui(⋅) from  νji(⋅),   i∈I,   in order to drive each subsystem      Si  from an arbitrary initial point to a predefined final state.Remark 2.7.
(2.5.1) If condition
(iv) is such that the matricesAi−1 and Mi−1,
some i∈I, the same extensions referenced in
Theorem 2.1(iii) are applicable (namely, see Remark 2.2).
(2.5.2) Theorem2.6(iv)
ensures the existence of an admissible control sequence which leads to the system outside the hyperplane of   bli+ji−1*T(i)xi(k)≠0            of insensitivity of the subsystem     Si  to the ui(⋅)-control
[17, 28]. A weak
additive signal can be added to the right-hand side of (2.3) at the kth 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
(3.1)JN=12∑k=0NxT(k)Q¯(k)x¯(k);Q¯(⋅)=Q¯T(⋅)≥0;N<∞,
where    “≥0”     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,
(3.2)f(k+1)=A¯(3)(k)f(k)+C(3)(k)ν(3)(k);k=0,1,2,….
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
(3.3)JN=12∑k=0N[fT(k)Q¯(3)(k)f(k)+ν(3)T(k)R¯(3)(k)ν(3)(k)];N<∞.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 Theorem2.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:
(3.4)νji*(k)=νji*(3)(k+1)−[CiT(k)Q¯(k+1)Ci(k)]−1CiT(k)Q¯(k+1)A¯x¯(k)−ejiTR¯(3)−1(k+1)W¯(2)T(k+1)A¯x¯(k);i∈I;k=0,1,2,…,N,
where the optimal equivalent redefined
input sequence is
(3.5)νji(3)*(k+1)=−ejiT[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1×C(3)T(k+1)P(k+2)A¯(3)(k+1)A¯x¯(k);k=0,1,2,…,N
with  P(⋅)  being the recursive Riccati
matrices; that is, (3.6)P(k)=Q¯(3)(k)+A¯(3)(k)P(k+1)×{I−C(3)(k)[R¯(3)(k)+C(3)(k)P(k+1)C(3)(k)]−1C(3)(k)P(k+1)}A¯(3)(k);P(N)=Q¯(3)(N);k=0,1,…,N.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 Theorem3.1, consider the sequences of scalars(3.7)u¯i(k)≡{[bjiT(i)xi(k)]−1νji*(k);allj∈Ni,ifji∈Ni−,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]−1νji*(k);allj∈Ni,ifji∈Ni+,0;ifbjiT(i)xi(k)=0,pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)=0∀j∈Ni;k=0,1,2,…,N
with νji*(k),k=0,1,2,…   being the jith       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 ui*(k)=ui(k),k=0,1,2,…  if
the interval  (−∞,∞)⊂ℝ is allowed as input definition domain, alli∈I.(b)
If|ui(k)|≤M<∞,  alli∈I, then the optimal input sequence with
respect to the loss function (3.1) is ui*=sat[u¯i(k)],  alli∈I,k=0,1,2,…; where the saturation function is defined by
(3.8)sat[g(x)]≡{Msgn[g(x)];if|g(x)|>M,g(x);otherwise.Proof.
It follows from Theorem3.1 and the fact
that, from (3.7), propositions (a)-(b) in Corollary 3.2 imply that  u(⋅)(⋅)  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 differentSi through the states of the remaining Sj,i,j∈I.
(3.1.2) Also, note that (2.7) can be rewritten as(3.9)x¯(k+1)=A¯x¯(k)+∑j=1pHj(k)uj(k);k=0,1,2,…,
where
(3.10)Hi(k)≡B¯ix¯(k);alli∈I;k=0,1,2,….
By applying optimization techniques to
(3.9), one could obtain, at first, an optimal u(⋅)(⋅)-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 thec(⋅)-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 Theorem3.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 Theorem3.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 theC(⋅)(⋅)-vectors
by their values by using each current value at the first sample of each
optimization horizon; that is, C(t)=C(k),  allk≤t≤k+N;k=0,1,2,….


#### 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   ξ(⋅)   on the  ν(⋅)(3)(⋅)-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, η(t)(costate)≡P(t)f(t)−ρ(t).
Thus, by taking into account that now f(k)≡A¯x¯(k−1)+ξ(k−1), 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 Theorem3.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(3.11)νji*(k)=νji(3)*(k+1)−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)−ejiTR¯(3)−1(k+1)W¯(2)T(k+1)[A¯x¯(k)+ξ(k)];alli∈I;k=0,1,2,…, where the optimal equivalent redefined input
sequence is
(3.12)νji*(3)(k)=ejiT[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)T(k)×[ρ(k+1)−P(k+1){A(3)(k)[A¯x¯(k−1)+ξ(k−1)]+ξ(k)}];alli∈I;k=1,2,3,….
The recursive Riccati-matrix sequence is now given by(3.13)P(k)=Q¯(3)(k)+A¯(3)(k)×{P(k+1)−P(k+1)C(3)(k)[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)(k)P(k+1)}A¯(3)(k);P(N)=Q¯(3)(N);k=0,1,2,…,N,(3.14)ρ(k+1)={I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C¯(3)(k)C¯(3)T(k)A¯(3)−T(k)}−A¯(3)−T(k){ρ(k)+[P(k)−Q¯(3)(k)]A¯(3)−1(k)ξ(k)};ρ(N)=0;k=0,1,2,…,N.Proof.
It follows immediately from applying the three transformations of variables of AppendixB
to (2.23) and (A.2), and from the fact that f(k)≡A¯x¯(k−1)+ξ(k−1),
or equivalently
(3.15)f(k+1)=A¯(3)(k)f(k)+C(3)(k)ν(3)(k)+ξ(k);k=0,1,2,…,
and the fact that the  C(3)(⋅)
vectors are now
(3.16)C(3)(k)=C(3)≡[A¯ej1+l1−1,A¯ej2+l2−1,…,A¯ejp+lp−1]withk=0,1,2,3,….Remark 3.5.
(3.2.1) Note that Corollary3.2
also applies here to obtain the optimal system input sequence under constraints
(see (3.7) and (3.8)).
(3.2.2) Theρ(⋅)-vector
(related to the deterministic disturbance in the modified Riccati
transformation is dependent on the   ν(3)(⋅)-sequence. To
solve such a dependence, the ρ(k+1)-vector
of (3.14) can be decomposed into (3.17)ρ(k+1)=m(k)ρ(k)+M¯(k)ξ(k);k=0,1,2,…, being M(⋅), M¯(⋅)∈R(n+a)x(n+a) matrices defined by(3.18)M(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1;k=0,1,2,3,…,M¯(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1×A¯(3)−T(k)[P(k)−Q¯(3)(k)]A¯(3)−1(k);k=0,1,2,3,….
Thus, the following result yields.Corollary 3.6 (useful implementability result for obtaining the optimal equivalent input).
If
the same assumptions from Theorem3.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.19)νji*(k)={1+ejiT{[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)P(k+2)A¯(3)(k+1)+R¯(3)−1(k+1)W¯(2)(k+1)F(k+1)+[ejiTQ¯(k+1)eji]−1}ejiQ¯(k+1)F(k+1)}−1×{ejiT[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)×{[M¯(k+1)−P(k+2)]ξ(k+1)+M(k+1)ρ(k+1)−P(k+2)A¯(3)(k+1)A¯x¯(k)}−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)A¯x¯(k)−ejiTR¯(3)−1(k+1)W¯(2)T(k)A¯x¯(k)};alli∈I;k=0,1,2,3,…,where the matrix F(⋅)∈R(n+a)×p
is partitioned as (3.20)F(k)≡[F1(k),F2(k),…,Fp(k)];k=0,1,2,3,…. So that each vector  Fi(k)∈Rn+a(3.21)Fi(k)≡∑i=1p∑j=1,j≠ini+ai{bji−niT(i)xi(k)δji+[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]δ¯ji}−1×{bj−niT(i)xi(k)δji+[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯ji};alli∈I,k=0,1,2,….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,[k,k+N],  1≤N<∞,   integers,(b)
to test the system controllability according
to Theorem2.6 and/or Corollary3.2,(c)
to estimate the state vector (suboptimal
modified scheme 1 of  Section3.2.1)
and the equivalent deterministic disturbance (suboptimal modified scheme 2 of Section
3.2.2) on the optimization horizon [k,k+N]. To do that, the estimated
system is chosen from (3.2) as (3.22)f^k+1=A¯^1,k+1(3)f¯^k+C^1,k(3)νk−1(3),f^k+j=A¯^1,k+j(3)f^k+j−1+C^k+j−1(3)νk+j−1(3);j∈[2,N],       and the
equivalent input estimates become (3.23)ν^k(3)[k]=ν^k−1(3)[k−1]=νk−1(3),(3.24)ν^k+j(3)[k]=ν^k+j(3)[k−1];allintegerk≠j,j∈[1,N],(3.25)ν^k+j(3)[k]=0;ifk=j,(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νji*(k),all i∈I,k=0,1,2,…       from (3.4) or (3.11), respectively.(f)
to obtain the suboptimal input
sequenceui*(k),alli∈I,k=0,1,2,3,…,N, 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 (3.26)fk+1=A¯k(3)fk+Ck(3)νk(3);k=0,1,2,…  instead of (3.2), and the subscript  “1”
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ν(3)(k−1)
equivalent input is applied in each optimization horizon.
(3.3.3) The usual asymptotic stability tests
viaLyapunov'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   N→∞ . 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 A¯ is a Hurwitz's  matrix and u(⋅) is
a bounded sequence, if a saturation type rule is applied when computing the
optimal redefined equivalent inputs  ν(⋅)(3)(⋅) in the various
schemes.


## 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,
(3.2)f(k+1)=A¯(3)(k)f(k)+C(3)(k)ν(3)(k);k=0,1,2,….
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
(3.3)JN=12∑k=0N[fT(k)Q¯(3)(k)f(k)+ν(3)T(k)R¯(3)(k)ν(3)(k)];N<∞.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 Theorem2.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:
(3.4)νji*(k)=νji*(3)(k+1)−[CiT(k)Q¯(k+1)Ci(k)]−1CiT(k)Q¯(k+1)A¯x¯(k)−ejiTR¯(3)−1(k+1)W¯(2)T(k+1)A¯x¯(k);i∈I;k=0,1,2,…,N,
where the optimal equivalent redefined
input sequence is
(3.5)νji(3)*(k+1)=−ejiT[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1×C(3)T(k+1)P(k+2)A¯(3)(k+1)A¯x¯(k);k=0,1,2,…,N
with  P(⋅)  being the recursive Riccati
matrices; that is, (3.6)P(k)=Q¯(3)(k)+A¯(3)(k)P(k+1)×{I−C(3)(k)[R¯(3)(k)+C(3)(k)P(k+1)C(3)(k)]−1C(3)(k)P(k+1)}A¯(3)(k);P(N)=Q¯(3)(N);k=0,1,…,N.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 Theorem3.1, consider the sequences of scalars(3.7)u¯i(k)≡{[bjiT(i)xi(k)]−1νji*(k);allj∈Ni,ifji∈Ni−,[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]−1νji*(k);allj∈Ni,ifji∈Ni+,0;ifbjiT(i)xi(k)=0,pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)=0∀j∈Ni;k=0,1,2,…,N
with νji*(k),k=0,1,2,…   being the jith       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 ui*(k)=ui(k),k=0,1,2,…  if
the interval  (−∞,∞)⊂ℝ is allowed as input definition domain, alli∈I.(b)
If|ui(k)|≤M<∞,  alli∈I, then the optimal input sequence with
respect to the loss function (3.1) is ui*=sat[u¯i(k)],  alli∈I,k=0,1,2,…; where the saturation function is defined by
(3.8)sat[g(x)]≡{Msgn[g(x)];if|g(x)|>M,g(x);otherwise.Proof.
It follows from Theorem3.1 and the fact
that, from (3.7), propositions (a)-(b) in Corollary 3.2 imply that  u(⋅)(⋅)  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 differentSi through the states of the remaining Sj,i,j∈I.
(3.1.2) Also, note that (2.7) can be rewritten as(3.9)x¯(k+1)=A¯x¯(k)+∑j=1pHj(k)uj(k);k=0,1,2,…,
where
(3.10)Hi(k)≡B¯ix¯(k);alli∈I;k=0,1,2,….
By applying optimization techniques to
(3.9), one could obtain, at first, an optimal u(⋅)(⋅)-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 thec(⋅)-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 Theorem3.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 Theorem3.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 theC(⋅)(⋅)-vectors
by their values by using each current value at the first sample of each
optimization horizon; that is, C(t)=C(k),  allk≤t≤k+N;k=0,1,2,….


### 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   ξ(⋅)   on the  ν(⋅)(3)(⋅)-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, η(t)(costate)≡P(t)f(t)−ρ(t).
Thus, by taking into account that now f(k)≡A¯x¯(k−1)+ξ(k−1), 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 Theorem3.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(3.11)νji*(k)=νji(3)*(k+1)−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)−ejiTR¯(3)−1(k+1)W¯(2)T(k+1)[A¯x¯(k)+ξ(k)];alli∈I;k=0,1,2,…, where the optimal equivalent redefined input
sequence is
(3.12)νji*(3)(k)=ejiT[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)T(k)×[ρ(k+1)−P(k+1){A(3)(k)[A¯x¯(k−1)+ξ(k−1)]+ξ(k)}];alli∈I;k=1,2,3,….
The recursive Riccati-matrix sequence is now given by(3.13)P(k)=Q¯(3)(k)+A¯(3)(k)×{P(k+1)−P(k+1)C(3)(k)[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)(k)P(k+1)}A¯(3)(k);P(N)=Q¯(3)(N);k=0,1,2,…,N,(3.14)ρ(k+1)={I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C¯(3)(k)C¯(3)T(k)A¯(3)−T(k)}−A¯(3)−T(k){ρ(k)+[P(k)−Q¯(3)(k)]A¯(3)−1(k)ξ(k)};ρ(N)=0;k=0,1,2,…,N.Proof.
It follows immediately from applying the three transformations of variables of AppendixB
to (2.23) and (A.2), and from the fact that f(k)≡A¯x¯(k−1)+ξ(k−1),
or equivalently
(3.15)f(k+1)=A¯(3)(k)f(k)+C(3)(k)ν(3)(k)+ξ(k);k=0,1,2,…,
and the fact that the  C(3)(⋅)
vectors are now
(3.16)C(3)(k)=C(3)≡[A¯ej1+l1−1,A¯ej2+l2−1,…,A¯ejp+lp−1]withk=0,1,2,3,….Remark 3.5.
(3.2.1) Note that Corollary3.2
also applies here to obtain the optimal system input sequence under constraints
(see (3.7) and (3.8)).
(3.2.2) Theρ(⋅)-vector
(related to the deterministic disturbance in the modified Riccati
transformation is dependent on the   ν(3)(⋅)-sequence. To
solve such a dependence, the ρ(k+1)-vector
of (3.14) can be decomposed into (3.17)ρ(k+1)=m(k)ρ(k)+M¯(k)ξ(k);k=0,1,2,…, being M(⋅), M¯(⋅)∈R(n+a)x(n+a) matrices defined by(3.18)M(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1;k=0,1,2,3,…,M¯(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1×A¯(3)−T(k)[P(k)−Q¯(3)(k)]A¯(3)−1(k);k=0,1,2,3,….
Thus, the following result yields.Corollary 3.6 (useful implementability result for obtaining the optimal equivalent input).
If
the same assumptions from Theorem3.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.19)νji*(k)={1+ejiT{[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)P(k+2)A¯(3)(k+1)+R¯(3)−1(k+1)W¯(2)(k+1)F(k+1)+[ejiTQ¯(k+1)eji]−1}ejiQ¯(k+1)F(k+1)}−1×{ejiT[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)×{[M¯(k+1)−P(k+2)]ξ(k+1)+M(k+1)ρ(k+1)−P(k+2)A¯(3)(k+1)A¯x¯(k)}−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)A¯x¯(k)−ejiTR¯(3)−1(k+1)W¯(2)T(k)A¯x¯(k)};alli∈I;k=0,1,2,3,…,where the matrix F(⋅)∈R(n+a)×p
is partitioned as (3.20)F(k)≡[F1(k),F2(k),…,Fp(k)];k=0,1,2,3,…. So that each vector  Fi(k)∈Rn+a(3.21)Fi(k)≡∑i=1p∑j=1,j≠ini+ai{bji−niT(i)xi(k)δji+[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]δ¯ji}−1×{bj−niT(i)xi(k)δji+[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯ji};alli∈I,k=0,1,2,….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.2.1. Centralized Suboptimal Modified Scheme 1

It consists of
the following variations on the optimal scheme from Theorem3.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 theC(⋅)(⋅)-vectors
by their values by using each current value at the first sample of each
optimization horizon; that is, C(t)=C(k),  allk≤t≤k+N;k=0,1,2,….


## 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   ξ(⋅)   on the  ν(⋅)(3)(⋅)-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, η(t)(costate)≡P(t)f(t)−ρ(t).
Thus, by taking into account that now f(k)≡A¯x¯(k−1)+ξ(k−1), 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 Theorem3.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(3.11)νji*(k)=νji(3)*(k+1)−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)−ejiTR¯(3)−1(k+1)W¯(2)T(k+1)[A¯x¯(k)+ξ(k)];alli∈I;k=0,1,2,…, where the optimal equivalent redefined input
sequence is
(3.12)νji*(3)(k)=ejiT[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)T(k)×[ρ(k+1)−P(k+1){A(3)(k)[A¯x¯(k−1)+ξ(k−1)]+ξ(k)}];alli∈I;k=1,2,3,….
The recursive Riccati-matrix sequence is now given by(3.13)P(k)=Q¯(3)(k)+A¯(3)(k)×{P(k+1)−P(k+1)C(3)(k)[R¯(3)(k)+C(3)T(k)P(k+1)C(3)(k)]−1C(3)(k)P(k+1)}A¯(3)(k);P(N)=Q¯(3)(N);k=0,1,2,…,N,(3.14)ρ(k+1)={I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C¯(3)(k)C¯(3)T(k)A¯(3)−T(k)}−A¯(3)−T(k){ρ(k)+[P(k)−Q¯(3)(k)]A¯(3)−1(k)ξ(k)};ρ(N)=0;k=0,1,2,…,N.Proof.
It follows immediately from applying the three transformations of variables of AppendixB
to (2.23) and (A.2), and from the fact that f(k)≡A¯x¯(k−1)+ξ(k−1),
or equivalently
(3.15)f(k+1)=A¯(3)(k)f(k)+C(3)(k)ν(3)(k)+ξ(k);k=0,1,2,…,
and the fact that the  C(3)(⋅)
vectors are now
(3.16)C(3)(k)=C(3)≡[A¯ej1+l1−1,A¯ej2+l2−1,…,A¯ejp+lp−1]withk=0,1,2,3,….Remark 3.5.
(3.2.1) Note that Corollary3.2
also applies here to obtain the optimal system input sequence under constraints
(see (3.7) and (3.8)).
(3.2.2) Theρ(⋅)-vector
(related to the deterministic disturbance in the modified Riccati
transformation is dependent on the   ν(3)(⋅)-sequence. To
solve such a dependence, the ρ(k+1)-vector
of (3.14) can be decomposed into (3.17)ρ(k+1)=m(k)ρ(k)+M¯(k)ξ(k);k=0,1,2,…, being M(⋅), M¯(⋅)∈R(n+a)x(n+a) matrices defined by(3.18)M(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1;k=0,1,2,3,…,M¯(k)=[I+[Q¯(3)(k)−P(k)]A¯(3)−1(k)R¯(3)−1(k)C(3)(k)C(3)T(k)A¯(3)−T(k)]−1×A¯(3)−T(k)[P(k)−Q¯(3)(k)]A¯(3)−1(k);k=0,1,2,3,….
Thus, the following result yields.Corollary 3.6 (useful implementability result for obtaining the optimal equivalent input).
If
the same assumptions from Theorem3.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.19)νji*(k)={1+ejiT{[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)P(k+2)A¯(3)(k+1)+R¯(3)−1(k+1)W¯(2)(k+1)F(k+1)+[ejiTQ¯(k+1)eji]−1}ejiQ¯(k+1)F(k+1)}−1×{ejiT[R¯(3)(k+1)+C(3)T(k+1)P(k+2)C(3)(k+1)]−1C(3)T(k+1)×{[M¯(k+1)−P(k+2)]ξ(k+1)+M(k+1)ρ(k+1)−P(k+2)A¯(3)(k+1)A¯x¯(k)}−[ejiTQ¯(k+1)eji]−1ejiTQ¯(k+1)A¯x¯(k)−ejiTR¯(3)−1(k+1)W¯(2)T(k)A¯x¯(k)};alli∈I;k=0,1,2,3,…,where the matrix F(⋅)∈R(n+a)×p
is partitioned as (3.20)F(k)≡[F1(k),F2(k),…,Fp(k)];k=0,1,2,3,…. So that each vector  Fi(k)∈Rn+a(3.21)Fi(k)≡∑i=1p∑j=1,j≠ini+ai{bji−niT(i)xi(k)δji+[pji−niT(i)zi(k)+∑l=1pmji−niT(i,l)y1(k)]δ¯ji}−1×{bj−niT(i)xi(k)δji+[pj−niT(i)zi(k)+∑l=1pmj−niT(i,l)y1(k)]δ¯ji};alli∈I,k=0,1,2,….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,[k,k+N],  1≤N<∞,   integers,(b)
to test the system controllability according
to Theorem2.6 and/or Corollary3.2,(c)
to estimate the state vector (suboptimal
modified scheme 1 of  Section3.2.1)
and the equivalent deterministic disturbance (suboptimal modified scheme 2 of Section
3.2.2) on the optimization horizon [k,k+N]. To do that, the estimated
system is chosen from (3.2) as (3.22)f^k+1=A¯^1,k+1(3)f¯^k+C^1,k(3)νk−1(3),f^k+j=A¯^1,k+j(3)f^k+j−1+C^k+j−1(3)νk+j−1(3);j∈[2,N],       and the
equivalent input estimates become (3.23)ν^k(3)[k]=ν^k−1(3)[k−1]=νk−1(3),(3.24)ν^k+j(3)[k]=ν^k+j(3)[k−1];allintegerk≠j,j∈[1,N],(3.25)ν^k+j(3)[k]=0;ifk=j,(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νji*(k),all i∈I,k=0,1,2,…       from (3.4) or (3.11), respectively.(f)
to obtain the suboptimal input
sequenceui*(k),alli∈I,k=0,1,2,3,…,N, 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 (3.26)fk+1=A¯k(3)fk+Ck(3)νk(3);k=0,1,2,…  instead of (3.2), and the subscript  “1”
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ν(3)(k−1)
equivalent input is applied in each optimization horizon.
(3.3.3) The usual asymptotic stability tests
viaLyapunov'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   N→∞ . 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 A¯ is a Hurwitz's  matrix and u(⋅) is
a bounded sequence, if a saturation type rule is applied when computing the
optimal redefined equivalent inputs  ν(⋅)(3)(⋅) in the various
schemes.


## 4. Numerical Simulation

### 4.1. Example 1

This first example
implemented deals with the following ninth-order discrete model:(4.1)S1:x1(k+1)≡[A1(T)+u1(k)B1(T)]x1(k)+D1(T)z1(k);u1∈R,y1(k)≡x1(k);x1(⋅)∈R2x1,H1:z1(k+1)≡[M1(T)+u1(k)P1(T)]z1(k)+L11(T)x1(k);z1(⋅)∈R2x1,H2:z2(k+1)≡M2(T)z2(k)+L21(T)x1(k);z2(⋅)∈R3x1,H3:z3(k+1)≡M3(T)z3(k)+L31(T)x1(k);z3(⋅)∈R2x1,allk≥0,k∈z. Table 1 displays
the entries of the constant discrete matrices related to the composite
homogeneous bilinear structures  S1, Hi,  i=1,2,3,
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.
A1≡[1−0.1557T001−0.1800T];  B1≡[00−0.0189T0];   D1≡[000(−0.0557+0.0043T)T(−0.0557+0.0043T)T00000]M1≡[1−100.0000T013.3333T1−85.8333T];   P1≡[0(80.0000−8000.0000T)T01066.6640T];
L11≡[00(−0.0062+0.5322T)T0]M2≡[1−28.6145T1950.3287T35.7332T1.9554T1−231.1652T−4.1281T0.0260T−0.0056T1+0.1706T];  L21≡[(35.7322−4536.9199T)T0(−4.1281+512.0844T)T0(−0.0056+0.0575T)T0]M3≡[1−2.0000T06.8570T1−75.2380T]; L21≡[00(−6.9022+44.5212T)T0]The
steady state of the signal trajectory under no control leads to the following
parameters:y1,st*=21.8, y2,st*=238.63, ust*=400,  Nst*=57     samples. The system performance is studied as a function of the optimization
horizon, N. Also, the plant settings are: Initial conditions, x1(0)=[20.0,250.0]T,  z1(0)=[0.0,0.0]T,  z2(0)=[0.0,0.0,0.0]T,    and z3(0)=[0.0,0.0]T;    sampling period T=0.03   seconds and working horizon WH=150   samples.The
obtained results are shown in Table2 where the deviations (percent) are
related to the normal operation signals (i.e.,  y1,st*,    y2,st*   and  ust*) 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.
Steady-state responseFirst componentSecond componentSteady-statey1,sty2,stcontrol responseWaveform plotOptimizaton horizonN (samples)Steady-state responseNst (samples)ValuePercent deviationValuePercent deviationUstPercent deviationConst functionJNCCT (seconds)—–14722.51+ 3.26237.51− 0.47405.38+ 1.3515392.761.132.....23522.33+ 2.43238.15− 0.2400.43+ 0.1212468.191.121---33222.14+ 1.56237.32− 0.55403.61+ 0.9210319.780.991xxx42821.87+ 0.32236.49− 0.9398.27− 0.4314215.841.001none53222.05+ 1.15238.01− 0.26399.15− 0.2113854.210.997none63321.14− 3.03238.18− 0.19382.54− 4.3714129.830.946—73520.25− 7.11225.06− 5.09375.82− 6.0515639.080.94none84120.03− 8.12210.87− 11.63366.39− 8.410378.510.984none94519.98− 8.35205.17− 14.02351.41− 9.6513227.831............104819.81− 9.13201.39− 15.61360.23− 9.9416863.710.945none155915.74− 23.21195.43− 18.1328.51− 17.8717459.290.954none207615.92− 26.97187.31− 21.51310.32− 22.4229243.620.935none1011513.61−37.57132.75− 44.37200.39− 49.958346.430.957Exhaustive 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 Figure2,
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


### 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 Section2.2, by (4.2)x(k+1)=[A(T)+u(k)N(T)]x(k)+E(T)u(k)+w(k),y(k)=Cx(k),  where x(k)≡[x1T(k),x2T(k),x3T(k)]T,   y(k)≡[y1T(k),0,0]T     and there exists a deterministic disturbance  w(k)≡[w1T(k),0,0]T
so that w1(k)≡S1(T)u2(k)+[S2(T)+u(k)S3(k)]q(k),
being T  the sampling period,  (⋅)T   denotes transpose
and the sampling instant is defined as  tk+1=tk+T.In (4.2), the
coefficients of the block matrices are composed of the following structures: (4.3)A(T)=[A1(T)0M1(T)P2(T)A2(T)R2(T)P3(T)R3(T)A3(T)];N(T)=[N1(T)0S2(T)000000];E(T)=[E1(T)E2(T)E3(T)];C=[I00]. 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.Figure4 plots a typical tracking response output of the system as a
function of the sampling period, for a fixed optimization horizon N=50   samples and  Q¯   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.


## 4.1. Example 1

This first example
implemented deals with the following ninth-order discrete model:(4.1)S1:x1(k+1)≡[A1(T)+u1(k)B1(T)]x1(k)+D1(T)z1(k);u1∈R,y1(k)≡x1(k);x1(⋅)∈R2x1,H1:z1(k+1)≡[M1(T)+u1(k)P1(T)]z1(k)+L11(T)x1(k);z1(⋅)∈R2x1,H2:z2(k+1)≡M2(T)z2(k)+L21(T)x1(k);z2(⋅)∈R3x1,H3:z3(k+1)≡M3(T)z3(k)+L31(T)x1(k);z3(⋅)∈R2x1,allk≥0,k∈z. Table 1 displays
the entries of the constant discrete matrices related to the composite
homogeneous bilinear structures  S1, Hi,  i=1,2,3,
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.
A1≡[1−0.1557T001−0.1800T];  B1≡[00−0.0189T0];   D1≡[000(−0.0557+0.0043T)T(−0.0557+0.0043T)T00000]M1≡[1−100.0000T013.3333T1−85.8333T];   P1≡[0(80.0000−8000.0000T)T01066.6640T];
L11≡[00(−0.0062+0.5322T)T0]M2≡[1−28.6145T1950.3287T35.7332T1.9554T1−231.1652T−4.1281T0.0260T−0.0056T1+0.1706T];  L21≡[(35.7322−4536.9199T)T0(−4.1281+512.0844T)T0(−0.0056+0.0575T)T0]M3≡[1−2.0000T06.8570T1−75.2380T]; L21≡[00(−6.9022+44.5212T)T0]The
steady state of the signal trajectory under no control leads to the following
parameters:y1,st*=21.8, y2,st*=238.63, ust*=400,  Nst*=57     samples. The system performance is studied as a function of the optimization
horizon, N. Also, the plant settings are: Initial conditions, x1(0)=[20.0,250.0]T,  z1(0)=[0.0,0.0]T,  z2(0)=[0.0,0.0,0.0]T,    and z3(0)=[0.0,0.0]T;    sampling period T=0.03   seconds and working horizon WH=150   samples.The
obtained results are shown in Table2 where the deviations (percent) are
related to the normal operation signals (i.e.,  y1,st*,    y2,st*   and  ust*) 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.
Steady-state responseFirst componentSecond componentSteady-statey1,sty2,stcontrol responseWaveform plotOptimizaton horizonN (samples)Steady-state responseNst (samples)ValuePercent deviationValuePercent deviationUstPercent deviationConst functionJNCCT (seconds)—–14722.51+ 3.26237.51− 0.47405.38+ 1.3515392.761.132.....23522.33+ 2.43238.15− 0.2400.43+ 0.1212468.191.121---33222.14+ 1.56237.32− 0.55403.61+ 0.9210319.780.991xxx42821.87+ 0.32236.49− 0.9398.27− 0.4314215.841.001none53222.05+ 1.15238.01− 0.26399.15− 0.2113854.210.997none63321.14− 3.03238.18− 0.19382.54− 4.3714129.830.946—73520.25− 7.11225.06− 5.09375.82− 6.0515639.080.94none84120.03− 8.12210.87− 11.63366.39− 8.410378.510.984none94519.98− 8.35205.17− 14.02351.41− 9.6513227.831............104819.81− 9.13201.39− 15.61360.23− 9.9416863.710.945none155915.74− 23.21195.43− 18.1328.51− 17.8717459.290.954none207615.92− 26.97187.31− 21.51310.32− 22.4229243.620.935none1011513.61−37.57132.75− 44.37200.39− 49.958346.430.957Exhaustive 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 Figure2,
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


## 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 Section2.2, by (4.2)x(k+1)=[A(T)+u(k)N(T)]x(k)+E(T)u(k)+w(k),y(k)=Cx(k),  where x(k)≡[x1T(k),x2T(k),x3T(k)]T,   y(k)≡[y1T(k),0,0]T     and there exists a deterministic disturbance  w(k)≡[w1T(k),0,0]T
so that w1(k)≡S1(T)u2(k)+[S2(T)+u(k)S3(k)]q(k),
being T  the sampling period,  (⋅)T   denotes transpose
and the sampling instant is defined as  tk+1=tk+T.In (4.2), the
coefficients of the block matrices are composed of the following structures: (4.3)A(T)=[A1(T)0M1(T)P2(T)A2(T)R2(T)P3(T)R3(T)A3(T)];N(T)=[N1(T)0S2(T)000000];E(T)=[E1(T)E2(T)E3(T)];C=[I00]. 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.Figure4 plots a typical tracking response output of the system as a
function of the sampling period, for a fixed optimization horizon N=50   samples and  Q¯   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.



---

*Source: 817063-2008-04-02.xml*