block = [draw, rectangle,
minimum height=3em, minimum width=4em]
sum = [draw, circle, node distance=1cm]
input = [coordinate]
output = [coordinate]
pinstyle = [pin edge=to-,thin,black]

Consensus-based Distributed Optimization for Multi-agent Systems over Multiplex Networks

Christian D. Rodríguez-Camargo,
Andrés F. Urquijo-Rodríguez,
and Eduardo Mojica-Nava
C. D. Rodríguez-Camargo is with The Atomic, Molecular, Optical and Positron Physics (AMOPP) group of the Department of Physics and Astronomy from the University College London and Programa de Investigación sobre Adquisición y Análisis de Señales (PAAS-UN) from Universidad Nacional de Colombia. E-mail<EMAIL_ADDRESS>
A. F. Urquijo-Rodríguez is with the Industrial Engineering Program and the Centro de Estudios Industriales y Logísticos para la Productividad (CEIL, MD) both from the Corporación Universitaria Minuto de Dios, and Grupo de Superconductividad y Nanotecnología and Grupo de Óptica e Información Cuántica both from the Department of Physics of Universidad Nacional de Colombia.

E. Mojica-Nava is with the Department of Electrical and Electronics Engineering and Programa de Investigación sobre Adquisición y Análisis de Señales (PAAS-UN) from Universidad Nacional de Colombia, Bogotá, Colombia.
E-mail<EMAIL_ADDRESS>Manuscript received MM DD, YYYY; revised MM DD, YYYY.

IEEE Trans. on Network Science and Engineering, Vol. XX, No. X, MM YYYY
Shell et al.: Bare Demo of IEEEtran.cls for Computer Society Journals

Multilayer networks refer to systems composed of several interacting networks. This approach provides a generalized framework for exploring real-world and engineering systems in more depth and with greater accuracy than the traditional single-layer network scheme. Despite the wide set of results in distributed optimization in single-layer networks, there is a lack of developments in multilayer systems. In this paper, we propose two algorithms for distributed optimization problems using the supra-Laplacian matrix of a multiplex network and its diffusion dynamics. The algorithms are a distributed primal-dual saddle-point algorithm and a distributed gradient descent algorithm. We establish a set of saddle-point equations based on the relation between consensus and diffusion dynamics via the graph Laplacian and impose an analogous restriction to the traditional Laplacian in our multiplex distributed optimization problem to obtain the multiplex supra-Laplacian matrix. We extend the distributed gradient descent algorithm for multiplex networks using this supra-Laplacian matrix. We propose several theoretical results to analyze the convergence of both algorithms, demonstrating that they converge to the unique optimizer. We validate the proposed algorithms using numerical examples. We also explore the existence of a phase transition in the consensus time due to the interlayer diffusion constant in the numerical examples. Finally, a coordinated dispatch for interdependent infrastructure networks (energy-gas) is presented to illustrate the application of the proposed framework to real engineering problems.

Distributed optimization, multiplex networks, saddle-point flow, diffusion.

§ INTRODUCTION

During the last two decades, the single-layer network representation of complex systems has proven to be a valuable tool for revealing the relationships between topological properties of different networked systems and their dynamics [1, 2]. The quantitative study of complex networks has generated a wide set of interesting results in the deeper understanding of complex systems throughout the biological, social, engineering, and physical sciences [3, 4].

However, as network science advances and more data is gathered from real-world networked systems, there is a growing need to understand real and engineered systems as networks of networks rather than isolated networks [5]. These systems are commonly referred to as multilayer networks [6].

A multilayer network refers to systems composed of several interacting networks [5, 6, 7]. The first introduction of the multilayer networks has been in the context of social science to characterize the different types of social relationships which exist between the nodes of a social network [7]. Currently, multilayer networks are investigated in many fields, including neuroscience [8], molecular biology [9, 10], ecology [11], transportation [12], and power networks [13]. Despite there are significant advances in phenomena such as percolation [14], cascades failures [15], diffusion dynamics [16], resilient consensus [17], adaptive control [18], epidemic spreading [19], there has been a lack of progress in generalizing single-layer network distributed optimization problems to multilayer networks.

In traditional distributed optimization problems, the objective function is typically a sum of convex functions, each of which is only available to individual nodes [20]. The goal of each node is to obtain the optimal solution for the global convex function through information exchange with neighboring nodes and local computation [21, 22]. Several algorithms have been proposed to solve this distributed optimization problem, from discrete-time gradient descent [21], [22] to recent continuous-time approaches [23], [24]. From a continuous-time perspective, the methods based on saddle-point pint dynamics [25], [26], gradient descent in one variable and gradient ascent in the multipliers, have emerged as an alternative to view optimization algorithms as dynamical systems [27], bringing up new possibilities to deal with complex real-world systems such as complex networks as it is discussed in the recently introduced concept of feedback-based optimization [28], [29].

In this paper, we extend control-based and consensus-based distributed optimization algorithms for single-layer networks to multiplex networks, which generalizes the relationship between consensus and diffusion dynamics in single-layer networks. A multiplex network is a collection of graphs comprised of $M$ distinct layers, where the same set of nodes $N$ can be connected to each other through links belonging to $M$ different classes or types [5]. Multiplex networks are a particular case of multilayer networks, where each node in each layer has a one-to-one correspondence. From another perspective, multiplex networks are composed of the same set of nodes connected by links that indicate different types of interactions [30]. The increasing interest in multiplex networks has highlighted the significant impact of their structure on the behavior of dynamic processes. Therefore, it is crucial to study the properties and dynamics of multiplex networks to comprehend the complex relationships and interdependencies that arise in the real world. From an engineering perspective, it is anticipated that multiplex networks will be utilized to model and manage energy in multienergy systems. A multienergy system is an integrated energy system that employs multiple sources of energy such as electricity, heat, and gas in a coordinated and optimized manner to provide efficient, reliable, and sustainable energy services [31]. The objective is to maximize the use of renewable and low-carbon energy sources, as well as to enhance energy efficiency and reduce greenhouse gas emissions [32]. Multienergy systems are becoming increasingly significant as a means of satisfying the growing demand for energy while decreasing the environmental impact of energy production and consumption [33]. A coordinated dispatch for interdependent infrastructure Networks (energy-gas) is presented to illustrate the application of the proposed framework to real engineering problems.

By recognizing the relationship and mathematical similarity between the consensus equation in multi-agent systems, which is utilized to solve distributed optimization problems, and the diffusion equation in statistical mechanics, it has been possible to propose connections between diffusion processes in multiplex networks and consensus-based distributed optimization algorithms. We expand upon these definitions by utilizing the supra-Laplacian matrix of a multiplex network. Various authors, such as [34, 35], among others, have studied the direct and inverse problems of the $\omega$-Laplace equation on graphs, which can be interpreted as a diffusion equation on electric networks and represents a significant problem in the field of network analysis [36].
In our proposed approach, we derive a distributed primal-dual saddle point algorithm for convex optimization in multiplex networks and a distributed gradient descent algorithm for multiplex networks using the relationship between the diffusion equation and consensus equation for multiplex networks. Through numerical calculations, we show the convergence of each node in each layer to the optimal value of the convex objective function. We demonstrate the accuracy of both algorithms. Notably, only knowledge of the neighborhood is required to achieve the local and global optimal value in multiplex networks. Moreover, we investigate the existence of critical phenomena in the consensus time by manipulating the diffusion constants.

This paper is organized as follows: In Section 2, we present the preliminaries of multiplex networks, the diffusion processes, and the control approach of consensus and diffusion dynamics. In Section 3, we establish the problem statement of distributed optimization in multiplex networks. In Section 4, we present the main result of this paper and a convergence analysis of our theorem related to the distributed primal-dual saddle-point algorithm for multiplex networks. In Section 5 we present the generalized gradient descent algorithm. In Section 6 numerical examples are presented with an additional real-life example in energy management systems. In addition, we evaluate the existence of critical phenomena in the function of the diffusion parameters. The conclusions are presented in Section 7

§ PRELIMINARIES

§.§ Consensus and Diffusion Dynamics From a Control Perspective

By a graph $G = G(V, E)$ we mean a finite set $V$ of vertices with a set $E$ of two-element subsets of $V$ (whose elements are called edges). We denote either $v \in V$ or $v \in G$ the fact that $v$ is a vertex in $G$. A graph $G$ is said to be simple if it has neither multiple edges nor loops, and $G$ is said to be connected if for every pair of vertices $v_r$ and $v_p$ there exists a sequence of vertices $v_r = v_{0} \sim v_{1} \sim v_{2} \sim \cdots \sim v_{n-1} \sim v_{n} = v_p$ such that $v_{j-1}$ and $v_{j}$ are connected by an edge (termed adjacent) for $j = 1, 2, ..., n$, where $v_i \sim v_j$ means that two vertices $v_i$ and $v_j$ are connected (adjacent) by an edge in $E$.

A weighted (undirected) graph is a graph associated with a weight function $\omega: V \times V \rightarrow \mathbb{R}^{+}$ satisfying:

(i)$\omega (v_i,v_j)=\omega (v_j,v_i)$, $v_i, v_j \in V$.

(ii) $\omega (v_i,v_j)=0$ if and only if $\{ v_i, v_j \}\not\in V$.

Here, $\{ v_i, v_j\}$ denotes the edge connecting the vertices $v_i$ and $v_j$.

The degree $d_{\omega}(v_i)$ of a vertex $v_i$ is defined to be

\begin{equation*}
d_{\omega}(v_i):=\sum _{j\in \mathcal{N}_i}\omega (v_i,v_j),
\end{equation*}

where $\mathcal{N}_i$ is the neighborhood of the node $v_i$. A function on a graph is understood to be a function defined just on the set of vertices.

For a measure related with the weighted degree of each vertex, the $\omega$-Laplacian operator $\Delta _{\omega}$ can be represented as a matrix defined by

