Next Article in Journal
Gottlieb Polynomials and Their q-Extensions
Next Article in Special Issue
Solution of Moore–Gibson–Thompson Equation of an Unbounded Medium with a Cylindrical Hole
Previous Article in Journal
Open Markov Type Population Models: From Discrete to Continuous Time
Previous Article in Special Issue
Integrated Structure-Control Design of a Bipedal Robot Based on Passive Dynamic Walking
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Communication

Efficient Algorithm for the Computation of the Solution to a Sparse Matrix Equation in Distributed Control Theory

Institute for Systems and Robotics, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(13), 1497; https://doi.org/10.3390/math9131497
Submission received: 29 May 2021 / Revised: 23 June 2021 / Accepted: 23 June 2021 / Published: 25 June 2021
(This article belongs to the Special Issue Numerical Analysis and Scientific Computing)

Abstract

:
In this short communication, an algorithm for efficiently solving a sparse matrix equation, which arises frequently in the field of distributed control and estimation theory, is proposed. The efficient algorithm stems from the fact that the sparse equation at hand can be reduced to a system of linear equations. The proposed algorithm is shown to require significantly fewer floating point operations than the state-of-the-art solution. The proposed solution is applied to a real-life example, which models a wide range of industrial processes. The experimental results show that the solution put forward allows for a significant increase in efficiency in relation to the state-of-the-art solution. The significant increase in efficiency of the presented algorithm allows for a valuable widening of the applications of distributed estimation and control.

1. Notation

The n × o null matrix and the null vector of dimension n are denoted by 0 n × o and 0 n , respectively. The i-th component of a vector v is denoted by [ v ] i , and the entry ( i , j ) of a matrix A is denoted by [ A ] i , j . The vectorization of a matrix A , denoted herein by vec ( A ) , returns a vector composed of the concatenated columns of A . The Kronecker product of two matrices A and B is denoted by A B . The cardinality of a set χ is denoted by | χ | .

2. Introduction

Sparse matrix equations arise in various real-life problems in a broad range of engineering fields. Thus, given their impact, they have been extensively studied for a long time [1]. In recent developments in the field of distributed estimation and control theory, a sparse matrix equation has arisen, for which, to the best of the knowledge of the authors, an efficient algorithm for its solution has not yet been proposed.
Consider the matrix equation
A X B C i , j = 0 , ( i , j ) χ X i , j = 0 , ( i , j ) χ ,
where A R n × n , B R o × o , and C R n × o are known and X R n × o is the unknown. Let E R n × o represent a known sparsity pattern, and define the set χ of integer pairs of the form ( i , j ) to index the nonzero entries of E as
( i , j ) χ , E i , j 0 ( i , j ) χ , otherwise , i = 1 , , n , j = 1 , , o .
Equation (1) arises frequently in the field of distributed control and estimation theory. In fact, it is the backbone of the state-of-the-art methods for distributed filter design [2] and controller synthesis [3] for large-scale dynamical systems. However, the solution to (1), presented in these papers, relies on an inefficient closed-form matrix computation. Define a matrix Z , such that the vector Z vec ( E ) contains all the nonzero elements of E . The solution to (1) presented in the literature is given by
vec ( X ) = Z T Z B A Z T 1 Z vec ( C ) ,
as proven in ([2], Appendix B). Note that (3) requires the multiplication of a matrix of dimension | χ | × n o by other dimension n o × n o , which requires O ( | χ | ( n o ) 2 ) floating point operations. The number of practical applications, where it is convenient to model the interconnection of complex systems as a whole network, has been steadily increasing [4]. They have emerged in a broad range of engineering fields such as irrigation networks [5,6], power systems networks [7,8], traffic networks [9,10,11,12], and industrial processes [13]. Nevertheless, given the substantial increase in memory and processing power needed to control such systems, the classical control solutions cannot be implemented in practice. That is why, over the past decades, an effort has been made towards the development of distributed and computationally efficient algorithms to handle such systems. In this short communication, an algorithm to obtain an efficient and exact solution to (1) is proposed, which significantly increases the applicability of distributed methods that address the control and estimation problems for large-scale dynamical systems.

