Modular Feedback Control of Networked Systems by Clustering: A Drinking Water Network Case Study

Abstract

This article presents a method based on linear matrix inequalities (LMIs) for designing a modular feedback control law, whose synthesis guarantees the system stability, while switching to different network topologies. Such stability is achieved by means of a common Lyapunov function to all network admissible configurations. Several mechanisms to relieve the computational burden of this methodology in large-scale systems are also presented. To assess its applicability, the modular controller is tested on a real case study, namely the Barcelona drinking water network (DWN), and its performance is compared with that of other control strategies, showing the effectiveness of the proposed approach.

1. Introduction

In recent years, distributed control architectures have gained relevance due to the multiple advances in information and communication technologies [1]. In general, distributed control systems are characterized by independent and interacting subsystems governed by controllers that exchange information to get additional performance by coordinating their control actions. This approach is useful in many practical problems that cannot be addressed from a centralized perspective, e.g., due to the sheer size of the system and/or limits in the information exchange between controllers [2].

In this context, we are interested in the design of feedback controllers that respect the constraints imposed by the communication topology. To this end, several methods have been proposed in the literature. For example, in [3], a gradient method for considering sparsity constraints in linear quadratic regulators is implemented by a game-theoretic algorithm. In [4], communication constraints are imposed via sparsity-promoting penalty functions on the cardinality of the communication links used in the control architecture. In [5], the notion of quadratic invariance of a set of sparsity or delay constraints on the feedback controller is introduced as a means to guarantee a convex design problem. Also, a similar strategy was recently used in [6] for large-scale systems. Linear matrix inequality (LMI)-based approaches for feedback controllers can also be found in the literature. For instance, a design method based on LMIs for distributed linear systems is proposed in the context of coalitional control first in [7] and later in [8]. Under this framework, communication links remain enabled as long as they provide a significant performance increase, and otherwise they are disconnected. As a consequence, controllers are grouped dynamically into disjoint cooperating sets, which are called coalitions. In this way, a trade-off between performance and communication burden is obtained at the cost of a more complex implementation. In particular, the design method seeks a feedback gain K with as many zero entries as possible. Some other works also propose LMI-based methodologies [9,10,11] to design sparse feedback matrices. Similarly, $H_{\infty }$ techniques are used in [12] to achieve a structured feedback control law. It is also remarkable that these design constraints can also be considered to compute other classical controllers by using a feedback gain, e.g., a proportional-integral controller for an irrigation canal is calculated in [13] in this manner.

Here, we apply the modular control foundations presented in [14], which rely on the principles introduced in [7,8]. The key idea of this methodology is that the feedback controllers used in the different topologies are created considering a common template composed of pieces associated with the different communication links in the network. Hence, if a link is disabled the corresponding elements in the controller are merely replaced by zeros without affecting the rest of the elements. As a consequence, given two different network topologies, the nonzero elements of the corresponding feedback controllers share the same values. Therefore, the feedback matrix shall easily change its internal configuration depending on the topology to be controlled. Furthermore, the modular control law guarantees the stability by means of a common Lyapunov function to all network topologies.

Undoubtedly, this control strategy can be of interest for the aforementioned framework of coalitional control, where the information structure of the system plays an essential role [15,16,17,18,19]. Natural application fields for this type of online partitioning approach are, e.g., traffic [20], water [21,22,23], cellular [24,25] and power networks [26], and renewable energy generation systems [27,28]. In particular, the approach presented here may be suitable to be combined with the game-theoretic methods in [29,30,31], which are also based on LMIs, and with plug-and-play control capabilities [32,33,34]. Also, controllers designed with the modular structure mentioned above are interesting for applications where communications between controllers can fail, e.g., due to packet losses and jamming attacks [35,36,37]. The rationale is simple: any missing information from a neighbor can be considered to be a disabled communication link, which allows applying the results of this article. Likewise, model predictive control (MPC)-based works could use the proposed controller as stabilizing terminal feedback, especially for coupling-dependant clustering architectures, e.g., in [38,39,40], and also including most of the coalitional works commented before.

The main contribution of this paper is a methodology to design a feedback gain suitable for the control of networked systems by clustering [14], which is applied here to complex large-scale systems. Certainly, modular controller synthesis for large-scale schemes is challenging due to the computational complexity. Different ways to tackle the computational explosion while guaranteeing the system stability are shown in this article. In particular, an approach based on clustering of agents has been implemented in the Barcelona drinking water network (DWN) as a case study, where the model deployed in [41,42,43] has been partitioned into a network of eight agents following [44]. In fact, the need for systematic methods to achieve the partitioning objective has gained importance recently, with partitioning schemes based on graph theory [40,42,44,45,46,47,48], states and inputs estimation [49], social network algorithms [50], genetic algorithms [51], and PageRank [22,52]. Specific partitioning techniques applied to large-scale water systems as the proposed case study can be found in [42,43,44,50,53,54,55]. Indeed, control applications to water distribution systems are becoming more common, with recent contributions in the monitoring and control of valves leakage [56], pumps speed [57], pressure management by clustering [23], or pump scheduling [58].

The outline of the rest of the paper is organized as follows. Section 2 presents the problem statement. Section 3 introduces the concept of modular control, providing an LMI-based design method together with some properties of interest. Section 4 proposes different techniques to tackle the computational explosion for large-scale systems. Section 5 introduces the application of modular feedback control to the Barcelona DWN case study. Finally, concluding remarks are given in Section 6.

2. Problem Formulation

We consider an overall discrete linear time-invariant system, which is composed of a set $\mathcal{N}=\{1,2,\dots ,N\}$ of subsystems characterized as

$$\begin{array}{ccc}\hfill x_i(k+1)& =& A_{ii}{x}_i\left(k\right)+B_{ii}{u}_i\left(k\right)+{\widehat{d}}_i\left(k\right),\hfill \\ \hfill {\widehat{d}}_i\left(k\right)& =& \sum _{j\ne i}{A}_{ij}{x}_j\left(k\right)+\sum _{j\ne i}{B}_{ij}{u}_j\left(k\right),\hfill \end{array}$$