\begin{equation*}
\Delta _{\omega}(v_i,v_j)=
\begin{cases}
1-\frac{\omega _{v_i,v_i}}{d_{\omega}(v_i)}, \quad \text{if}\, v_i=v_j\, \text{and}\, d_{\omega}(v_i)\neq 0,\\
-\frac{\omega (v_i,v_j)}{d_{\omega}(v_i)}, \quad \text{if}\, v_i\sim v_j, \\
0,\quad \text{otherwise}.
\end{cases}
\end{equation*}

The matrix representation of $\Delta _{\omega}$ is equivalent to the Laplacian matrix $L$ given by

\begin{equation*}
L=D^{1/2}\Delta _{\omega}D^{-1/2},
\label{laplacian001}
\end{equation*}

where $D$ is a diagonal matrix with entries $D(v_i,v_j)=d_{\omega}(v_i)$. $L$ is a non-negative definite symmetric matrix, and its eigenvalues are given by $0<\lambda _{0}\leq \lambda _{1}\leq \lambda _{2}\leq \cdots \leq \lambda _{N-1}$,
and corresponding eigenfunctions $\varphi _{0},\varphi _{1}, \varphi _{2},...,\varphi _{N-1}$,
which are orthonormal in the sense that for each pair of distinct $i$ and $j$: $\sum _{v_i\in V}\varphi_{i}(v_i)\varphi_{j}(v_i)=\delta _{ij}$. It is easy to show that $\lambda _{0}=0$, $\lambda _{1}>0$.

Now, we can introduce the traditional consensus equation for single-layer networks of multiagent systems using the Laplacian matrix. First, consider the dynamic for each node as

\begin{equation}\label{consenagent}
\dot{x}_{i}(t)= -\sum _{j \in \mathcal{N}_{i}}\omega _{ij}(t)(x_{i}(t)-x_{j}(t)),
\end{equation}

where $x_{i}(t)$ is the information state of an agent in node $i$, $\mathcal{N}_{i}$ is the set of agents whose information is available to an agent in node $i$, and $\omega _{ij}(t)$ denotes a positive time-varying weighting factor. Equation <ref> can be presented in compact form as a diffusion equation

\begin{equation} \label{consen}
\dot{x}=-Lx,
\end{equation}

where $L$ is the graph Laplacian.  Using the consensus equation (<ref>), several continuous-time saddle point dynamics have proposed to solve a distributed optimization problem in multiagent systems [25, 26, 23, 24]

\begin{equation}\label{opprob1}
\min _{x \in \mathbb{R}^{n}}\tilde{f}(x)=\sum _{i=1}^{n}f_{i}(x)
\end{equation}

subject to

\begin{equation}\label{cons1}
Lx = 0,
\end{equation}

where $f(x)$ is the global objective function to be minimized subject to the communication constraints between agents represented by the graph Laplacian. The corresponding saddle-point dynamics for optimization problem <ref>-<ref>

\begin{equation*}
\dot{x} = -L x - \frac{\partial \tilde{f}(x)}{\partial x}-L\lambda,
\end{equation*}

\begin{equation*}
\dot{\lambda}= Lx,
\end{equation*}

where $\lambda$ is the vector of Lagrange multipliers [37]. From a control point of view, the saddle-point dynamics have also been studied. Consider the distributed gradient flow with Laplacian coupling

\begin{equation*}\label{ngf001}
K_{1}\dot{y}=-\rho Ly -\frac{\partial f(y)}{\partial y},
\end{equation*}

which comes from the regular gradient flow of the optimization problem within the soft penalty method and control gain $K_1$:

\begin{equation}\label{opprob2}
\min _{y\in \mathbb{R}^{n}}\sum _{i=1}^{n}f_{i}(y_{i})+\frac{\rho}{2}y^{T}Ly,
\end{equation}

where the soft penalty $\frac{\rho}{2}y^{T}Ly$ is tasked with enforcing a Laplacian consensus-constraint.

Several references have suggested the use of a control approach to interpret and improve diffusion and consensus dynamics in situations where control elements can be employed [28, 27, 23] By leveraging this approach, new distributed optimization problems can be obtained. For instance, in [23], the authors deal with the situation where the objective of each node is to find the solution of the following optimization problem

\begin{equation}
\min \sum _{i=1}^{n}f_{i}(y);\quad \text{st.}\, y\in \mathbb{R}^{n},
\label{problemopt001}
\end{equation}

where $f_{i}:\mathbb{R}^{n}\rightarrow \mathbb{R}$ is a convex function which is only available to node $i$. To solve problem (<ref>), in [24] it is proposed a continuous time dynamic model based on controlling the sum of subgradients to force the state to the optimal solution set.

Let $y_{i}\in \mathbb{R}^{n}$ and $z_{i}\in \mathbb{R}^{n}$ be the state and auxiliary state of node $i$. It is assumed that the dynamics of each node are governed by the following ordinary differential equations:

\begin{eqnarray}\label{dyndist}
\dot{y}_{i}(t)=&&\sum _{j\in N_{i}}a_{ij}(y_{j}(t)-y_{i}(t))\\
\nonumber
&+&\sum _{j\in N_{i}}a_{ij}(z_{j}(t)-z_{i}(t))-g_{i}(y_{i}(t))\\
\nonumber
\dot{z}_{i}(t)=&&\sum _{j\in N_{i}}a_{ij}(y_{j}(t)-y_{i}(t)) ,
\end{eqnarray}

where $a_{ij}=a_{ji}>0$ are scalar weights chosen by node $i$, $g_{i}(t)$ is the subgradient of the function $f_{i}$ at the point $y_{i}(t)$.

The dynamical system <ref> can be written compactly as

\begin{eqnarray}\label{ccm001}
\dot{y}(t)&=&-\mathbf{L}y(t)-\mathbf{L}z(t)-G(y(t))\\
\nonumber
\dot{z}(t)&=&\mathbf{L}y(t),
\end{eqnarray}

where $\mathbf{L} = L \otimes I_{n}$ being $L$ the Laplacian matrix, and $G(t)$ is a concatenation of functions $g_{i}(t)$.

In the next section, we explore the diffusion dynamics in multiplex networks, in order to obtain the object which drives these dynamics on these physical systems, such that, we can use the control approach to generate the distributed-multilayer algorithms of interest.

§.§ Multiplex Networks

A multiplex network structure is a simplified multilayer network where nodes in different layers are mapped one-to-one, and interlinks exclusively connect to corresponding replica nodes [5]. Typically, a multiplex network is used to model interactions between the same set of nodes, with each layer representing a distinct type of interaction. However, multiplex networks can also represent interactions between different sets of nodes, as long as the nodes are mapped one-to-one.

A multiplex network can be represented in two different ways. In the first representation, there is no distinction between the identity of corresponding replica nodes, and interlinks are not explicitly used. In the second representation, corresponding replica nodes are treated as distinguishable entities, and interlinks are explicitly described. Although nodes in different layers of a multiplex network are typically mapped one-to-one, indicating the same identity, in some cases, it may still be useful to distinguish the identity of nodes in each layer.

In this paper, we focus on the second representation, which includes interlinks and treats corresponding replica nodes as distinguishable entities. This approach allows for a more explicit and detailed description of the multiplex network structure and facilitates the analysis of complex interactions between nodes in different layers. In this representation, a multiplex network is formed by $N$ nodes $i=1, 2, ..., N$ and $M$ layers, each node $i$ admits $M$ replica nodes $(i, \alpha)$ with $\alpha = 1, 2, ..., M$ indicating the identity of node $i$ in layer $\alpha$. Starting from the set of nodes $V$ of the multiplex network

\begin{equation*}
V = \{ i | i \in \{ 1, 2, ..., N \} \}.
\end{equation*}

We construct $M$ sets of nodes $V_{\alpha}$, each one representing the replicas of the nodes in $V$ in layer $\alpha$: $V_{\alpha}= \{ (i,\alpha)| i \in \{ 1, 2, ..., N \} \}$. Distinguishing between different replica nodes, the multiplex network can be represented as a multilayer network $\mathcal{M}$ comprising the triplet [5]

\begin{equation*}
\mathcal{M}=(Y, \mathbf{G}, \mathcal{G})
\end{equation*}

being $\mathbf{G}$ given by $\mathbf{G} = (G_{1}, G_{2}, ..., G_{\alpha}, ..., G_{M})$, and $G_{\alpha}=(V_{\alpha}, E_{\alpha})$ is the network in layer $\alpha$ and the $M\times M$ list $\mathcal{G}$ of bipartite networks $\mathcal{G}$ characterizes the interactions across pairs of different layers and has elements $\mathcal{G}_{\alpha,\beta}$ given by $\mathcal{G}_{\alpha ,\beta}=(V_{\alpha},V_{\beta},E_{\alpha ,\beta})$, where for each $\alpha < \beta$ and $\alpha,\, \beta \in \{ 1, 2,..., M \}$. $\mathcal{G}_{\alpha ,\beta}$ indicates the bipartite network with node sets $V_{\alpha}$ and $V_{\beta}$ and link set $E_{\alpha ,\beta}$. The links of the networks $\mathcal{G}_{\alpha,\beta}$ are called interlinks and connect the nodes of layer $\alpha$ to the nodes of a different layer $\beta$. In this particular representation, the interlinks in the set $E_{\alpha ,\beta}$ connecting the replica nodes in layer $\alpha$ with the replica nodes in layer $\beta$ are exclusively of the type $[(i,\alpha),(i,\beta)]$, i.e. they connect each replica
node $(i, \alpha)$ on layer $\alpha$ with its replica node $(i, \beta)$ on layer $\beta$,

\begin{equation*}
E_{\alpha ,\beta} = \{ [(i,\alpha),(i,\beta)]| \; i \in \{ 1, 2, ..., N \}\} .
\end{equation*}

In this multiplex network representation, both the intralinks and the interlinks can be indicated by the $N\cdot M \times N\cdot M$ supra-adjacency matrix $\mathcal{A}$ given by [5]

\begin{equation}
\mathcal{A}_{i\alpha,j\beta} =
\begin{cases}
a_{ij}^{[\alpha]}\quad \text{if} \quad \alpha = \beta, \\
\delta (i,j)\quad \text{if} \quad \alpha \neq \beta,
\end{cases}
\end{equation}