3. Materials and Methods

The following result is the basis for the algorithm proposed in this short communication.
Theorem 1.
The system of Equation (1) can be reduced to a system of linear equations of dimension | χ | ,
S x = r ,
where S R | χ | × | χ | , r R | χ | , and x R | χ | is the vector of the ordered nonzero entries of vec ( X ) .
Proof. 
For every integer pair ( i , j ) χ , the first equation of (1) is given by
k = 1 n l = 1 o A i , k X k , l B l , j C i , j = 0 .
Given the second equation of (1), one can rewrite (5) as
( k , l ) χ A i , k X k , l B l , j C i , j = 0 .
Denote by x R | χ | the vector of the ordered nonzero entries of vec ( X ) , and let ϕ = [ ϕ 1 ϕ 2 ] T : N 1 N 2 denote the resulting one-to-one map between the entries of x and the nonzero entries of X . Defining p : = ϕ 1 ( i , j ) and q : = ϕ 1 ( k , l ) , it follows that satisfying (6) for every integer pair ( i , j ) χ is equivalent to satisfying S x = r , where S R | χ | × | χ | and r R | χ | are defined such that
S p , q = A ϕ 1 ( p ) , ϕ 1 ( q ) B ϕ 2 ( q ) , ϕ 2 ( p ) ,
p = 1 , , | χ | , q = 1 , , | χ | ; and
r p = C ϕ ( p ) , p = 1 , , | χ | ;
respectively. It then follows that the solution to (1) is
X i , j = x ϕ 1 ( i , j ) , ( i , j ) χ 0 , ( i , j ) χ ,
where x is the solution of (4). □
Corollary 1.
It follows immediately from Theorem 1 that (1) admits: (i) one and only one solution if S is full rank; (ii) infinitely many solutions if (4) is consistent and S is not full rank; and (iii) no solution if (4) is inconsistent.
Proposition 1.
Algorithm 1 can be used to obtain the solution to (1), with O ( | χ | 3 ) floating point operations.
Proof. 
Algorithm 1 computes S R | χ | × | χ | and r R | χ | according to (7) and (8), respectively. It, then, follows from Theorem 1 that (1) can be reduced to the solution to S x = r . The LU decomposition of S can be carried out with O ( | χ | 3 ) floating point operations, using Gaussian elimination. The type of solution is then assessed according to Corollary 1. If the solution is unique, it follows from Theorem 1 that it is given by (9), which is implemented in Algorithm 1. □
Algorithm 1 Algorithm to efficiently compute the exact solution to (1)
Mathematics 09 01497 i001
Remark 1.
It is important to point out that each of the three main steps of Algorithm 1 can be carried out in parallel. The first and last steps for loops of Algorithm 1 (which are commented on) can be carried out in parallel, for each j. In this case, the variable p has to be initialized, for each j, with the number of nonzero entries of E in the first j 1 columns. Moreover, the LU decomposition can also be carried out in parallel [14,15].
Remark 2.
There are asymptotically faster algorithms for computing the LU decomposition of a square matrix. For instance, the Strassen algorithm [16] can perform the LU decomposition of S in O ( | χ | log 2 7 ) floating point operations. Moreover, the Coppersmith–Winograd algorithm [17] requires only O ( | χ | 2.38 ) floating point operations. However, an improvement in relation to the Strassen algorithm can only be noticed for large values of | χ | . For more details, see ([18], Section 1.4).
Remark 3.
In the field of distributed estimation and control theory, | χ | is the sum of the number of sensor output signals available to each node of the network, and n is the order of the state–space dynamics of the whole system network. Given that, for most applications, each node has a similar order and a similar number of sensor output signals available, then | χ | c n , where c N is a constant. The system detailed in Section 4 illustrates this claim. It, thus, follows that the ratio of the number of floating points operations required by the proposed algorithm and the number required by the state-of-the-art methods is O ( o 2 ) .

4. Results