(1)

where $x_i\in {\mathbb{R}}^{n_{x_i}}$ and $u_i\in {\mathbb{R}}^{n_{u_i}}$, with $i=1,\dots ,n$, are the states and inputs of each subsystem $i\in \mathcal{N}$, respectively, with $n_{x_i}$ and $n_{u_i}$ denoting the size of vectors $x_i$ and $u_i$. Matrices $A_{ii}\in {\mathbb{R}}^{n_{x_i}\times n_{x_i}}$ and $B_{ii}\in {\mathbb{R}}^{n_{x_i}\times n_{u_i}}$ refer to the state and input-to-state matrices, and ${\widehat{d}}_i\in {\mathbb{R}}^{n_{x_i}}$ represents the influence of the neighboring states and inputs in the update of $x_i$. Finally, matrices $A_{ij}\in {\mathbb{R}}^{n_{x_i}\times n_{x_j}}$ and $B_{ij}\in {\mathbb{R}}^{n_{x_i}\times n_{u_j}}$ map the states and inputs of subsystem $j\in \mathcal{N}$ into the states of subsystem i, respectively.

The goal of the subsystems in $\mathcal{N}$ is to minimize the following stage cost:

$${\ell }_i\left(k\right)=x_i^{\mathrm{T}}\left(k\right)Q_i{x}_i\left(k\right)+u_i^{\mathrm{T}}\left(k\right)R_i{u}_i\left(k\right),$$

(2)

where $Q_i\in {\mathbb{R}}^{n_{x_i}\times n_{x_i}}$ and $R_i\in {\mathbb{R}}^{n_{u_i}\times n_{u_i}}$ are positive semi-definite and definite constant weighting matrices, respectively.

From a global viewpoint, the overall dynamics are simply described by

$$x_{\mathcal{N}}(k+1)=A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+B_{\mathcal{N}}{u}_{\mathcal{N}}\left(k\right),$$

(3)

where subscript $\mathcal{N}$ emphasizes that all system vectors and matrices come from the aggregation of local subsystems, i.e., $x_{\mathcal{N}}={\left[x_i\right]}_{i\in \mathcal{N}}$, $u_{\mathcal{N}}={\left[u_i\right]}_{i\in \mathcal{N}}$, $A_{\mathcal{N}}={\left[A_{ij}\right]}_{i,j\in \mathcal{N}}$, and $B_{\mathcal{N}}={\left[B_{ij}\right]}_{i,j\in \mathcal{N}}$. For convenience, we will respectively denote by $n_x$ and $n_u$ the number of states and controls of the overall system. Please note that in the global model there are no neighbors disturbances, because mutual interactions are already included in (3). Likewise, the stage cost of the overall system can be expressed as a function of states and inputs of the corresponding subsystems, i.e.,

$${\ell }_{\mathcal{N}}\left(k\right)=x_{\mathcal{N}}^{\mathrm{T}}\left(k\right)Q_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+u_{\mathcal{N}}^{\mathrm{T}}\left(k\right)R_{\mathcal{N}}{u}_{\mathcal{N}}\left(k\right),$$

(4)

where $Q_{\mathcal{N}}=\mathrm{diag}{\left(Q_i\right)}_{i\in \mathcal{N}}$ and $R_{\mathcal{N}}=\mathrm{diag}{\left(R_i\right)}_{i\in \mathcal{N}}$.

2.1. Modular Control Law and Communication Constraints

A linear feedback controller is proposed here to minimize the cost of the system while steering it towards the origin. Also, it must be designed taking into account constraints in the information flows due to the communication network that connects the subsystems. In particular, we assume that the network is described by a directed graph ($\mathcal{N}$, $\mathcal{L}$), where $\mathcal{N}$ is a set of subsystems and $\mathcal{L}$ is a set of unidirectional links given by $\mathcal{L}\subseteq {\mathcal{L}}^{\mathcal{N}}=\left\{\{i,j\}\right|\{i,j\}\subseteq \mathcal{N},i\ne j\}$. For convenience, we will define link $\{i,j\}$, or simply $l_{ij}$, as an arrow that goes from j to i and vice versa, with $i,j\in \mathcal{N}$, to stress that i receives information from j. Likewise, we assume that links allow only direct communication, i.e., agent i receives information from j only if they are directly connected by $l_{ij}$, although this assumption can be relaxed if needed. For representation simplicity, we symbolize two links in opposite directions with a double arrow. In this work, it is also considered that some communication links might be disabled to reduce the communication burden, due to, e.g., jamming attacks. To this end, let us define $\mathsf{\Lambda }$ as the topology given by the set of active links in $\mathcal{L}$. Given that there are $L=\left|\mathcal{L}\right|$ links, it is possible to define a set $\mathcal{T}=\{{\mathsf{\Lambda }}_{\mathrm{DC}},{\mathsf{\Lambda }}_1,{\mathsf{\Lambda }}_2,\dots ,{\mathsf{\Lambda }}_{2^{L-2}},{\mathsf{\Lambda }}_{\mathcal{L}}\}$ composed of $2^L$ topologies. Please note that we have introduced a slightly different notation for two special topologies, namely ${\mathsf{\Lambda }}_{\mathrm{DC}}$ and ${\mathsf{\Lambda }}_{\mathcal{L}}$, which correspond to the decentralized topology (all links are disabled) and the full communication topology (all links are enabled), respectively.

Since the feedback controller must be designed taking into account the constraints in the information flows imposed by topology $\mathsf{\Lambda }$, a superscript will be added to stress this fact so that $u_{\mathcal{N}}=K_{\mathsf{\Lambda }}{x}_{\mathcal{N}}$, with $\mathsf{\Lambda }\in \mathcal{T}$. More specifically, the control law must be suitable for any $\mathsf{\Lambda }\in \mathcal{T}$. Unlike [8], where different feedback gains are calculated for each possible topology, here we propose a controller composed of blocks associated with the links of the communication network [14]. To illustrate this idea, Figure 1 shows a system composed of four agents that communicate using 12 directed links, which are associated with the nondiagonal elements of the corresponding modular feedback controller for the full communication topology ${\mathsf{\Lambda }}_{\mathcal{L}}$, which is described by