where $a_{ij}^{[\alpha]}$ are the matrix elements of the $M$ distinct adjacency matrices $A^{[\alpha]}$ of the network in layer $\alpha$. For example, for weighted, undirected multiplex networks elements $a_{ij}^{[\alpha]}$ are given by

\begin{equation}
a_{ij}^{[\alpha]} = \begin{cases}
w _{ij}^{[\alpha]}\quad \text{if node $i$ is linked to node $j$ in layer $\alpha$}\\
\; \qquad\text{with weight }\; \; w _{ij}^{[\alpha]}\\
0 \qquad \text{otherwise.}
\end{cases}
\end{equation}

In Fig. <ref>, it is shown a multiplex network with $M=2$ layers and $N=5$ nodes in each layer.

Multiplex network with $M=2$ layers and $N=5$ nodes in each layer.

§.§.§ Diffusion Processes on Multiplex Networks

Diffusion processes within each layer and across a different layer can be described by a dynamical diffusion state $y_{i}^{[\alpha]}(t)$ associated with each replica node $(i,\alpha)$ of the multiplex network with $i=1,2,...,N$ and $\alpha = 1,2,...,M$. We assume that both intralinks and interlinks can sustain the diffusion dynamics.

The interlayer diffusion is determined by the diffusion constant $D^{[\alpha]}$. Similarly, the interlayer diffusion is dictated by the diffusion constants $D^{[\alpha,\beta]}=D^{[\beta,\alpha]}$. The diffusion across different links of the same layer can be modulated by a weight $w _{ij}^{[\alpha]}$ associated with each undirected link from the replica node $(i,\alpha)$ to the replica node $(j,\alpha)$. With this notation, the general diffusion equation in a multiplex network is given by [38]

\begin{eqnarray}
\dot{y}_{i}^{[\alpha]} =  D^{[\alpha]}\sum _{j=1}^{N}w _{ij}^{[\alpha]} \left( y_{j}^{[\alpha]} - y_{i}^{[\alpha]} \right)  \nonumber\\
\qquad \qquad +\sum _{\beta =1}^{M} D^{[\alpha,\beta]} \left( y_{i}^{[\beta]} - y_{i}^{[\alpha]} \right),
\label{diffusion1}
\end{eqnarray}

where the first term and the second term on the right-hand side account respectively for the intralayer and the interlayer diffusion.

Equation (<ref>) can be written as a general diffusion equation in an $(N, M)$ dimensional space

\begin{equation}\label{diffMP}
\dot{y}=-\mathcal{L}y,
\end{equation}

where $\mathcal{L}$ is an $(N\cdot M)\times (N\cdot M)$ matrix called the supra-Laplacian matrix and $y$ is an $N\cdot M$ column vector encoding the dynamical state of each replica node of the multiplex network. In a multiplex network of $M$ layers the supra-Laplacian matrix $\mathcal{L}$ is given by

\begin{equation}\label{supraL}
\mathcal{L} = \Lambda + \Delta
\end{equation}

where matrices $\Lambda$ and $\Delta$ in (<ref>) are given by

\begin{equation}
\Lambda = \left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c,}
D^{[1]}L^{[1]} & 0 & \cdots & 0\\
0 & D^{[2]}L^{[2]} & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & D^{[M]}L^{[M]}%
\end{IEEEeqnarraybox*}\right)
\label{lambda1}
\end{equation}

\begin{equation}\label{delta1}
\begin{tiny}
\Delta = \left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c,}
\sum _{\beta \neq 1}D^{[1,\beta]}I_{N\times N} & -D^{[1,2]}I_{N\times N} & \cdots & -D^{[1,M]}I_{N\times N}\\
-D^{[2,1]}I_{N\times N} & \sum _{\beta \neq 2}D^{[2,\beta]}I_{N\times N} & \cdots & -D^{[1,M]}I_{N\times N}\\
\vdots & \vdots & \ddots & \vdots \\
-D^{[M,1]}I_{N\times N} & -D^{[M,2]}I_{N\times N} & \cdots & \sum _{\beta \neq M}D^{[M,\beta]}I_{N\times N}%
\end{IEEEeqnarraybox*}\right) ,
\end{tiny}
\end{equation}

respectively, where $I_{N\times N}$ is the $N\times N$ identity matrix. In matrix (<ref>) $L^{[\alpha]}$ indicates the Laplacian matrix in each layer whose elements are $L _{ij}^{[\alpha]} = s _{i}^{[\alpha]}\delta _{ij}-w _{ij}^{[\alpha]}$,
being $s _{i}^{[\alpha]}$ the strength of the replica node $(i,\alpha)$ defined by $s _{i}^{[\alpha]} = \sum _{j}w _{ij}^{[\alpha]}$. Similarly, the dynamical vector $y$ can be written as $y = \left(y^{[1]} \; y^{[2]}\; \ldots  y^{[M]}\right)^{\top}$, where the $y^{[\alpha]}$ indicates the $N$ column vector of elements $y^{[\alpha]}_{i}$ with $i=1,2,...,N$. In an undirected multiplex network, where the diffusion constant $D^{[\alpha ,\beta ]}=D^{[\beta ,\alpha ]}$, the supra-Laplacian matrix $\mathcal{L}$ is symmetric and semi-positive definite. Therefore, its eigenvalues $\Lambda _{n}$ are real and non-negative [38]. Here we shall assume that the multiplex network, including its interlinks, is connected. We assume that from a replica node, it is possible to reach any other node by following a combination of intralinks and interlinks.

§ PROBLEM STATEMENT

We consider the following distributed optimization problem

\begin{equation}
\min _{x \in \mathbb{R}}  f(x)= \sum _{\alpha = 1}^{M} \sum _{i=1}^{N}f_{i}^{[\alpha]}(x),
\label{opproblem1}
\end{equation}

where $f_{i}^{[\alpha]}:\mathbb{R} \rightarrow \mathbb{R}$ is a private cost function known only to agent $(i,\alpha)$ in the layer $\alpha$. We assume that each $f_{i}^{[\alpha]}$ is strictly convex and twice continuously differentiable. Note that each agent's cost is a function of the global variable $x$ across all layers, and it is the reason that agents need to coordinate to solve problem (<ref>).

To find a distributed approach for solving the optimization problem (<ref>), we associate a local estimate $y_{i}^{[\alpha]} \in \mathbb{R}$ of the global variable $x \in \mathbb{R}$ to every agent and solve the equivalent problem

\begin{eqnarray}
&&\min _{y\in \mathbb{R}^{N\cdot M}}\tilde{f}(y)= \sum _{\alpha = 1}^{M} \sum _{i=1}^{N}f_{i}^{[\alpha]}(y_{i}^{[\alpha]}) \nonumber \\
&&\text{subject to}\quad \mathcal{C}(y) ,
\label{opproblem2}
\end{eqnarray}

where $\mathcal{C}(y)$ is a consensus constraint that assures that the local estimates of all processors coincide. There are different choices for the consensus constraint. We shall use the supra-Laplacian constraint $0 _{N\cdot M}=\mathcal{L}y$.

We begin our analysis of the distributed optimization problem (<ref>) by considering the supra-Laplacian-based consensus constrained problem

\begin{eqnarray}\label{opproblem3}
&&\min _{y\in \mathbb{R}^{N\cdot M}}\tilde{f}(y)= \sum _{\alpha = 1}^{M} \sum _{i=1}^{N}f_{i}^{[\alpha]}(y_{i}^{[\alpha]}) \nonumber \\
&&\text{subject to}\quad \mathcal{L}y =  0 _{N\cdot M}.
\end{eqnarray}

In the next sections, we propose two different saddle-point algorithms to solve the distributed optimization problem with supra-Laplacian-based constraints.

§ DISTRIBUTED PRIMAL-DUAL SADDLE-POINT ALGORITHM

For problem (<ref>), we propose the following augmented Lagrangian function $\mathcal{L}_{a}: \mathbb{R}^{N\cdot M}\times \mathbb{R}^{N\cdot M}\rightarrow \mathbb{R}$

\begin{equation}
\mathcal{L}_{a}(y,\lambda) = \tilde{f}(y)+\lambda ^{T}\mathcal{L}y +\frac{1}{2}y^{T}\mathcal{L}y ,
\label{augmented1}
\end{equation}

and the associated saddle-point flow yields

\begin{equation}
\dot{y}=-\frac{\partial \mathcal{L}_{a}(y,\lambda)}{\partial y}=-\frac{\partial \tilde{f}(y)}{\partial y}-\mathcal{L}y-\mathcal{L}\lambda
\label{saddleflow1}
\end{equation}

\begin{equation}
\dot{\lambda}=\frac{\partial \mathcal{L}_{a}(y,\lambda)}{\partial \lambda}= \mathcal{L}y .
\label{saddleflow2}
\end{equation}

In Fig. <ref> we depict the block-diagram representation of the saddle-point flow (<ref>) - (<ref>) including further gains for the augmentation, primal, and dual dynamics. The first three upper blocks are referring to the primal gradient dynamics. The following two, are taking the multiplex damping via augmentation $\rho$. The final two bottom blocks are the dual integrator dynamics. The supra-Laplacian blocks, in the input and output of the dual integrator dynamics, are the skew-symmetric multiplex coupling.

The saddle points of the augmented Lagrangian (<ref>) and the equilibrium of the saddle-point flow (<ref>)-(<ref>) are characterized by the following lemma.

Let $\mathcal{L}=\mathcal{L}^{T}\in \mathbb{R}^{N\cdot M \times N\cdot M}$ be a symmetric supra-Laplacian associated with an undirected, connected, and weighted multiplex network. Consider the augmented Lagrangian function (<ref>), where each $f_{i}^{[\alpha]}$ is strictly convex and twice continuously differentiable for all $(i,\alpha)$ where $i\in \{ 1,...,N \}$ and $\alpha\in \{ 1,...,M \}$. Then