In this section, the solution to the sparse matrix equation (1) is used to design a distributed filter for a large-scale network of N tanks. The computational efficiency of the synthesis with the proposed solution is compared to the state-of-the-art synthesis procedure. This network is widely studied [3,19,20] because it is analogous to a broad range of industrial processes. Consider N interconnected tanks, as shown in Figure 1, where N is an even integer. The water level of tank i is denoted by h i . The network is actuated by N / 2 pumps, which are controlled by the lower tanks, whose inputs are denoted by u i for i = 1 , , N / 2 , in accordance with the schematic. Each pump is connected to a three-way valve that regulates the fraction of the flow, held constant, that goes to each of the tanks supplied by the pump. Each tank has a sensor, which measures its water level. Making use of mass balances and Bernoulli’s law, the system dynamics, in the absence of noise, are given by
A i h ˙ i ( t ) = a i 2 g h i ( t ) + a N 2 + i 2 g h N 2 + i ( t ) + γ i k i u i ( t ) , i = 1 , , N / 2 A i h ˙ i ( t ) = a i 2 g h i ( t ) + ( 1 γ i N 2 1 ) k i N 2 1 u i N 2 1 ( t ) , i = N 2 + 2 , , N A N 2 + 1 h ˙ N 2 + 1 ( t ) = a N 2 + 1 2 g h N 2 + 1 ( t ) + ( 1 γ N 2 ) k N 2 u N 2 ( t ) , ,
where A i and a i are the cross sections of tank i and of its outlet hole, respectively; the constant γ i represents the fraction of the flow that passes through the valve i to the lower tanks; k i is the constant of proportionality between the mass flow and the input of pump i; and g denotes the acceleration of gravity.
The values of the constant parameters of the network are presented in Table 1. The dynamics of the network are nonlinear, thus the system is linearized about its operational point, which is selected to be h ( i ) = 15 cm, i = 1 , , N . The distributed filter is synthesized with the finite-horizon method ([2], Section 5), which is shown to have better state estimation performance in relation to other distributed filter synthesis methods. Each iteration k of this method computes a new gain, K ( k ) R n × o , solving
Λ ( k + 1 ) K ( k ) S ( k ) Λ ( k + 1 ) P ( k | k 1 ) C T i , j = 0 , ( i , j ) χ K ( k ) i , j = 0 , ( i , j ) χ ,
where Λ ( k + 1 ) R n × n , S ( k ) R o × o , and P ( k | k 1 ) R n × n are defined in [2], with n = o = N . Moreover, C = I N is the output matrix of the network and, given the configuration of the network, the sparsity pattern that defines the set χ is given by E = I N . Note that (11) has the same form as (1), for which a novel efficient solution is proposed in this short communication.
The goal is to compare the computational efficiency of the gain synthesis using the solution to (11) proposed in this short communication and compare it to the state-of-the-art solution (3). In that sense, the gain synthesis is performed for both alternatives, for a range of N. Figure 2 depicts the elapsed wall-clock time for the distributed filter synthesis, resulting from the average of three simulations of a MATLAB implementation on a single thread of a server with 24GB RAM and 24 Intel Xeon hexa-core CPUs at 2.40GHz. A MATLAB implementation of the distributed filter synthesis, of the application to the N tanks network, and of Algorithm 1, can be found in the DECENTER toolbox available at https://decenter2021.github.io (accessed on 23 June 2021). Figure 2 also depicts the linear regression in the logarithmic scale of the asymptotic evolution, which is considered to be achieved for N > 40 . First, if the solution of (11) is the bottleneck of the gain synthesis procedure, one would expect that the synthesis wall-clock time grows asymptotically with O ( N 5 ) for the state-of-the-art solution, and O ( N 3 ) for the novel solution proposed in this short communication. It is possible to conclude from Figure 2 that this is, in fact, the case. The experimental complexity is slightly lower than what was projected, since MATLAB makes use of efficient LU decomposition algorithms like the ones pointed out in Remark 2. Second, it is also possible to confirm that the use of the solution proposed in this short communication leads to a ratio of wall-clock synthesis time between the proposed solution and the state-of-the-art solution that evolves with O ( N 2.035 ) . This improvement significantly widens the applicability of various distributed methods to large-scale systems.

5. Conclusions