$$K_{{\mathsf{\Lambda }}_{\mathcal{L}}}=\left[\begin{array}{cccc}{K}_{11}& K_{12}& K_{13}& K_{14}\\ K_{21}& K_{22}& K_{23}& K_{24}\\ K_{31}& K_{32}& K_{33}& K_{34}\\ K_{41}& K_{42}& K_{43}& K_{44}\end{array}\right].$$

(5)

However, if the topology changes and becomes that of Figure 2, say ${\mathsf{\Lambda }}_{\mathrm{ex}}$, with five links omitted with respect to ${\mathsf{\Lambda }}_{\mathcal{L}}$, the control law becomes

$$K_{{\mathsf{\Lambda }}_{\mathrm{ex}}}=\left[\begin{array}{cccc}{K}_{11}& 0& K_{13}& K_{14}\\ K_{21}& K_{22}& 0& 0\\ 0& K_{32}& K_{33}& 0\\ K_{41}& K_{42}& K_{43}& K_{44}\end{array}\right],$$

(6)

where all nonzero entries have the same value that they had in $K_{{\mathsf{\Lambda }}_{\mathcal{L}}}$. Hence, modular feedback controller $K_{{\mathsf{\Lambda }}_{\mathcal{L}}}$ provides us with a family of control laws that can be adapted to different communication topologies by simply making zero the elements that correspond to disabled links. Please note that this modular structure adds robustness to the feedback laws designed through the proposed approach even if the system partitioning has been poorly performed because it rearranges the partitioning online starting from the set of atomic components in which the system is divided. Hence, an inadequate initial partitioning is not an issue because if two subsystems need to cooperate, the controller will promote the cooperation as long as the benefits expected are greater than the corresponding cooperation efforts.

Finally, some definitions regarding parents-children topologies are introduced below and illustrated in Figure 3.

Definition 1.

A set of ascendant topologies coming from a given topology ${\mathsf{\Lambda }}_i\in \mathcal{T}$, with $|{\mathsf{\Lambda }}_i|<L$, can be defined as

$$S_{{\mathsf{\Lambda }}_i}=\{{\mathsf{\Lambda }}_j\in \mathcal{T}\mid {\mathsf{\Lambda }}_i\subset {\mathsf{\Lambda }}_j,|{\mathsf{\Lambda }}_j|\ge |{\mathsf{\Lambda }}_i|+1\}.$$

(7)

A topology ${\mathsf{\Lambda }}_j\in S_{{\mathsf{\Lambda }}_i}$ is a parent of topology ${\mathsf{\Lambda }}_i$ if $|{\mathsf{\Lambda }}_j|=|{\mathsf{\Lambda }}_i|+1$ holds.

Definition 2.

A set of descendant topologies coming from a given topology ${\mathsf{\Lambda }}_i\in \mathcal{T}$, with $|{\mathsf{\Lambda }}_i|>1$ can be defined as

$${\overline{S}}_{{\mathsf{\Lambda }}_i}=\{{\mathsf{\Lambda }}_j\in \mathcal{T}\mid {\mathsf{\Lambda }}_j\subset {\mathsf{\Lambda }}_i,|{\mathsf{\Lambda }}_j|\le |{\mathsf{\Lambda }}_i|-1\}.$$

(8)

A topology ${\mathsf{\Lambda }}_j\in {\overline{S}}_{{\mathsf{\Lambda }}_i}$ is a child of topology ${\mathsf{\Lambda }}_i$ if $|{\mathsf{\Lambda }}_j|=|{\mathsf{\Lambda }}_i|-1$ is fulfilled.

2.2. Stability

Before addressing the controller design procedure, it is necessary to give some remarks regarding stability. Since the network topology might change at any time, there is a need for guaranteeing the stability of the closed-loop system despite the switchings between the corresponding control laws $K_{\mathsf{\Lambda }}$, with $\mathsf{\Lambda }\in \mathcal{T}$. To deal with this issue, a common Lyapunov function $f\left(x_{\mathcal{N}}\left(k\right)\right)=x_{\mathcal{N}}^{\mathrm{T}}\left(k\right)P{x}_{\mathcal{N}}\left(k\right)$ is designed for all feedback controllers $K_{\mathsf{\Lambda }}$. In particular, let $P\in {\mathbb{R}}^{n_x\times n_x}$ be a positive definite matrix that will also provide us with a bound on the cost-to-go of the closed loop system, i.e.,

$$x_{\mathcal{N}}^{\mathrm{T}}\left(k\right)P{x}_{\mathcal{N}}\left(k\right)\ge \sum _{t=k}^{\infty }{\ell }_{\mathcal{N}}\left(t\right).$$

(9)

In [8], it is shown that the Lyapunov function is a bound on the cost-to-go of the closed loop system if the following inequality holds:

$$x_{\mathcal{N}}^{\mathrm{T}}\left(k\right)P{x}_{\mathcal{N}}\left(k\right)\ge {\ell }_{\mathcal{N}}\left(k\right)+x_{\mathcal{N}}^{\mathrm{T}}(k+1)P{x}_{\mathcal{N}}(k+1).$$

(10)

In particular, (9) can be derived from (10) by telescope summation.

3. Modular Controller Design

In this section, we provide a theorem and a lemma for the modular controller design with a stage cost defined by $Q_{\mathcal{N}}$ and $R_{\mathcal{N}}$ following (4). To this end, we consider that the system is decomposed into a set $\mathcal{N}$ of subsystems connected by means of a set $\mathcal{L}$ of directed communication links that give rise to a set $\mathcal{T}$ of different topologies.