1. If $(y^{\star},\lambda ^{\star})\in \mathbb{R}^{N\cdot M}\times \mathbb{R}^{N\cdot M}$ is a saddle point of (<ref>), then so is $(y^{\star},\lambda ^{\star}+\gamma 1_{N\cdot M})$ for any $\gamma \in \mathbb{R}$.

2. If $(y^{\star},\lambda ^{\star})\in \mathbb{R}^{N\cdot M}\times \mathbb{R}^{N\cdot M}$ is a saddle point of (<ref>), then $y^{\star} = x^{\star}1_{N\cdot M}$ where $x^{\star} \in \mathbb{R}$ is a solution of the original optimization problem (<ref>).

3. If $x^{\star} \in \mathbb{R}$ is a solution of the original optimization problem (<ref>), then there are $\lambda ^{\star}\in \mathbb{R}^{N\cdot M}$ and $y^{\star} = x^{\star}1_{N\cdot M}$ satisfying

\begin{equation*}
\mathcal{L}\lambda ^{\star}+\frac{\partial \tilde{f}(y^{\star})}{\partial y}=0_{N\cdot M}
\end{equation*}

so that $(y^{\star},\lambda ^{\star})$  is a saddle point of (<ref>).

It follows from the Karush-Kuhn Tucker conditions applied to the Lagrangian function (<ref>) and the definition of a saddle-point in convex optimization [37].

Now, it is presented the convergence analysis of the proposed saddle-point dynamics for multiplex networks (<ref>)-(<ref>). The main result for convergence analysis of the proposed primal-dual saddle-point dynamics is presented in the following theorem.