In this short communication, an algorithm for efficiently solving the sparse matrix equation (1), which arises frequently in the field of distributed control and estimation theory, is proposed. The state-of-the-art methods for distributed filter design and controller synthesis rely on a very inefficient computation of the solution to (1), requiring O ( | χ | ( n o ) 2 ) floating point operations. In this letter, an algorithm is proposed to compute the solution to this equation with O ( | χ | 3 ) floating point operations. For distributed estimation and control theory applications, the ratio of the number of floating points operations required by the proposed algorithm and the number required by the state-of-the-art methods is O ( o 2 ) . Furthermore, the algorithm proposed in this short communication was applied to an example representative of real-life large-scale industrial processes, and was compared to the state-of-the-art solution. It is shown that the solution put forward allows for a decrease of the computational complexity of the gain synthesis from O ( N 4.805 ) , which is obtained for the state-of-the-art solution, to O ( N 2.770 ) , where N is the dimension of the network. The significantly greater efficiency achieved with the proposed algorithm allows for the state-of-the-art distributed control and estimation algorithms to be implemented in an online framework and for large-scale dynamical systems, which would otherwise be unfeasible.

Author Contributions

Conceptualization, L.P.; methodology, L.P. and P.B.; software, L.P.; validation, L.P. and P.B.; formal analysis, L.P.; writing—original draft preparation, L.P.; writing—review and editing, L.P. and P.B; visualization, L.P.; supervision, P.B.; project administration, P.B.; funding acquisition, P.B. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Fundação para a Ciência e a Tecnologia (FCT) through LARSyS-FCT Project UIDB/50009/2020 and through the FCT project DECENTER [LISBOA-01-0145-FEDER-029605], funded by the Programa Operacional Regional de Lisboa 2020 and PIDDAC programs.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Duff, I.S. A survey of sparse matrix research. Proc. IEEE 1977, 65, 500–535. [Google Scholar] [CrossRef]
  2. Viegas, D.; Batista, P.; Oliveira, P.; Silvestre, C. Discrete-time distributed Kalman filter design for formations of autonomous vehicles. Control Eng. Pract. 2018, 75, 55–68. [Google Scholar] [CrossRef]
  3. Viegas, D.; Batista, P.; Oliveira, P.; Silvestre, C. Distributed controller design and performance optimization for discrete-time linear systems. Optim. Control. Appl. Methods 2020, 1–18. [Google Scholar] [CrossRef]
  4. Šiljak, D.D.; Zečević, A. Control of large-scale systems: Beyond decentralized feedback. Annu. Rev. Control 2005, 29, 169–179. [Google Scholar] [CrossRef]
  5. Conde, G.; Quijano, N.; Ocampo-Martinez, C. Modeling and control in open-channel irrigation systems: A review. Annu. Rev. Control 2021, 51, 153–171. [Google Scholar] [CrossRef]
  6. Cantoni, M.; Weyer, E.; Li, Y.; Ooi, S.K.; Mareels, I.; Ryan, M. Control of large-scale irrigation networks. Proc. IEEE 2007, 95, 75–91. [Google Scholar] [CrossRef]
  7. Chen, C.; Wang, J.; Kishore, S. A distributed direct load control approach for large-scale residential demand response. IEEE Trans. Power Syst. 2014, 29, 2219–2228. [Google Scholar] [CrossRef]
  8. Bumiller, G.; Lampe, L.; Hrasnica, H. Power line communication networks for large-scale control and automation systems. IEEE Commun. Mag. 2010, 48, 106–113. [Google Scholar] [CrossRef]
  9. Tan, T.; Bao, F.; Deng, Y.; Jin, A.; Dai, Q.; Wang, J. Cooperative deep reinforcement learning for large-scale traffic grid signal control. IEEE Trans. Cybern. 2019, 50, 2687–2700. [Google Scholar] [CrossRef]
  10. Keyvan-Ekbatani, M.; Yildirimoglu, M.; Geroliminis, N.; Papageorgiou, M. Multiple concentric gating traffic control in large-scale urban networks. IEEE Trans. Intell. Transp. Syst. 2015, 16, 2141–2154. [Google Scholar] [CrossRef]
  11. Carlson, R.C.; Papamichail, I.; Papageorgiou, M.; Messmer, A. Optimal mainstream traffic flow control of large-scale motorway networks. Transp. Res. Part C Emerg. Technol. 2010, 18, 193–212. [Google Scholar] [CrossRef]
  12. Li, Y.; Yang, L.; Yang, G. Network-based coordinated motion control of large-scale transportation vehicles. IEEE/ASME Trans. Mechatron. 2007, 12, 208–215. [Google Scholar] [CrossRef]
  13. Vadigepalli, R.; Doyle Iii, F.J. Structural analysis of large-scale systems for distributed state estimation and control applications. Control Eng. Pract. 2003, 11, 895–905. [Google Scholar] [CrossRef]
  14. Buoni, J.J.; Farrell, P.A.; Ruttan, A. Algorithms for LU decomposition on a shared memory multiprocessor. Parallel Comput. 1993, 19, 925–937. [Google Scholar] [CrossRef]
  15. Liu, Z.; Cheung, D. Efficient parallel algorithm for dense matrix LU decomposition with pivoting on hypercubes. Comput. Math. Appl. 1997, 33, 39–50. [Google Scholar] [CrossRef] [Green Version]
  16. Strassen, V. Gaussian elimination is not optimal. Numer. Math. 1969, 13, 354–356. [Google Scholar] [CrossRef]
  17. Coppersmith, D.; Winograd, S. Matrix multiplication via arithmetic progressions. In Proceedings of the Nineteenth Annual ACM Symposium on Theory of Computing, New York, NY, USA, 25–27 May 1987; pp. 1–6. [Google Scholar]
  18. Pan, V. Complexity of algorithms for linear systems of equations. In Computer Algorithms for Solving Linear Algebraic Equations; Springer: Berlin/Heidelberg, Germany, 1991; pp. 27–56. [Google Scholar]
  19. Johansson, K.H. The quadruple-tank process: A multivariable laboratory process with an adjustable zero. IEEE Trans. Control Syst. Technol. 2000. [Google Scholar] [CrossRef] [Green Version]
  20. Casavola, A.; Garone, E.; Tedesco, F. A distributed multi-agent command governor strategy for the coordination of networked interconnected systems. IEEE Trans. Autom. Control 2014, 59, 2099–2112. [Google Scholar] [CrossRef]