Theorem 1.

Let a system be described by (3) with discrete-time linear dynamics given by $A_{\mathcal{N}}$ and $B_{\mathcal{N}}$. If there exist matrices $W=W^{\mathrm{T}}=\mathrm{diag}{\left(W_i\right)}_{i\in \mathcal{N}}$, where $W_i\in {\mathbb{R}}^{n_{x_i}\times n_{x_i}}$, and $Y\in {\mathbb{R}}^{n_u\times n_x}$ such that the following constraint is satisfied, for all $\mathsf{\Lambda }\in \mathcal{T}$:

$$\begin{array}{c}\left[\begin{array}{cccc}W& W{A}_{\mathcal{N}}^{\mathrm{T}}+Y_{\mathsf{\Lambda }}^{\mathrm{T}}{B}_{\mathcal{N}}^{\mathrm{T}}& W{Q}_{\mathcal{N}}^{1/2}& Y_{\mathsf{\Lambda }}^{\mathrm{T}}{R}_{\mathcal{N}}^{1/2}\\ A_{\mathcal{N}}W+B_{\mathcal{N}}{Y}_{\mathsf{\Lambda }}& W& 0& 0\\ Q_{\mathcal{N}}^{1/2}W& 0& I& 0\\ R_{\mathcal{N}}^{1/2}{Y}_{\mathsf{\Lambda }}& 0& 0& I\end{array}\right]>0,\hfill \end{array}$$

(11)

with $Y_{\mathsf{\Lambda },ij}=Y_{ij}$ if link $l_{ij}$ is activated, i.e., if $l_{ij}\in \mathsf{\Lambda }$, and $Y_{\mathsf{\Lambda },ij}=0$ otherwise, then there exists a modular controller that provides a family of feedback control laws $K_{\mathsf{\Lambda }}=Y_{\mathsf{\Lambda }}{W}^{-1}$, which stabilizes the system for topologies $\mathsf{\Lambda }\in \mathcal{T}$. Also, a common Lyapunov function $f\left(x_{\mathcal{N}}\left(k\right)\right)=x_{\mathcal{N}}^{\mathrm{T}}\left(k\right)P{x}_{\mathcal{N}}\left(k\right)$ that provides a bound on the cost-to-go is generated, with $P=W^{-1}$.

Proof. 

The iterative application of the Schur complement [59] in a backwards manner together with the proposed variable change allows us to transform LMI (11) into (10), which guarantees that stability for the cost-to-go must decrease at each time step. Then, a telescope summation of this inequality from $t=k$ to infinity gives (9), which allows us to use the Lyapunov function to get a bound on the cost-to-go. The constraints imposed on $Y_{\mathsf{\Lambda }}$ and W guarantee that $K_{\mathsf{\Lambda }}$ satisfies the communication restrictions imposed by topology $\mathsf{\Lambda }$, which require $K_{\mathsf{\Lambda },ij}=0,\forall i,j$ such that $l_{ij}\notin \mathsf{\Lambda }$. Here, let us recall that $K_{\mathsf{\Lambda }}=Y_{\mathsf{\Lambda }}{W}^{-1}$. Since $Y_{\mathsf{\Lambda },ij}=0$ if $l_{ij}\notin \mathsf{\Lambda }$, and W has a diagonal block structure that is to be inherited by its inverse matrix, i.e., $W_{ij}^{-1}=0$ for any $i\ne j$, then by the properties of matrix multiplication it holds $K_{\mathsf{\Lambda },ij}=0$ if link $l_{ij}\notin \mathsf{\Lambda }$. Finally, as $Y_{\mathsf{\Lambda },ij}$ always has the same value for topologies $\mathsf{\Lambda }\in \mathcal{T}$ when $l_{ij}=1$, otherwise being zero, this fact forces block $K_{\mathsf{\Lambda },ij}$ to be the same for all these topologies, hence providing the desired modular structure. □

Remark 1.

Matrix W must be a block sparse/diagonal matrix to preserve the modular features. for instance, let us assume that block $K_{\mathsf{\Lambda },34}$ is associated with link $l_{34}$, being $K_{\mathsf{\Lambda },34}$ the product of the 3rd-row of $Y_{\mathsf{\Lambda }}$ and the 4th-column of $W^{-1}$. The only way to assure that $K_{\mathsf{\Lambda },34}$ is non-null for a given $Y_{\mathsf{\Lambda }}$, when link $l_{34}$ is enabled, is to shape W as a block diagonal matrix, i.e., $W=W^{\mathrm{T}}=\mathrm{diag}{\left(W_i\right)}_{i\in \mathcal{N}}$.

Solving the set of LMIs defined by Theorem 1 provides us with a matrix $P=W^{-1}$ and matrices $K_{\mathsf{\Lambda }}=Y_{\mathsf{\Lambda }}{W}^{-1}$, where $K_{\mathsf{\Lambda },ij}$ maps the contribution of the state of agent j to the control action of agent i. Hence, if a link $l_{ij}\in \mathcal{L}$ is activated, block $K_{\mathsf{\Lambda },ij}\ne 0$; otherwise $K_{\mathsf{\Lambda },ij}=0$. It must be noticed that the set of LMIs is solved for matrices W and $Y=Y_{{\mathsf{\Lambda }}_{\mathcal{L}}}$, i.e., there are only two unknown matrices to be found. This fact represents a major difference with respect to [8], where different $W_{\mathsf{\Lambda }}$ and $Y_{\mathsf{\Lambda }}$ were obtained for each topology $\mathsf{\Lambda }$. It is also important to remark that all feedback matrices $K_{\mathsf{\Lambda }}$ are associated with the same Lyapunov matrix P, what means that the closed-loop stability is guaranteed for any topology $\mathsf{\Lambda }\in \mathcal{T}$, so that switchings between the different control laws can be performed without compromising the stability of the closed-loop system.

A necessary and sufficient condition for Theorem 1 to have a solution are given by the lemma below.

Lemma 1.

If there exists a feasible solution of (11) for topology ${\mathsf{\Lambda }}_i$, then topologies ${\mathsf{\Lambda }}_j$ in $S_{{\mathsf{\Lambda }}_i}\subseteq \mathcal{T}$, i.e., the set of ascendants of ${\mathsf{\Lambda }}_i$ (see Definition 1), admit a modular solution.

Corollary 1.

The family of LMIs of Theorem 1 has a feasible solution if and only if W and $Y_{\mathsf{\Lambda }}$ in (11) can be found for ${\mathsf{\Lambda }}_{\mathrm{DC}}$, i.e., for the decentralized case.

Remark 2.

The condition in Corollary 1 is necessary because any solution satisfying Theorem 1 must provide a feedback law for the decentralized topology, i.e., $K_{{\mathsf{\Lambda }}_{\mathrm{DC}}}$. Hence, without a solution for ${\mathsf{\Lambda }}_{\mathrm{DC}}$ there is no solution for the overall problem. Furthermore, it is a sufficient condition because W and $Y_{{\mathsf{\Lambda }}_{\mathrm{DC}}}$ constitute a feasible solution for the rest of topologies, which solve the very same LMI constraint with additional degrees of freedom due to the additional nonzero elements in the corresponding $Y_{\mathsf{\Lambda }}$. Hence, the more demanding constraints are those of the decentralized topology, which requires stable global dynamics although each subsystem uses only local state information. Consequently, any set of local feedback controllers leading to overall stable dynamics could generate a modular controller applying Theorem 1.

Design Method

To design the controller, we solve

$$\underset{W,Y_{\mathsf{\Lambda }}}{max}\mathrm{tr}\left(W\right),$$

(12)

subject to (11) for all $\mathsf{\Lambda }\in \mathcal{T}$. Then, it is enough to take $K_{\mathsf{\Lambda }}=Y_{\mathsf{\Lambda }}{W}^{-1}$. Please note that the maximization of the trace of W is an indirect manner of minimizing that of $P=W^{-1}$, hence minimizing the cost-to-go of the closed-loop system.

Once the modular controller is designed, it is straightforward to find new bounds on the cost-to-go tailored to each topology, i.e., $P_{\mathsf{\Lambda }}$, if needed, e.g., for a rapid topology selection in a coalitional control system. The fact that a common Lyapunov function exists guarantees that switchings can be performed without endangering the system stability.

4. Dealing with Computational Burden

LMI (11) is solved simultaneously for all network topologies $\mathsf{\Lambda }\in \mathcal{T}$. Notice that there exists a $Y_{\mathsf{\Lambda }}$ declared for every network topology, where each block $Y_{\mathsf{\Lambda },ij}$ corresponds to an enabled link $l_{ij}$. Therefore, the number of decision variables is defined by

$$\sum _{i\in \mathcal{N}}\frac{n_{x_i}\times (n_{x_i}+1)}{2}+n_u\times n_x,$$

(13)

being independent of the number of topologies considered. Nevertheless, the computational burden of solving LMIs does not scale linearly with the number of topology constraints [60], which might render the problem unfeasible for practical use. For example, in our experiments, we have been able to apply this method to systems with thousands of topologies, which in the most conservative case corresponds to a modest number of links, e.g., in the range 10–20. To overcome this issue, we provide some plausible strategies:

• Exploiting Convexity: Since constraint (11) is convex, any convex combination of solutions is also a solution. Given that W is common for all topologies, new solutions are generated by simply combining the results for matrices $Y_{\mathsf{\Lambda }}$ corresponding to different topologies. Likewise, it is straightforward to check that the same holds for $K_{\mathsf{\Lambda }}$.
  With this idea in mind, the problem can be simplified by solving only a subset of topologies that can be used to generate the rest of the topologies. For instance, it is possible to solve the decentralized topology, and then topologies with only one active link, i.e., instead solving $2^L$ LMIs the problem is reduced to the resolution of $L+1$ LMIs, thus avoiding the combinatorial explosion. Interestingly, the resulting feedback solution will preserve the values for the common elements in the feedback gains combined. Also, the values of the noncommon elements are easy to calculate due to the block structure of the modular controller (recall that blocks corresponding to a link become zero when the link is disabled).
  This strategy can be very convenient for hierarchical control, where an upper control layer can compute a convex combination of $K_{\mathsf{\Lambda }}$ to generate a feedback gain for the desired topology, e.g., by searching for the combination that maximizes the trace of W (or minimizes the trace of P).
• Ascendants replaced by descendants: Any topology provides a feasible solution for all its ascendants (see Definition 1). Thus, it is possible to remove topologies from set $\mathcal{T}$ and simply use one of its descendants (Definition 2) instead to reduce the number of constraints of the problem.
• Branch-and-bound-like approaches: Those topologies that degrade the most the performance of the controller are removed, e.g., those that decrease the most the trace of W. Please note that given a topology, its trace of W, i.e., when the LMI is solved for this specific topology, is a lower bound for all its descendants. This fact is used to compute a modular controller only for those topologies that provide the best performance according to the aforementioned criterion.
• Topology clustering: For large-scale systems, clustering the different agents/subsystems in super agents can be an interesting approach to handle the computational burden. For instance, the agents clustering could be achieved by different algorithms in the line of those described in [28,40,42,44], where the partitioning is performed by pre-selecting those sets of agents that are not highly interrelated with other agents and may separately work well offline prior to proceed with the control law implementation.

5. Case Study: Barcelona DWN

In this article, we analyze the Barcelona drinking water network (DWN), which is a complex large-scale system managed by the public entity Aguas de Barcelona (AGBAR), S.A. This system supplies water to the metropolitan area of Barcelona and it is fed by the Ter and Llobregat rivers using regulated dams with an overall capacity of 600 hm${}^3$. Besides the rivers, some additional underground wells also contribute to an overall inflow of around 7 ${\mathrm{m}}^3/\mathrm{s}$, which becomes potable by four drinking water treatment plants.

The Barcelona DWN can be broken down in two layers. The first (upper) layer consists of a transport network, which connects the water treatment plants with reservoirs distributed across the metropolitan area. The second (lower) layer is the distribution network, which in turn is subdivided in subnetworks that guarantee the water supply from the reservoirs to the consumers.

We have selected this case study to prove that the proposed approach can be applied to a real large-scale problem. In particular, we focus here on the transport network. Hence, the subnetworks within the distribution network are considered as demand sectors that will be characterized by a scheduled pattern and taken as disturbances by the control system.

5.1. Barcelona DWN Description

The Barcelona DWN model [41,43] is depicted in Figure 4 and consists of $n_x=63$ tanks, $n_u=114$ actuators (75 pumps and 39 valves), $n_n=17$ junction nodes and $n_d=88$ sectors of water demand, which are considered as known disturbances. The system can be modeled using flow-based differential-algebraic equations, which interrelate water levels x at tanks, controlled pipe flows u, and demands d. In particular, we have:

• Water tanks differences equations

  $$x_i(k+1)=x_i\left(k\right)+\Delta t\left(\sum _i{q}_{\mathrm{in},i}\left(k\right)-\sum _j{q}_{\mathrm{out},j}\left(k\right)\right),$$

  (14)

  where $x_i$ is the water level in tank i, and $q_{\mathrm{in},i}$ and $q_{\mathrm{out},j}$ are respectively the i-th and j-th inflows and outflows in m${}^3/s$. The aggregation of all differences equations allows us to reformulate the problem as

  $$x_{\mathcal{N}}(k+1)=A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+B_{\mathcal{N}}{u}_{\mathcal{N}}\left(k\right)+B_p{d}_{\mathcal{N}}\left(k\right),$$

  (15)

  with $d_{\mathcal{N}}\in {\mathbb{R}}^{n_d}$ being the disturbances vector, $B_p\in {\mathbb{R}}^{n_x\times n_d}$, and where the rest of variables are defined as introduced in Section 2.
• Mass-balance constraints imposed by the nodes

  $$\sum _i{q}_{in,i}\left(k\right)=\sum _j{q}_{out,j}\left(k\right),$$

  (16)

  with $q_{\mathrm{in},i}$ and $q_{\mathrm{out},j}$ defined as before. Equation (16) can be rewritten in matrix form and considering the known disturbances as

  $$E_u{u}_{\mathcal{N}}\left(k\right)+E_d{d}_{\mathcal{N}}\left(k\right)=0,$$

  (17)

  where $E_u\in {\mathbb{R}}^{n_n\times n_u}$ and $E_d\in R^{n_n\times n_d}$ respectively deal with the flows associated with the control variables and those corresponding to the water demands.
• Bounds on inputs, i.e.,

  $$u_{\min }\le u_{\mathcal{N}}\left(k\right)\le u_{\max },$$

  (18)

  where the values $u_{\min }$ and $u_{\max }$ are the upper and lower limits of the different actuators at the DWN, respectively.
• Bounds on states at tanks, i.e.,

  $$x_{\min }\le x_{\mathcal{N}}\left(k\right)\le x_{\max },$$

  (19)

  being $x_{\min }$ and $x_{\max }$ respectively the minimum and maximum levels at the water tanks.

5.2. Control Variables Parameterization

Equation (17) relates the control variables with the measured disturbances in the system nodes. Based on [61], let us assume that $\mathrm{rank}\left(E_u\right)=\mathrm{rank}\left(E_d\right)=n_n$, with $n_n\le n_d\le n_u$, i.e., some components of $u_{\mathcal{N}}\left(k\right)$ are not longer independent and can be parameterized as a function of the known disturbances. To this end, let us recast (17), as

$$\left[\begin{array}{cc}{E}_u& E_d\end{array}\right]\left[\begin{array}{c}{u}_{\mathcal{N}}\left(k\right)\\ d_{\mathcal{N}}\left(k\right)\end{array}\right]=0.$$

(20)

At this point, we consider a linear transformation $\widehat{P}\in {\mathbb{R}}^{n_u\times n_u}$ that allows us to perform the Gauss-Jordan elimination, i.e.,

$$E_u\widehat{P}=\left[\begin{array}{cc}{I}_{n_n}& M_1\end{array}\right],M_1\in {\mathbb{R}}^{n_n\times n_u-n_n},$$

(21)

which yields

$$\left[\begin{array}{cc}{E}_u& E_d\end{array}\right]P=\left[\begin{array}{ccc}{I}_{n_n}& M_1& M_2\end{array}\right],M_2\in {\mathbb{R}}^{n_n\times n_d},\mathrm{with}P=\left[\begin{array}{cc}\widehat{P}& 0\\ 0& I_{n_n}\end{array}\right].$$

(22)

Now, it is possible to utilize P to reformulate (20), obtaining

$$\left[\begin{array}{cc}{E}_u& E_d\end{array}\right]P{P}^T\left[\begin{array}{c}{u}_{\mathcal{N}}\left(k\right)\\ d_{\mathcal{N}}\left(k\right)\end{array}\right]=0\Rightarrow \left[\begin{array}{cc}{E}_u\widehat{P}& E_d\end{array}\right]\left[\begin{array}{c}{\widehat{P}}^T{u}_{\mathcal{N}}\left(k\right)\\ d_{\mathcal{N}}\left(k\right)\end{array}\right]=0.$$

(23)

Then, in order to separate the effects of the dependent and independent control inputs, we define

$$v_{\mathcal{N}}\left(k\right)={\widehat{P}}^T{u}_{\mathcal{N}}\left(k\right)=\left[\begin{array}{c}{\overline{u}}_{\mathcal{N}}\left(k\right)\\ {\widehat{u}}_{\mathcal{N}}\left(k\right)\end{array}\right],$$

(24)

reaching the following expression equivalent to (20):

$$\left[\begin{array}{ccc}{I}_{n_n}& M_1& M_2\end{array}\right]\left[\begin{array}{c}{\overline{u}}_{\mathcal{N}}\left(k\right)\\ {\widehat{u}}_{\mathcal{N}}\left(k\right)\\ d_{\mathcal{N}}\left(k\right)\end{array}\right]=0\Rightarrow {\overline{u}}_{\mathcal{N}}\left(k\right)=-M_1{\widehat{u}}_{\mathcal{N}}\left(k\right)-M_2{d}_{\mathcal{N}}\left(k\right).$$

(25)

Next, using the definition of $v_{\mathcal{N}}\left(k\right)$

$$v_{\mathcal{N}}\left(k\right)=\left[\begin{array}{c}{\overline{u}}_{\mathcal{N}}\left(k\right)\\ {\widehat{u}}_{\mathcal{N}}\left(k\right)\end{array}\right]=\left[\begin{array}{c}-M_1{\widehat{u}}_{\mathcal{N}}\left(k\right)-M_2{d}_{\mathcal{N}}\left(k\right)\\ {\widehat{u}}_{\mathcal{N}}\left(k\right)\end{array}\right]=\underset{\widetilde{M_1}}{\underbrace{\left[\begin{array}{c}-M_1\\ I_{n_u-n_n}\end{array}\right]}}{\widehat{u}}_{\mathcal{N}}\left(k\right)+\underset{\widetilde{M_2}}{\underbrace{\left[\begin{array}{c}-M_2\\ 0_{n_d}\end{array}\right]}}{d}_{\mathcal{N}}\left(k\right).$$

(26)

Equation (26) enables the parameterization of $u_{\mathcal{N}}\left(k\right)$

$$u_{\mathcal{N}}\left(k\right)=\widehat{P}{v}_{\mathcal{N}}\left(k\right)=\widehat{P}\widetilde{M_1}{\widehat{u}}_{\mathcal{N}}\left(k\right)+\widehat{P}\widetilde{M_2}{d}_{\mathcal{N}}\left(k\right).$$

(27)

Substituting (27) into (15) yields

$$\begin{array}{ccc}\hfill x_{\mathcal{N}}(k+1)& =& A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+B_{\mathcal{N}}{u}_{\mathcal{N}}\left(k\right)+B_p{d}_{\mathcal{N}}\left(k\right)\hfill \\ \hfill & =& A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+B_{\mathcal{N}}\left(\widehat{P}\widetilde{M_1}{\widehat{u}}_{\mathcal{N}}\left(k\right)+\widehat{P}\widetilde{M_2}{d}_{\mathcal{N}}\left(k\right)\right)+B_p{d}_{\mathcal{N}}\left(k\right).\hfill \end{array}$$

(28)

Finally, reorganizing terms, we have

$$\begin{array}{ccc}\hfill x_{\mathcal{N}}(k+1)& =& A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+\underset{{\tilde{B}}_{\mathcal{N}}}{\underbrace{(B_{\mathcal{N}}\widehat{P}\widetilde{M_1})}}{\widehat{u}}_{\mathcal{N}}\left(k\right)+\underset{\widetilde{B_p}}{\underbrace{(B_p+B_{\mathcal{N}}\widehat{P}\widetilde{M_2}}})d_{\mathcal{N}}\left(k\right)\hfill \\ \hfill & =& A_{\mathcal{N}}{x}_{\mathcal{N}}\left(k\right)+{\tilde{B}}_{\mathcal{N}}{\widehat{u}}_{\mathcal{N}}\left(k\right)+\widetilde{B_p}{d}_{\mathcal{N}}\left(k\right).\hfill \end{array}$$

(29)

Summing up, this procedure allows us to reduce the size of the system, while forcing the fulfillment of the node equations in (17) into the aggregated state-space equation. More specifically, the number of control actions is reduced from $n_u=114$ to $n_u-n_d=97$, which in turn reduces the size of matrices $Y_{\mathsf{\Lambda }}$ in LMI (11). Finally, note that terms I and 0 have been used along this procedure to symbolize, respectively, the identity and null matrices of the corresponding dimensions.

5.3. Modular Controller for the Barcelona DWN

The states, control actions and disturbances of the Barcelona DWN are divided according to the eight agent partitioning proposed in [44], which in turn is based on that of [42]. Please note that this partitioning was made assuring that node equations are implicitly satisfied. From the agents viewpoint, the system considered is depicted in Figure 5, where agents are represented by blue circles, connections by green lines, and they are symbolized by Arabic and Roman numerals, respectively.

Since the purpose of this large-scale example is to illustrate the applicability of the proposed method, the modular controller synthesis has been performed considering coalitions of agents, i.e., network topology constraints have not been taken into account. In this way, it is assumed that each pair of agents can communicate, which leads to the simultaneous resolution of $2^8$ cooperation scenarios. Nevertheless, the resulting solution satisfies any situation of communication that can stem from the links depicted in Figure 5, and it is less demanding from a computational viewpoint, because less LMI constraints are imposed on the problem.

Again, since this case study is considered for illustration purposes, we have used the flow-based model of the Barcelona DWN depicted in Figure 4 with a sampling time of 1 h. In this way, each tank becomes an integrator and control actions and disturbances represent how many cubic meters of water are transferred through the corresponding pipe at each time step. Likewise, stage cost ${\ell }_{\mathcal{N}}\left(k\right)$ introduced in (4) is defined using unit matrices of the corresponding size for simplicity, i.e., $Q_{\mathcal{N}}=I_{63}$ and $R_{\mathcal{N}}=I_{97}$. This cost has been extended to account for coordination costs by simply adding a penalty consisting on the number of nonzero elements in the feedback gain employed times a weight ${\gamma }_{\mathrm{c}}={10}^{-5}$ tuned by a trial-and-error procedure. Hence, the modular controller must find a trade-off between cooperation burden and performance. Finally, note that mass-balance constraints have not explicitly been taken into account in the controller design, although the parameterization of the control variables presented in Section 5.2 allows us to guarantee that they hold. As for bounds on tanks and flows, there are tools within the LMI framework to impose them if necessary, e.g., see [62].

Taking into account the aforementioned parameters values, the simulations have been implemented using the Matlab^® toolbox, in a 3.6 GHz Intel^® Octa Core^TM/32 GB RAM computer. The LMI problem executed to obtain the modular controller has around 6000 variables and $2^8=256$ constraints and it is solved in a few hours by the computer. This time is expected to increase nonlinearly as more constraints and variables are added, which means that larger problems might need to resort on additional relaxations as those proposed in Section 4. Nevertheless, since the problem can be solved offline, it can be admissible to have computation times in the order of days.

The controllers considered for assessing the performance of the proposed control method are:

• Modular controller (MOD), obtained using the procedure described in Section 3. Here, the topologies are assessed every three time steps, so that the feedback controller providing the minimum expected cost before the next topology change instant is selected. Please note that system stability is guaranteed despite topology changes due to the existence of a common Lyapunov function.
• Linear-quadratic regulator (LQR), designed for the centralized system.
• Decentralized controller (DEC), which has been obtained following the coalitional approach of [8] for the decentralized communication topology.

The second and third controllers are provided as a means to illustrate the maximum performance that can be expected of a centralized and a decentralized feedback, respectively. Please note that the coordination costs are constant for both LQR and DEC because they cannot change their structure.

The performance of the previously mentioned controllers will be compared using accumulated stage costs, with and without disturbances. In particular, the simulation starts from a random state and the system goes undisturbed during 50 time steps to assess the controllers from a pure regulation viewpoint. After that, the system is fed with disturbances representing the water demands in the different DWN sectors for 450 time instants. Since the dynamics is linear, note that the origin of the linear system can be interpreted as the operation point of the real system.

Accumulated costs are presented in Figure 6, where we can quickly examine the regulation capability of each controller by looking at the costs during the first 50 instants. As expected, the best performance is that of LQR and the worst one is that of DEC, with the performance of the modular controller lying in between. Please note that during this short interval, coordination costs are small and do not affect significantly the sheer control costs of each method. Once demands come into play, the setup is different from the classical regulation problem. As a consequence, the performance of the controllers change, with LQR being the worst one in this part of the simulation. As counterintuitive as this fact can be, LQR is only optimal in very specific conditions that do not hold during this period. Also, coordination efforts penalize LQR, which at around 300 instants is outperformed by MOD. As observed by the end of the simulation, also DEC eventually outperforms LQR due to its lower coordination costs. Nevertheless, its performance never increases enough to improve the results of MOD.

The evolution of some of the system water levels around a fixed operation point is shown in Figure 7. Notice that the oscillating behavior observed is caused by the periodicity of the water demands of the network, which implies the water volumes into the tanks should follow the corresponding oscillatory profile towards reaching the control objective.

Finally, Figure 8 shows the evolution of the topology used by the modular controller and the corresponding coordination efforts. Please note that the blue line indicates the topology being used at each time step (left y-axis). In particular, the topology ranges between 0 and 255 and when converted into binary digits provides us with the agents that are cooperating, i.e., if the i-th binary digit is ‘1’, then the i-th agent is communicating so as to coordinate its actions with other agents. Likewise, the discontinuous red line indicates the evolution of coordination costs (right y-axis) and clearly shows how the control architecture reacts to the disturbances the overall system receives. As can be seen, MOD can adjust its structure and achieves a better trade-off between performance and coordination efforts, which is key to understand its superior efficiency. Notice that the initial random state is farther away than the net effect of the disturbances on the state. Hence, higher coordination costs are incurred to steer the system towards the origin at the very beginning. As the system reaches the origin, coordination costs reach the minimum (that of a decentralized control architecture). After that, coordination efforts are regularly adjusted as a response to the pattern of the demand.

As the obtained results show, the implementation of the proposed control strategy not only provides with an approach able to reach the control objectives while satisfying the physical and operational constraints of the network, but also takes into account two key factors that make it suitable to be implemented in these large-scale water systems:

• The former is related to the reliable modularity the control strategy confers since the overall system gains certain upper-level of robustness against fault events that might occur. Notice that the coordination and noncentralized features of the proposed approach make that, once a problem takes place, the system keeps working by isolating the affected part while the rest is self-adjusted to provide water to demand sectors.
• On the other hand, the latter factor, also operationally related to the former, relies on the fact that the approach acts as a fast decision-maker given its offline design and low online computational burden. This feature allows the system operator to promptly react facing possible problems (caused by a fault event), avoiding scenarios of water supply lacking.

6. Conclusions

A modular feedback design method has been presented to generate a family of feedback controllers suitable for different topologies in a networked control system. The proposed method leverages linear matrix inequalities (LMIs) to attain the modular structure adding a certain upper-level of robustness against fault events, while guaranteeing stability despite topology switchings. In particular, the method starts from a linear model of the system and requires to solve an optimization problem subject to as many LMI constraints as possible cooperation scenarios can be defined using the available communication resources. Even when computational burden does not increase linearly with the number of topologies considered, it is possible to apply the method in large-scale systems as the one chosen as case study: the Barcelona drinking water network (DWN), where in a scenario with disturbances and communication costs the modular control manages to stabilize the system outperforming other approaches in terms of accumulated costs.

Moreover, the resulting controller has significant advantages that can be useful in networking control applications where packet losses might occur, for these events can be interpreted as topology switchings. Also, plug-and-play and coalitional control strategies can benefit from this method because of the much simpler implementation of the modular controller. In any case, the results of the proposed scheme are sensitive to the tuning parameters employed, and, particularly, to the penalty for cooperation efforts. The controller might tend to either decentralized or centralized configurations if this penalty is not properly adjusted, thus deteriorating its performance, especially when compared with specifically designed feedback gains for every topology. Finally, future work should include the development of distributed synthesis techniques for this family of feedback controllers.