[auto, node distance=2cm,>=latex']

[sum, right of=no, node distance=0.7cm](sumo);
[block, right of=sumo, node distance = 1.5cm] (ko) $K_{1}$;
[coordinate, right of=ko, node distance=1.5cm](nt);
[block, right of=nt, node distance=1.5cm](ss)$\frac{1}{s}$;
[coordinate, right of=ss, node distance=1.5cm](ntr);
[coordinate, right of=ntr, node distance=0.7cm](nf);

[coordinate, below of=no, node distance=2cm](nfi);
[sum, below of=sumo, node distance=2cm](sumt);
[coordinate, below of=ko, node distance = 2cm] (ns) ;
[block, below of=nt, node distance=2cm](grad)$\frac{\partial \tilde{f}}{\partial y}$;
[coordinate, below of=ss, node distance=2cm](ne);
[coordinate, below of=ntr, node distance=2cm](nn);
[coordinate, below of=nf, node distance=2cm](nten);

[coordinate, below of=nfi, node distance=2cm](nel);
[coordinate, below of=sumt, node distance=2cm](ntw);
[block, below of=ns, node distance = 2cm] (rho) $\rho$;
[coordinate, below of=grad, node distance=2cm](ntrh);
[block, below of=ne, node distance=2cm](suplapo)$\mathcal{L}$;
[coordinate, below of=nn, node distance=2cm](nfrth);
[coordinate, below of=nten, node distance=2cm](nfith);

[block, below of=nel, node distance=2cm](suplapt)$-\mathcal{L}^{T}$;
[coordinate, below of=ntw, node distance=2cm](nsith);
[coordinate, below of=rho, node distance = 2cm] (nseth) ;
[coordinate, below of=ntrh, node distance=2cm](neith);
[coordinate, below of=suplapo, node distance=2cm](ninth);
[coordinate, below of=nfrth, node distance=2cm](ntwww);
[block, below of=nfith, node distance=2cm](suplaptr)$\mathcal{L}$;

[coordinate, below of=suplapt, node distance=2cm](ntwo);
[coordinate, below of=nsith, node distance=2cm](nttr);
[block, below of=nseth, node distance = 2cm] (sso) $\frac{1}{s}$;
[coordinate, below of=neith, node distance=2cm](ntfr);
[block, below of=ninth, node distance=2cm](kt)$K_{2}$;
[coordinate, below of=ntwww, node distance=2cm](ntfi);
[coordinate, below of=suplaptr, node distance=2cm](ntsix);

[->] (no) – node  (sumo);
[->] (sumo) – node  (ko);
[-] (ko) – node  (nt);
[->] (nt) – node $\dot{y}$ (ss);
[-] (ss) – node $y$ (ntr);
[-] (ntr) – node  (nf);
[->] (ns) – node[pos=0.99] $+$ node [near end]  (sumt);
[-] (grad) – node  (ns);
[->] (ne) – node  (grad);
[-] (nn) – node  (ne);
[->] (sumt) – node[pos=0.99] $-$ node [near end]  (sumo);
[-] (ntr) – node  (nn);
[-] (rho) – node  (ntw);
[->] (ntw) – node  (sumt);
[->] (suplapo) – node  (rho);
[->] (nfrth) – node  (suplapo);
[-] (nfrth) – node  (nn);
[-] (suplapt) – node  (no);
[->] (nf) – node  (suplaptr);
[-] (suplaptr) – node  (ntsix);
[->] (ntsix) – node  (kt);
[->] (kt) – node $\dot{\lambda}$ (sso);
[-] (sso) – node $\lambda$ (ntwo);
[->] (ntwo) – node  (suplapt);

Block-diagram representation of the saddle-point flow (<ref>) - (<ref>) including gains for the augmentation, primal, and dual dynamics.

Let $\tilde{f}: \mathbb{R}^{N\cdot M} \rightarrow \mathbb{R}$ be twice continuously differentiable and strictly convex. Then each solution $t\mapsto (y(t), \lambda (t))$ of the distributed saddle-point flow (<ref>)-(<ref>) converges to a saddle point $(y^{\star}, \lambda^{\star})$ uniquely characterized by

1. $y^{\star} = x^{\star}1_{N\cdot M}$ where $x^{\star} \in \mathbb{R}$ is a solution of the original optimization problem (<ref>).

2. $\lambda^{\star} = \bar{\lambda}+\text{average}(\lambda _{0})1_{N\cdot M}$ where $\bar{\lambda}\perp 1_{N\cdot M}$ satisfies

\begin{equation*}
\mathcal{L}\bar{\lambda}+\frac{\partial \tilde{f}(y^{\star})}{\partial y}=0_{N\cdot M}
\end{equation*}

For the general convergence to the set of saddle points $(y^{\star}, \lambda^{\star})$, let $\tilde{y}=y-y^{\star}$ and $\tilde{\lambda}=\lambda -\lambda$.  In the new variables, the saddle-point flow (<ref>)-(<ref>) yields

\begin{equation}
\dot{\tilde{y}}=-\frac{\partial \tilde{f}(y)}{\partial y}+ \frac{\partial \tilde{f}(y^{\star})}{\partial y}-\mathcal{L}\tilde{y}-\mathcal{L}\tilde{\lambda},
\label{saddleflow3}
\end{equation}

\begin{equation}
\dot{\tilde{\lambda}}=\mathcal{L}\tilde{y} .
\label{saddleflow4}
\end{equation}

Consider the following quadratic candidate Lyapunov function

\begin{equation}
\label{lyapunov1}
\end{equation}

Then, the temporal derivative of (<ref>) along the trajectories of (<ref>)-(<ref>) is

\begin{equation}
\dot{V}(\tilde{y},\tilde{\lambda})=-\tilde{y}^{T}\frac{\partial \tilde{f}(y)}{\partial y}+\tilde{y}^{T}\frac{\partial \tilde{f}(y^{\star})}{\partial y}-\tilde{y}^{T}\mathcal{L}\tilde{y}.
\end{equation}

From the global under-estimator property of the gradient, we have that

\begin{equation}
\tilde{f}(y^{\star})\geq \tilde{f}(y)+(y^{\star}-y)^{T}\frac{\partial \tilde{f}(y)}{\partial y} .
\label{ineqproof}
\end{equation}

From inequality (<ref>), and if $\tilde{f}$ is strictly convex, we have

\begin{equation*}
-\tilde{y}^{T}\frac{\partial \tilde{f}(y)}{\partial y}+\tilde{y}^{T}\frac{\partial \tilde{f}(y^{\star})}{\partial y} < 0.
\end{equation*}

Now, since the supra-Laplacian matrix is semi-positive definite, we obtain that

\begin{equation*}
\dot{V}(\tilde{y},\tilde{\lambda})=-\tilde{y}^{T}\frac{\partial \tilde{f}(y)}{\partial y}+\tilde{y}^{T}\frac{\partial \tilde{f}(y^{\star})}{\partial y}-\tilde{y}^{T}\mathcal{L}\tilde{y} <0.
\end{equation*}

However, $\dot{V}(\tilde{y},\tilde{\lambda})=0$ on the set $E=\{ (\tilde{y},\tilde{\lambda})|\tilde{y}=0 \}$. To show that the equilibrium point $(y^{\star},\lambda ^{\star})$  is globally asymptotically stable, we invoke LaSalle invariant principle and show that there are no trajectories in $E$ besides the equilibrium point $\tilde{y}=0, \, \tilde{\lambda}=0$.

When $\tilde{y}=0$ the system (<ref>)-(<ref>) reduces to

\begin{equation*}
\dot{\tilde{y}}=-\mathcal{L}\tilde{\lambda}
\end{equation*}

\begin{equation*}
\dot{\tilde{\lambda}}=0.
\end{equation*}

It follows that since the supra-Laplacian matrix $\mathcal{L}$ is semi-positive definite and by Lemma<ref> the solution is $(0, \gamma 1_{N\cdot M})$ or $(0,0)$. Furthermore, from Lemma <ref>, the saddle points $(y^{\star}, \lambda^{\star})$ are identified as $y^{\star}=x^{\star}1_{N\cdot M}$ where $x^{\star}\in \mathbb{R}$ is the solution of the original optimization problem (<ref>) and $\lambda^{\star} = \bar{\lambda}+\gamma 1_{N\cdot M}$ with $\gamma \in \mathbb{R}$ and $\bar{\lambda}\perp 1_{N\cdot M}$ satisfying

\begin{equation*}
\mathcal{L}\bar{\lambda}+\frac{\partial \tilde{f}(y^{\star})}{\partial y}=0_{N\cdot M}.
\end{equation*}

From the dual dynamics (<ref>), we have that $1_{N\cdot M}\dot{\lambda}=0$ and therefore $\text{average}(\lambda (t))=\text{average}(\lambda _{0})$ $\forall t \geq 0$. It follows that the saddle-point dynamics (<ref>)-(<ref>) converge to a unique saddle-point $(y^{\star}, \lambda ^{\star})$ satisfying $y^{\star}=x^{\star}1_{N\cdot M}$ and $\lambda^{\star} = \bar{\lambda}+ \text{average}(\lambda _{0}) 1_{N\cdot M}$

§ DISTRIBUTED GRADIENT DESCENT FOR MULTIPLEX NETWORKS

The distributed gradient flow with supra-Laplacian coupling reads

\begin{equation}
K_{1}\dot{y}=-\frac{\partial \tilde{f}(y)}{\partial y}-\rho \mathcal{L}y .
\label{distributedgflow1}
\end{equation}

From the soft-penalty method, the distributed gradient flow (<ref>) can be viewed as the regular gradient flow of the following optimization problem

\begin{equation}
\min _{y\in \mathbb{R}^{N\cdot M}} \sum _{\alpha =1}^{M} \sum _{i=1}^{N} f_{i}^{[\alpha]}(y_{i}^{[\alpha]}) + \frac{1}{2}\rho y^{T}\mathcal{L}y ,
\label{opproblem5}
\end{equation}

where the soft-penalty $\frac{1}{2}\rho y^{T}\mathcal{L}y$ is tasked with enforcing a supra-Laplacian consensus-constraint. Despite the minimizer of (<ref>) being different from that of the distributed optimization problem (<ref>), the following result shows that the flow (<ref>) converges to the correct solution of the distributed optimization problem (<ref>) provided that the gains are time-varying. In Fig. <ref> we show the block-diagram representation of the distributed gradient flow (<ref>). Notice that it is a constitutive part of the control scheme presented in Fig. <ref>.

Block-diagram representation of the gradient flow  (<ref>).

Consider the time-varying flow

\begin{equation}
\dot{y}=-\varsigma (t)\frac{\partial \tilde{f}(y)}{\partial y}-\mathcal{L}y ,
\label{distgflow2}
\end{equation}

where the time-varying positive gain $\varsigma:\mathbb{R}_{\geq 0}\rightarrow \mathbb{R}_{>0}$ satisfies the persistence condition

\begin{equation}
\varsigma (t)>0, \quad \int _{0}^{\infty}\varsigma (\tau)d\tau = \infty ,\quad \text{and}\, \lim _{t\rightarrow \infty}\varsigma (t) = 0,
\label{persistencec1}
\end{equation}

that is, over time the supra-Laplacian coupling in the distributed gradient flow (<ref>) dominates the gradient descent with a decaying yet non-integrable gain $\varsigma (t)$.

The convergence of the distributed gradient flow for multiplex networks is given by the following theorem.

Consider the gradient flow (<ref>) where $\varsigma: \mathbb{R}_{\geq 0}\rightarrow \mathbb{R}_{>0}$ satisfies the persistence condition (<ref>), and $\tilde{f}:\mathbb{R}^{N\cdot M}\rightarrow \mathbb{R}$ is assumed to be twice continuously differentiable, strictly convex, and radially unbounded. Then each solution $t\mapsto y(t)$ of the distributed gradient flow (<ref>) converges to the unique optimizer $y^{\star}=x^{\star}1_{N\cdot M}$.

Let us first study the boundedness of the state $y$ by considering the following Lyapunov function

\begin{equation}
V(y)=\frac{1}{2}\| y - y^{*}\| _{2}^{2} .
\label{lyap51}
\end{equation}

The time derivative of (<ref>) along the trajectories of (<ref>) is given by

\begin{equation*}
\dot{V}(y)=-\varsigma(t)(y-y^{*})^{T}\frac{\partial \tilde{f}(y)}{\partial y}-(y-y^{*})^{T}\mathcal{L}y.
\end{equation*}

Taking advantage of the of under estimator property (<ref>), and that $\mathcal{L}y = \mathcal{L}x^{*}1_{N\cdot M}=0$, the last expression yields

\begin{equation*}
\dot{V}(y)=-\varsigma(t)(\tilde{f}(y)-\tilde{f}(y^{*}))-(y-y^{*})^{T}\mathcal{L}(y-y^{*}).
\end{equation*}

Thus, $\dot{V}(y)\leq 0$, and then, the state $y$ is bounded. Therefore, the distance between the state $y$ and the optimizer is given by

\begin{equation}
\| y(t) - y^{*}\| \leq \| y _{0} - y^{*}\| .
\label{bound51}
\end{equation}

Within this result, we shall evaluate the asymptotic consensus and optimality via the following coordinate transformation

\begin{equation}
\begin{matrix}
\left[\begin{IEEEeqnarraybox*}[][c]{,c,}
\psi \\
\eta %
\end{IEEEeqnarraybox*}\right]
\end{matrix} =
\begin{matrix}
\left[\begin{IEEEeqnarraybox*}[][c]{,c,}
U ^{T} \\
1_{N\cdot M}^{T}/(N\cdot M) %
\end{IEEEeqnarraybox*}\right] y
\end{matrix}
\end{equation}

Being $U$ an orthogonal matrix such that

\begin{equation}
U^{T}\mathcal{L}U = -\text{diag}(\lambda _{2},..., \lambda _{N\cdot M}) .
\end{equation}

The dynamics of the $\psi (t)$ function reads

\begin{equation}
\dot{\psi} = - \text{diag}(\lambda _{2},..., \lambda _{N\cdot M}) \psi (t) + \xi (t)
\label{psi111}
\end{equation}

being $\xi (t)=-\varsigma (t) U^{T}\frac{\partial \tilde{f}(y)}{\partial y}$.

The input signal (or source term) in (<ref>) is bounded due to (<ref>) and it is decaying to zero by the properties established for $\varsigma (t)$. Then, the dynamics driven by (<ref>) are an exponentially stable linear system subject to the input $\xi(t)$. Therefore, 	$\psi (t) \rightarrow 0$ as $t \rightarrow \infty$, and $y(t)\rightarrow \eta $, that is, $\psi (t\rightarrow \infty)\in \text{span}(1_{N\cdot M})$, achieving asymptotically the consensus.

With the remaining coordinate $\eta = \text{average}(y)$, we shall study the asymptotic optimality. The dynamics driving this coordinate are given by

\begin{eqnarray}
\dot{\eta} &=& -\frac{\varsigma (t)}{N\cdot M}1_{N\cdot M}^{T}\frac{\partial \tilde{f}(y)}{\partial y} \nonumber
\\
&=&  -\frac{\varsigma (t)}{N\cdot M} \sum _{i = 1}^{N\cdot M}\frac{\partial f _{i}(y_{i})}{\partial y_{i}} \nonumber
\\
&=&  -\frac{\varsigma (t)}{N\cdot M} \sum _{i = 1}^{N\cdot M}\frac{\partial f _{i}}{\partial y_{i}}((U\psi)_{i}+\eta) .
\end{eqnarray}

Considering the following time scaling

\begin{equation*}
\tau = \int _{0}^{t}\varsigma (t')dt'\, \Rightarrow \, \frac{d\tau}{dt}= \varsigma (t),
\end{equation*}

we have that $\tau$ is strictly monotonically increasing in $t$ due to $\varsigma (t) > 0$. On the other hand, the integral persistence condition makes that $t \in [0, \infty]$ is mapped to $\tau = [0, \infty]$. Thus, the time-coordinate transformation is invertible. In the $\tau$-time-scale, the dynamics which are driving the coordinate $\eta$ are given by

\begin{equation}
\frac{d}{d\tau}\eta = -\frac{1}{n} \sum _{\alpha=1}^{M}\sum _{i=1}^{N}\frac{\partial f _{i}^{[\alpha]}}{\partial y_{i}^{[\alpha ]}}((U\psi)_{i}+\eta) .
\end{equation}

With this in mind, let us consider the following Lyapunov function

\begin{equation}
W(\eta) = \frac{n}{2}(\eta -\eta ^{*})^{2} ,
\label{lyapeta51}
\end{equation}

being $\eta ^{*}=x^{*}$ is the global minimizer of $f(x)=\sum _{\alpha=1}^{M}\sum _{i=1}^{N}f_{i}^{[\alpha]}(x)$.

The temporal-scaled derivative of (<ref>) is

\begin{eqnarray*}
\frac{d}{d\tau}W(\eta) = &-&(\eta - \eta ^{*})\sum _{\alpha=1}^{M}\sum _{i=1}^{N}\frac{\partial f_{i}^{[\alpha]}}{\partial y_{i}^{[\alpha]}}(\eta) \nonumber
\\
&+& (\eta ^{*}-\eta)\sum _{\alpha=1}^{M}\sum _{i=1}^{N}\left[ \frac{\partial f_{i}^{[\alpha]}}{\partial y_{i}^{[\alpha]}}((U\psi)_{i}^{[\alpha]}+\eta) \right.
\\
&-&\left. \frac{\partial f_{i}^{[\alpha]}}{\partial y_{i}^{[\alpha]}}(\eta) \right] .
\end{eqnarray*}

By the under estimator property, the first term of $\frac{d}{d\tau}W(\eta)$ is upper-bounded by $f(\eta ^{*})-f(\eta)$. The second term of $\frac{d}{d\tau}W(\eta)$ can be upper-bounded as follows,

\begin{eqnarray*}
(\eta ^{*}-\eta)\sum _{\alpha=1}^{M}\sum _{i=1}^{N}\left[ \frac{\partial f_{i}^{[\alpha]}}{\partial y_{i}^{[\alpha]}}((U\psi)_{i}^{[\alpha]}+\eta) \right. \nonumber
\\
-\left. \frac{\partial f_{i}^{[\alpha]}}{\partial y_{i}^{[\alpha]}}(\eta) \right] \leq | \eta ^{*}-\eta | \sum _{\alpha=1}^{M}\sum _{i=1}^{N} L_{i}^{[\alpha]}|(U\psi)_{i}^{[\alpha]} |.
\label{ineq511}
\end{eqnarray*}

It can be done due to the boundedness of the state $y$. In the inequality (<ref>) the constant $L_{i}^{[\alpha]}$ is the maximal Lipschitz constant of $\partial f _{i}^{[\alpha]}\partial y$ evaluated in the compact invariant set $\{ y \in \mathbb{R}^{N\cdot M} | \| y - y^{*}\| \leq \| y_{0} - y^{*}\| \}$, defined by

\begin{equation}
L_{i}^{[\alpha]}=\max _{\| y - y^{*}\| \leq \| y_{0} - y^{*}\|} \text{Hess}f_{i}^{[\alpha]}(y) .
\end{equation}

Due to the asymptotic consensus, for each $\epsilon \geq 0$, there is a $\tau _{\epsilon}$ so that for all $\tau \geq \tau _{\epsilon}$ we have

\begin{equation*}
\sum _{\alpha=1}^{M}\sum _{i=1}^{N} L_{i}^{[\alpha]}|(U\psi)_{i}^{[\alpha]}|\leq \epsilon.
\end{equation*}

That is, for each $\tau \geq \tau _{\epsilon}$

\begin{equation}
\frac{d}{d\tau}W(\eta) \leq f(\eta ^{*})-f(\eta)+\epsilon |\eta ^{*}-\eta | \leq 0.
\end{equation}

Since $\tilde{f}$ (and also $f$) is radially unbound and strictly convex, we have that there is a $|\eta _{\epsilon}|$ such that the temporal-scaled derivative of (<ref>) is strictly negative for sufficiently large $|\eta ^{*}-\eta | \geq |\eta _{\epsilon}|$ and therefore $W(\eta (\tau \rightarrow \infty ))\leq W(\eta _{\epsilon})$. Additionally, $\eta _{\epsilon} \rightarrow 0$ as $\epsilon \rightarrow 0$. Thus, for each $\tilde{\epsilon}\geq 0$, there exists an $\epsilon$ sufficiently small so that for $\tau \geq \tau _{\epsilon}$, $W(\eta (\tau))\leq \tilde{\epsilon}$. Therefore, $W(\eta (\tau))$ converges to $0$, or, $\eta (\tau)$ converges to the optimal value $\eta ^{*}$ as $\tau \rightarrow \infty$. Also occurs for $\eta (t)$ as $t\rightarrow \infty$.

In the next section, we illustrate the effectiveness of the proposed saddle-point algorithms in several numerical experiments.

§ NUMERICAL EXPERIMENTS

§.§ Multiplex Network with Two Layers and Three Nodes in Each Layer

We consider a multiplex network with $M=2$ layers and $N=3$ nodes in each layer. See Fig. <ref>.

Multiplex network topology with $M=2$ layers and $N=3$ nodes in each layer for the distributed optimization problem (<ref>).

We study the distributed optimization problem

\begin{equation}
\min _{x \in \mathbb{R}}  f(x)= \sum _{\alpha = 1}^{2} \sum _{i=1}^{3}f_{i}^{[\alpha]}(x),
\label{opproblemex3}
\end{equation}

via the primal-dual saddle-point algorithm to solve the following optimization problem

\begin{eqnarray}
&&\min _{y\in \mathbb{R}^{3\cdot 2}}\tilde{f}(y)=\sum _{i=1}^{3}f^{[1]}(y_{i}^{[1]})+\sum _{i=1}^{3}f^{[2]}(y_{i}^{[2]})\nonumber \\
&&\text{subject to}\quad \mathcal{L}y =  0 _{3\cdot 2},
\label{opproblem3b}
\end{eqnarray}

where $f^{[1]}(y_{i}^{[1]})=\frac{1}{2} (y_{i}^{[1]})^{2}+iy_{i}^{[1]}$, $f^{[2]}(y_{i}^{[2]})=\frac{1}{2} (y_{i}^{[2]})^{2}+(i+3)y_{i}^{[2]}$, and the supra-Laplacian matrix, with $D^{[1]}=D^{[2]}=D^{[1,2]}=1$ reads

\begin{equation}
\mathcal{L} = \left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c/c/c,}
3 & -1 & -1 & -1 & 0 & 0 \\
-1 & 3 & -1 & 0 & -1 & 0 \\
-1 & -1 & 3 & 0 & 0 & -1 \\
-1 & 0 & 0 & 3 & -1 & -1 \\
0 & -1 & 0 & -1 & 3 & -1 \\
0 & 0 & -1 & -1 & -1 & 3 %
\end{IEEEeqnarraybox*}\right) .
\label{lapex3}
\end{equation}

The results of the optimal consensus are shown in Fig. <ref>. We can observe a rapid convergence to the optimal $y^{\star}$ due to the value of the interdiffusion constant $D^{[1,2]}$ and the number of nodes in each layer.

Simulation results for the distributed optimization problem (<ref>) with supra-Laplacian matrix (<ref>).

Now, to the same problem, we apply the distributed gradient descent considering the time-varying flow (<ref>) and a time-varying positive gain $\varsigma (t,\theta)$ defined by

\begin{equation*}
\varsigma (t,\theta)=\frac{1}{\theta + t}.
\end{equation*}

The results of the time-varying flow (<ref>) associated with the optimization problem (<ref>) with supra-Laplacian (<ref>) are depicted in Fig. <ref>. Simulations are showing an interesting result, as $\theta \gg 1$ the optimal value is getting closer to the one obtained by the distributed primal-dual saddle-point algorithm; even, with a huge $\theta$ the consensus times are also getting similar.

Numerical results for the distributed optimization problem (<ref>) resolved by distributed gradient descent considering the time-varying flow (<ref>) and a time-varying positive gain: a) $\varsigma(t,0.1)$, b) $\varsigma(t,1)$, c) $\varsigma(t,10)$, d) $\varsigma(t,100)$, e) $\varsigma(t,1000)$, and f) $\varsigma(t,10000)$. The conventions are the same as presented in Fig. <ref>.

§.§ Multiplex Network with Two Layers and Five Nodes in Each Layer

We consider a multiplex network with $M=2$ layers and $N=5$ nodes in each layer as it is shown in Fig. <ref>. We study the distributed optimization problem

\begin{equation}
\min _{x \in \mathbb{R}}  f(x)= \sum _{\alpha = 1}^{2} \sum _{i=1}^{5}f_{i}^{[\alpha]}(x),
\label{opproblemex5}
\end{equation}

via the primal-dual saddle-point algorithm for the optimization problem

\begin{eqnarray}
&&\min _{y\in \mathbb{R}^{5\cdot 2}}\tilde{f}(y)=\sum _{i=1}^{5}f^{[1]}(y_{i}^{[1]})+\sum _{i=1}^{5}f^{[2]}(y_{i}^{[2]})\nonumber \\
&&\text{subject to}\quad \mathcal{L}y =  0 _{5\cdot 2},
\label{opproblem5b}
\end{eqnarray}

where $f^{[1]}(y_{i}^{[1]})=\frac{1}{2} (y_{i}^{[1]})^{2}$, $f^{[2]}(y_{i}^{[2]})=\frac{1}{2} (y_{i}^{[2]})^{2}$, and the supra-Laplacian matrix with $D^{[1]}=D^{[2]}=D^{[1,2]}=1$ can be obtained similar to (<ref>).

The results of the optimal consensus are shown in Fig. <ref>. We can observe a rapid convergence to the optimal $y^{\star}$ due to the value of the interdiffusion constant $D^{[1,2]}$ and the number of nodes in each layer.

Simulation results for the distributed optimization problem (<ref>).

Now, we shall study the behavior of the flow solutions when $D^{[1,2]}=0.1$.  In this case, the supra-Laplacian reads

\begin{equation}
\mathcal{L} = -\left(\begin{IEEEeqnarraybox*}[][c]{,c/c,}
\mathcal{L}_{5}^{[1]} & 0.1I_{5\times 5}\\
0.1I_{5\times 5} & \mathcal{L}_{5}^{[2]}
\end{IEEEeqnarraybox*}\right) .
\label{lapex501}
\end{equation}

In Fig. <ref> we depict the solutions of the distributed optimization problem (<ref>) applying the saddle-point equations given by (<ref>) and (<ref>). As is expected, the convergence takes more time than when we have $D^{[1,2]}=1$. However, we can observe a particular issue: As $D^{[1,2]}$ goes to zero, each layer is getting its own consensus, and after that, they search for the global optimal solution. From Fig. <ref>a, we can evidence that the more connected is the layer, the faster the intralayer consensus is reached. Comparing both figures and the result obtained with the topology of 3 nodes in each layer, the topology and number of nodes are showing a special relationship within the interlayer diffusion constant in the optimal consensus dynamics. Similar behavior can be obtained by manipulating the intralayer diffusion constants $D^{[1]}$ and $D^{[2]}$. Therefore, the parameters $D^{[1]}$, $D^{[2]}$, and $D^{[1,2]}$ serve as control parameter for the convergence.

Simulation results for the distributed optimization problem (<ref>).

§.§ Consensus Convergence to the Optimal $y^{\star}$ Time $t_{c}$

These results, together with the fact that the interlayer diffusion constant $D_{x}$ is driving a critical phenomenon referred to as the second eigenvalue of the supra-Laplacian matrix studied (connected with diffusion times), motivate us to study the behavior of the consensus convergence to the optimal $y^{\star}$ times $t_{c}$ as a function of the interlayer diffusion constant $D_{x}$. In Fig. <ref> we depict the time at which each $y_{i}^{[\alpha]}$ does reach the optimal for a two-layer multiplex network, wherein the first layer we have a ring topology, and in the second one we have a full-connected network. We can observe that, as it was expected: First, the nodes will reach the optimal faster as $D_{x}\rightarrow \infty$. Second, we have found values of $D_{x}$ where there is present a discontinuity. That is to say, we have encountered a situation that is validating our hypothesis about the possible presence of critical phenomena driven by $D_{x}$ and its relationship with the diffusion dynamics, which at the same time, is related to the consensus and optimization processes in the multiplexes. As we can evidence in Fig. <ref>, there are points where there are sudden changes that do not obey certain continuity but are presenting an optimization time process.

Simulation results for the distributed optimization problem (<ref>).

§.§ Coordinated Dispatch for Multienergy Network

As a final example, let us explore an applied situation where 2 + 1 restrictions are present. We study the case of three layers of micro-grids of seven power generators (each one) which are connected between them. Such is the case of coordinate dispatch for interdependent infrastructure networks. The multilayer dispatch problem seeks to minimize the total cost of operating power and natural gas networks. The global objective function and power and gas constraints are defined as follows

\begin{equation}
\min _{x \in \mathbb{R}}\sum _{\alpha = 1}^{3}\sum _{i = 1}^{7}f_{i}^{[\alpha]}(x) ,
\end{equation}

which translated to primal-dual saddle-point language, gives us the problem of distributed power optimization

\begin{equation*}
\min_{q} \left( \sum_{\tau \in \mathcal{T}}C_{\tau}^{[\mathcal{T}]}(p_{\tau}) + \sum_{k \in \mathcal{K}}C_{k}^{[\mathcal{K}]}(p_k) + \sum_{j \in \mathcal{G}}C_{j}^{[\mathcal{G}]}(g_j)\right) ,
\end{equation*}

subject to

\begin{equation*}
\mathcal{L}q = 0,
\end{equation*}

\begin{equation*}
\sum_{\tau \in \mathcal{T}}p_{\tau} + \sum_{k \in \mathcal{K}}p_k = P_{D},
\end{equation*}

\begin{equation*}
\sum_{j \in \mathcal{G}}g_j - \sum_{k \in \mathcal{K}}\phi _{k}p_k = G_{D},
\end{equation*}

\begin{equation*}
\;  \underline{p}_i \leq p_i \leq \overline{p}_i, \; 0 \leq g_j \leq \overline{g}_j,
\end{equation*}

where $q$ is the vector which contains: the conventional power generation $p_{\tau}$ (whose generators are in layer $\mathcal{T}$), the gas-fired generation $p_{k}$ (whose generators are in layer $\mathcal{K}$), and the gas supply $g_{j}$ (whose generators are in layer $\mathcal{G}$). $\mathcal{L}$ is the supra-Laplacian matrix of the three-layer multiplex network (See Fig. <ref>). $P_{D}$ is the electrical power demand, $G_{D}$ is the inelastic gas demand, and $\phi_k$ is the fuel conversion factor. $C_{\tau}^{[\mathcal{T}]}$, $C_{k}^{[\mathcal{K}]}$, and $C_{j}^{[\mathcal{G}]}$ are the cost functions of power production and gas suppliers in each layer.

Schematics for the three-layer micro-grid power generation and dispatch.

The fuel consumption of natural gas-fired generators is translated into a nodal time-varying gas demand $\phi_k p_k$ via the fuel consumption factor.

Therefore, the associated saddle-point flow yields

\begin{equation}
\dot{q}=-\frac{\partial C (q)}{\partial q}-\mathcal{L}q-\mathcal{L}\lambda -\mu \frac{\partial g(q)}{\partial q},
\end{equation}

\begin{equation}
\dot{\lambda}=\mathcal{L}q
\end{equation}

\begin{equation}
\dot{\mu}=g(q)
\end{equation}

\begin{equation*}
g (q) = \left(\begin{IEEEeqnarraybox*}[][c]{,c,}
\sum_{\tau \in \mathcal{T}}p_{\tau} + \sum_{k \in \mathcal{K}}p_k - P_{D}  \\
\sum_{j \in \mathcal{G}}g_j - \sum_{k \in \mathcal{K}}\phi _{k}p_k - G_{D}
\end{IEEEeqnarraybox*}\right) .
\end{equation*}

Let us assume a cost function defined by

\begin{equation*}
C(q)   =\begin{pmatrix}
\vdots \\
C_{\tau }^{[\mathcal{T}]}\\
\vdots \\
C_{\kappa }^{[ K ]}\\
\vdots \\
\vdots
\end{pmatrix} =\begin{pmatrix}
\vdots \\
\frac{1}{2} p_{\tau }^{2}\\
\vdots \\
\frac{1}{2} p_{\kappa }^{2}\\
\vdots \\
\frac{1}{2} g_{j}^{2}\\
\vdots
\end{pmatrix} .
\end{equation*}

For this situation, the supra-Laplacian matrix reads

\begin{equation*}
\mathcal{L}=\begin{pmatrix}
\mathcal{L}_{\mathcal{T}} & -D_{12} I_{7\times 7} & 0\\
-D_{12} I_{7\times 7} & \mathcal{L}_{K} & -D_{23} I_{7\times 7}\\
0 & -D_{23} I_{7\times 7} & \mathcal{L}_{\mathcal{G}}
\end{pmatrix}\\ ,
\end{equation*}

being $\mathcal{L}_{\mathcal{T}}$, $\mathcal{L}_{\mathcal{K}}$, and $\mathcal{L}_{\mathcal{G}}$ the Laplacian matrix for the layer $\mathcal{T}$, $\mathcal{K}$, and $\mathcal{G}$ respectively.

For the numerical implementation, we have that

\begin{align*}
-\frac{\partial C(q)}{\partial q} & =-\begin{pmatrix}
\vdots \\
p_{\tau }\\
\vdots \\
p_{\kappa }\\
\vdots \\
\vdots
\end{pmatrix} =-q .
\end{align*}

Therefore, the temporal evolution of the vector $q$ is given by (with $n=7$)

\begin{align*}
\ddot{q} & =-( I_{3n\times 3n} +\mathcal{L})\dot{q} -\mathcal{L}\dot{\lambda } -\dot{\mu }\frac{\partial g(q)}{\partial q}\\
& =-( I_{3n\times 3n} +\mathcal{L})\dot{q} -\mathcal{L}(\mathcal{L}q) -g(q)\frac{\partial g(q)}{\partial q}\\
& =-( I_{3n\times 3n} +\mathcal{L})\dot{q} -\mathcal{L}^{2}q -g(q)\frac{\partial g(q)}{\partial q} ,
\end{align*}

which could be written as

\begin{align}
\dot{q} & =v\\
\dot{v} & =-( I_{3n\times 3n} +\mathcal{L})v -\mathcal{L}^{2}q -g(q)\frac{\partial g(q)}{\partial q}.
\label{eq:eq2}
\end{align}

In Fig. <ref> we depict the consensus dynamics for the following parameters: $D_{12} = D_{23} = D_x = 0.6$, $D_{ \mathcal{T}} = D_{\mathcal{G}} = 0.2$, $D_{K}=0.8$, $\phi_k=0.7$, $P_D= G_{D}=100$. We can evidence how the different nodes in each layer reach the consensus by reaching values which are minimizing the total cost of operating power and natural gas microgrids. Notice that, in this situation, conventional generators are converging to the same value of gas supply. Gas-fired generators are going to elevate cost value. We can observe, how the totally connected graph, is going to the same value as the less connected structure.

Simulation for the coordinated dispatch for multienergy system with two-layers: energy and gas.

§ CONCLUSIONS

Currently, considering recent advances in network science, statistical mechanics, and big data, it has become evident that different networked systems are part of larger structures that define network behavior. Multilayer networks have emerged as an effective tool to model these situations and understand how networks interact with other networks. Multilayer systems exhibit various phenomena seen in traditional monolayer networks, including percolation, phase transitions, diffusion, and epidemic spreading. Despite this, there is a lack of results generalizing distributed control and optimization from monolayer networks to multilayer networks. In this paper, we have obtained two algorithms for distributed optimization in a particular case of multilayer networks with a one-by-one relation between the nodes in each layer. By leveraging the relationship between diffusion and consensus dynamics and employing a control approach, we developed a distributed primal-dual saddle point algorithm that exhibits fast convergence speed and robustness. We found that intra- and interlayer diffusion constants act as control parameters for optimal-consensus time. We observed that each layer achieves a local consensus before seeking a global consensus. Additionally, we discovered a critical phenomenon linked to the relationship between the consensus time and the interlayer diffusion constant. This finding complements previous work on determining diffusion dynamics for multilayer networks, which identified a phase transition related to the second eigenvalue of the supra-Laplacian matrix. Finally, we developed a distributed gradient descent algorithm for multilayer networks and observed that a time-varying positive gain plays an essential role in the convergence dynamics of the system.

Possible extensions to this work involve obtaining a diffusion operator for more general multilayer systems and applying the methodology presented here. Recent advancements in the tensor representation of multilayer networks have enabled the creation of such objects [39]. Additionally, recent research has developed Laplacian-like operators for hypergraphs and simplicial complexes [40, 41], which could be useful in developing distributed optimization algorithms for complex systems described using these general models [39]. Other related reaction-diffusion operators, such as the Dirac operator [42], could also be explored. Another potential extension involves using graph zeta functions to investigate diffusion dynamics in higher-order networks [43, 44]. These topics are currently being researched by the authors.

§ ACKNOWLEDGMENTS

§ ACKNOWLEDGMENT

This paper was partially supported by the VIII Convocatoria para el Desarrollo y Fortalecimiento de los Grupos de Investigación en Uniminuto with code C119-173, the Convocatoria de investigación para prototipado de tecnologías que promueven el cuidado o la restauración del medioambiente with code CPT123-200-5220, and Industrial Engineering Program from the Corporación Universitaria Minuto de Dios (Uniminuto, Colombia). We thank the Engineering and Physical Sciences Research Council (EPSRC) (Grants No. EP/R513143/1 and No. EP/T517793/1) for financial support.

[1]
A. Barabási, Network Science.1em plus 0.5em minus 0.4emCambridge University Press, Cambridge, 2016.

[2]
M. Newman, Networks: An Introduction, New York, NY, USA: Oxford
Univ.1em plus 0.5em minus 0.4emPress, 2010.

[3]
R. Albert and A.-L. Barabási, “Statistical mechanics of complex
networks,” Reviews of modern physics, vol. 74, no. 1, p. 47, 2002.

[4]
S. N. Dorogovtsev and J. F. Mendes, Evolution of networks: From
biological nets to the Internet and WWW.1em plus 0.5em minus
0.4emOxford university press, 2003.

[5]
G. Bianconi, Multilayer networks: structure and function.1em plus
0.5em minus 0.4emOxford university press, 2018.

[6]
M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A.
Porter, “Multilayer networks,” Journal of complex networks, vol. 2,
no. 3, pp. 203–271, 2014.

[7]
S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes,
M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, “The structure and
dynamics of multilayer networks,” Physics reports, vol. 544, no. 1,
pp. 1–122, 2014.

[8]
M. Vaiana and S. F. Muldoon, “Multilayer brain networks,” Journal of
Nonlinear Science, vol. 30, no. 5, pp. 2147–2169, 2020.

[9]
N. A. Kiani, D. Gomez-Cabrero, and G. Bianconi, Networks of Networks in
Biology: Concepts, Tools and Applications.1em plus 0.5em minus
0.4emCambridge University Press, 2021.

[10]
G. Mangioni, G. Jurman, and M. De Domenico, “Multilayer flows in molecular
networks identify biological modules in the human proteome,” IEEE
Transactions on Network Science and Engineering, vol. 7, no. 1, pp.
411–420, 2018.

[11]
S. Pilosof, M. A. Porter, M. Pascual, and S. Kéfi, “The multilayer nature
of ecological networks,” Nature Ecology & Evolution, vol. 1, no. 4,
p. 0101, 2017.

[12]
W.-B. Du, X.-L. Zhou, M. Jusup, and Z. Wang, “Physics of transportation:
Towards optimal capacity using the multilayer network framework,”
Scientific reports, vol. 6, no. 1, pp. 1–8, 2016.

[13]
V. Toro, E. Mojica-Nava, and N. Rakoto-Ravalontsalama, “Multiplex centrality
measurements applied to islanded microgrids,” International Journal of
Control, Automation and Systems, vol. 19, no. 1, pp. 449–458, 2021.

[14]
A. Santoro and V. Nicosia, “Optimal percolation in correlated multilayer
networks with overlap,” Physical Review Research, vol. 2, no. 3, p.
033122, 2020.

[15]
M. Turalska, K. Burghardt, M. Rohden, A. Swami, and R. M. D'Souza, “Cascading
failures in scale-free interdependent networks,” Physical Review E,
vol. 99, no. 3, p. 032308, 2019.

[16]
M. Perc, “Diffusion dynamics and information spreading in multilayer networks:
An overview,” The European Physical Journal Special Topics, vol. 228,
no. 11, pp. 2351–2355, 2019.

[17]
Y. Shang, “Resilient consensus for robust multiplex networks with asymmetric
confidence intervals,” IEEE Transactions on Network Science and
Engineering, vol. 8, no. 1, pp. 65–74, 2020.

[18]
Z. Guo, M. Lian, S. Wen, and T. Huang, “An adaptive multi-agent system with
duplex control laws for distributed resource allocation,” IEEE
Transactions on Network Science and Engineering, vol. 9, no. 2, pp.
389–400, 2021.

[19]
M. Salehi, R. Sharma, M. Marzolla, M. Magnani, P. Siyari, and D. Montesi,
“Spreading processes in multilayer networks,” IEEE Transactions on
Network Science and Engineering, vol. 2, no. 2, pp. 65–83, 2015.

[20]
T. Yang, X. Yi, J. Wu, Y. Yuan, D. Wu, Z. Meng, Y. Hong, H. Wang, Z. Lin, and
K. H. Johansson, “A survey of distributed optimization,” Annual
Reviews in Control, vol. 47, pp. 278–305, 2019.

[21]
A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multi-agent
optimization,” IEEE Transactions on Automatic Control, vol. 54,
no. 1, pp. 48–61, 2009.

[22]
A. Nedić and J. Liu, “Distributed optimization for control,” Annual
Review of Control, Robotics, and Autonomous Systems, vol. 1, pp. 77–103,

[23]
J. Wang and N. Elia, “A control perspective for centralized and distributed
convex optimization,” in 2011 50th IEEE conference on decision and
control and European control conference.1em plus 0.5em minus
0.4emIEEE, 2011, pp. 3800–3805.

[24]
P. Lin, W. Ren, and J. A. Farrell, “Distributed continuous-time optimization:
nonuniform gradient gains, finite-time convergence, and convex constraint
set,” IEEE Transactions on Automatic Control, vol. 62, no. 5, pp.
2239–2253, 2016.

[25]
A. Cherukuri, B. Gharesifard, and J. Cortes, “Saddle-point dynamics:
conditions for asymptotic stability of saddle points,” SIAM Journal on
Control and Optimization, vol. 55, no. 1, pp. 486–511, 2017.

[26]
D. Feijer and F. Paganini, “Stability of primal–dual gradient dynamics and
applications to network optimization,” Automatica, vol. 46, no. 12,
pp. 1974–1981, 2010.

[27]
M. Colombino, E. Dall’Anese, and A. Bernstein, “Online optimization as a
feedback controller: Stability and tracking,” IEEE Transactions on
Control of Network Systems, vol. 7, no. 1, pp. 422–432, 2019.

[28]
A. Hauswirth, S. Bolognani, G. Hug, and F. Dörfler, “Optimization
algorithms as robust feedback controllers,” arXiv preprint
arXiv:2103.11329, 2021.

[29]
D. Krishnamoorthy and S. Skogestad, “Real-time optimization as a feedback
control problem-a review,” Computers & Chemical Engineering, p.
107723, 2022.

[30]
J. Iacovacci and G. Bianconi, “Extracting information from multiplex
networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science,
vol. 26, no. 6, p. 065306, 2016.

[31]
P. Mancarella, “Mes (multi-energy systems): An overview of concepts and
evaluation models,” Energy, vol. 65, pp. 1–17, 2014.

[32]
E. Guelpa, A. Bischi, V. Verda, M. Chertkov, and H. Lund, “Towards future
infrastructures for sustainable multi-energy systems: A review,”
Energy, vol. 184, pp. 2–21, 2019.

[33]
M. Chertkov and G. Andersson, “Multienergy systems,” Proceedings of the
IEEE, vol. 108, no. 9, pp. 1387–1391, 2020.

[34]
N. Biggs, N. L. Biggs, and B. Norman, Algebraic graph theory.1em
plus 0.5em minus 0.4emCambridge university press, 1993, no. 67.

[35]
E. B. Curtis and J. A. Morrow, “The dirichlet to neumann map for a resistor
network,” SIAM Journal on Applied Mathematics, vol. 51, no. 4, pp.
1011–1029, 1991.

[36]
A. C. Gilbert, J. G. Hoskins, and J. C. Schotland, “Diffuse scattering on
graphs,” Linear Algebra and its Applications, vol. 496, pp. 1–35,

[37]
D. Bertsekas, Convex optimization theory.1em plus 0.5em minus
0.4emAthena Scientific, 2009, vol. 1.

[38]
S. Gomez, A. Diaz-Guilera, J. Gomez-Gardenes, C. J. Perez-Vicente, Y. Moreno,
and A. Arenas, “Diffusion dynamics on multiplex networks,” Physical
review letters, vol. 110, no. 2, p. 028701, 2013.

[39]
F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G.
Young, and G. Petri, “Networks beyond pairwise interactions: structure and
dynamics,” Physics Reports, vol. 874, pp. 1–92, 2020.

[40]
D. Horak and J. Jost, “Spectra of combinatorial laplace operators on
simplicial complexes,” Advances in Mathematics, vol. 244, pp.
303–336, 2013.

[41]
S. Majhi, M. Perc, and D. Ghosh, “Dynamics on higher-order networks: A
review,” Journal of the Royal Society Interface, vol. 19, no. 188, p.
20220043, 2022.

[42]
L. Calmon, M. T. Schaub, and G. Bianconi, “Dirac signal processing of
higher-order topological signals,” arXiv preprint arXiv:2301.10137,

[43]
G. Chinta, J. Jorgenson, and A. Karlsson, “Heat kernels on regular graphs and
generalized ihara zeta function formulas,” Monatshefte für
Mathematik, vol. 178, no. 2, pp. 171–190, 2015.

[44]
A. Saldivar, N. F. Svaiter, and C. A. Zarro, “Functional equations for
regularized zeta-functions and diffusion processes,” Journal of
Physics A: Mathematical and Theoretical, vol. 53, no. 23, p. 235205, 2020.

Christian David Rodríguez-Camargo
Received a B.S. degree in physics from Universidad Nacional de Colombia in 2014, an M.Sc. degree in physics from the Brazilian Center for Research in Physics (CBPF) in 2016. He does research in relativistic quantum information, quantum biology, quantum field theory in curved spaces, smart grits and statistical field theory in network science within researching collaboration at The Atomic, Molecular, Optical and Positron Physics (AMOPP) group of the Department of Physics and Astronomy from the University College London and Programa de Investigación sobre Adquisición y Análisis de Señales (PAAS-UN) from Universidad Nacional de Colombia.

Andrés F. Urquijo-Rodríguez
Received a B.S. and M.Sc. degree in physics from Universidad Nacional de Colombia, and he is currently pursuing his Ph.D. from Universidad Nacional de Colomba. He does research in network science, quantum information and condensed matter physics in low dimensional systems. He is an assistant professor at Corporación Universitaria Minuto de Dios, Bogotá, Colombia.

Eduardo Mojica-Nava
He received a B.S. degree in Electronics Engineering from Universidad Industrial de Santander in 2002, an M.Sc. degree in Electronics and Computer Science Engineering from Universidad de Los Andes, and a Ph.D. degree in Automatique et Informatique Industrielle from École des Mines de Nantes, Nantes, France and also Universidad de Los Andes in 2010. From 2011 to 2012, he was a Post-Doctoral Researcher at Universidad de Los Andes. He has been visiting professor at Université de Mons, Mons, Belgium, and Politecnico de Milano, Milan, Italy. Currently, he is a full professor at Universidad Nacional de Colombia, Bogotá, Colombia.