Figure 1. Schematic of the N tanks network.
Figure 1. Schematic of the N tanks network.
Mathematics 09 01497 g001
Figure 2. Synthesis wall-clock time of a distributed filter for the N tanks network, as a function of the dimension of the network.
Figure 2. Synthesis wall-clock time of a distributed filter for the N tanks network, as a function of the dimension of the network.
Mathematics 09 01497 g002
Table 1. Values of the physical constants of the N tanks network.
Table 1. Values of the physical constants of the N tanks network.
ConstantValue
A i , i odd 28 cm 2
A i , i even 32 cm 2
a i , i odd N / 2 0.071 cm 2
a i , i even N / 2 0.057 cm 2
a i , i > N / 2 0.040 cm 2
g 981 cm s 2
k i 3.33 cm 3 s 1 V 1
γ i , i odd 0.7
γ i , i even 0.6
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Pedroso, L.; Batista, P. Efficient Algorithm for the Computation of the Solution to a Sparse Matrix Equation in Distributed Control Theory. Mathematics 2021, 9, 1497. https://doi.org/10.3390/math9131497

AMA Style

Pedroso L, Batista P. Efficient Algorithm for the Computation of the Solution to a Sparse Matrix Equation in Distributed Control Theory. Mathematics. 2021; 9(13):1497. https://doi.org/10.3390/math9131497

Chicago/Turabian Style

Pedroso, Leonardo, and Pedro Batista. 2021. "Efficient Algorithm for the Computation of the Solution to a Sparse Matrix Equation in Distributed Control Theory" Mathematics 9, no. 13: 1497. https://doi.org/10.3390/math9131497

APA Style

Pedroso, L., & Batista, P. (2021). Efficient Algorithm for the Computation of the Solution to a Sparse Matrix Equation in Distributed Control Theory. Mathematics, 9(13), 1497. https://doi.org/10.3390/math9131497

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop