Automatica 159 (2024) 111309

Contents lists available at ScienceDirect

Automatica
journal homepage: www.elsevier.com/locate/automatica

Brief paper

Distributed continuous-time proximal algorithm for nonsmooth
resource allocation problem with coupled constraints✩
∗

Yi Huang a,b , Ziyang Meng c , Jian Sun a,b,d , , Gang Wang a,b,d
a

School of Automation, Beijing Institute of Technology, Beijing 100081, China
National Key Lab of Autonomous Intelligent Unmanned Systems, Beijing 100081, China
c
Department of Precision Instrument, Tsinghua University, Beijing 100084, China
d
Beijing Institute of Technology Chongqing Innovation Center, Chongqing 401120, China
b

article

info

Article history:
Received 10 November 2022
Received in revised form 9 July 2023
Accepted 23 August 2023
Available online 14 October 2023
Keywords:
Distributed optimization
Nonsmooth resource allocation
Proximal splitting
Coupled constraints

a b s t r a c t
This paper studies the distributed resource allocation problem with nonsmooth local cost functions
subject to the coupled equality and inequality constraints. In particular, each local cost function is
expressed as the sum of a differentiable function and two nonsmooth functions. By using the operator
splitting and primal–dual method, a continuous-time distributed proximal algorithm is developed,
which can be applied to more general local cost functions that are convex but not necessarily smooth.
In addition, the proposed algorithm is fully distributed in the sense that the gain parameter can be
determined locally and does not require any global information of the network. By applying Lyapunov
stability analysis and convex optimization theory, it is shown that the decision variables of all the
agents converge to an optimal solution. Finally, a simulation example is carried out to demonstrate
the effectiveness of the proposed algorithm.
© 2023 Published by Elsevier Ltd.

1. Introduction
During the last decades, distributed optimization over multiagent networks has received much attention due to its widespread
applications in sensor networks, machine learning, cooperative
localization, power systems and resource scheduling (Boyd, Parikh,
Chu, et al., 2011; Buehrer, Wymeersch, & Vaghefi, 2018; Chen,
Sun, & Wang, 2022; Liu, Zheng, & Yu, 2021; Wang, Giannakis, &
Chen, 2019; Wang, Zeng, Zhao, & Hong, 2023; Yang et al., 2019;
Yu, Liu, Zheng, & Zhu, 2021). In general, most distributed optimization problems in the literature can be categorized into two
types: the optimal consensus problem and resource allocation
problem. The main difference between these two problems is that
each agent of the first one involves a common decision variable,
while all the agents of the second one own the local decision
✩ This work has been supported in part by the National Natural Science
Foundation of China under Grants 62103223, 61925303, 62088101, 62273195,
62173034 and U19B2029. The material in this paper was not presented at
any conference. This paper was recommended for publication in revised form
by Associate Editor Giacomo Como under the direction of Editor Christos G.
Cassandras.
∗ Corresponding author at: School of Automation, Beijing Institute of
Technology, Beijing 100081, China.
E-mail addresses: yihuang@bit.edu.cn (Y. Huang),
ziyangmeng@mail.tsinghua.edu.cn (Z. Meng), sunjian@bit.edu.cn (J. Sun),
gangwang@bit.edu.cn (G. Wang).
https://doi.org/10.1016/j.automatica.2023.111309
0005-1098/© 2023 Published by Elsevier Ltd.

variables but they are coupled via a global equality constraint,
i.e, resource constraint.
Recently, there are a number of works concerning the resource
allocation problem over multi-agent systems. For instance, some
distributed resource allocation algorithms were proposed with an
assumption that the initial points of all the agents satisfy the coupled equality constraint (Chen & Yang, 2019; Cherukuri & Cortes,
2015; Li, Yu, Huang, & He, 2017). To remove the initialization
condition, some initialization-free distributed continuous-time
algorithms were proposed in Cherukuri and Cortes (2016), Yi,
Hong, and Liu (2016), Yun, Shim, and Ahn (2019). In addition,
besides the coupled equality constraint, the distributed resource
allocation problems subject to local constraints (e.g., local inequality/equality constraints and set constraints) were considered in Bai, Ye, Sun, and Hu (2017), Deng, Liang, and Hong (2017),
Deng, Nian, and Hu (2019), Li, Lu, and Huang (2018)
Nonsmooth optimization is a general problem in practical
applications, where the objective function contains some nonsmooth terms due to regularization constraint or penalty function
for some constraints. The subgradient-based method, which is
derived by replacing the gradient with an arbitrary subgradient, is an effective way for solving the nonsmooth optimization
problem. Chen, Yang, Song, and Lewis (2020), Deng et al. (2019),
Zeng, Yi, Hong, and Xie (2018) proposed some subgradient-based
distributed algorithms to solve the nonsmooth resource allocation problems. Due to the discontinuous nature, the subgradientbased distributed algorithms may result in violent vibrations

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

(iii) Compared with the nonsmooth analysis methods in Chen
et al. (2020), Deng et al. (2019), Liang et al. (2017), Zeng
et al. (2018), the convergence analysis of the proposed
algorithm is more intuitive since the proposed algorithms
are Lipschitz continuous and then Lyapunov stability theory can be applied. In addition, the proposed distributed
algorithm and convergence analysis can be easily extended
to the multiple nonsmooth optimization problem where
each local objective function contains multiple nonsmooth
convex functions.

of system states and cause difficulties in convergence analysis
since the traditional Lyapunov stability theory cannot be used.
Note that the proximal operator of a nonsmooth convex function is nonexpansive and Lipchitz continuous. Then, it provides
a powerful tool for designing a Lipchitz continuous algorithm
for nonsmooth optimization problem. Based on this attractive
property, some continuous-time distributed proximal algorithms
were developed in Wang, Chen, Zeng, and Xin (2020), Zhu, Wen,
Yu, and Yu (2021a, 2021b) to solve the nonsmooth resource
allocation problem. Also, some discrete-time distributed proximal
algorithms were proposed in Aybat, Wang, Lin, and Ma (2017),
Shi, Ling, Wu, and Yin (2015).
Note that the above distributed proximal algorithms cannot
be applied directly to the nonsmooth optimization problem when
two or multiple nonsmooth functions are contained in the local
objective function since the sum of two proximable nonsmooth
functions may be not proximable. Recently, Wei, Fang, Zeng, and
Chen (2019), Wei et al. (2022) proposed two distributed proximal
algorithms to solve two and multiple nonsmooth optimization
problems. However, the proposed algorithms rely on the strong
convexity of the local objective function. Moreover, the algorithms of Wei et al. (2019, 2022) are not fully distributed since
the selection criteria of gain parameters depend on the global
information of network, i.e., the eigenvalues of the Laplacian
matrix. In addition, besides the global resource constraint that
is a coupled linear equality constraint, the coupled inequality
constraint is also an important type of constraints in the resource allocation problem. In particular, it is more challenging
that the coupled inequality constraint has a nonlinear structure.
In Ge, Mei, Jiang, Qiu, and Yu (2022), Le, Chen, Yan, and Xi
(2017), Liang and Yin (2021), Liang, Zeng, and Hong (2017), some
continuous-time and discrete-time distributed algorithms were
proposed to solve the optimal resource allocation problems with
coupled inequality constraint. Note that the results of Le et al.
(2017), Liang and Yin (2021) are only applied to smooth local
objective function. Even though nonsmooth resource problems
were considered in Ge et al. (2022), Liang et al. (2017), the
proposed algorithms are based on the subgradient method and
are discontinuous. Moreover, the algorithm of Ge et al. (2022)
involves two-hop neighbors’ information.
Motivated by the above discussions, this paper aims to design
a continuous-time distributed proximal algorithm to solve the
nonsmooth resource allocation problem subject to the coupled
equality constraint and inequality constraints. Compared with the
existing results, the contributions of this paper are three-fold.

The rest of this paper is organized as follows. In Section 2,
some preliminaries are provided and the considered nonsmooth
resource allocation problem is formulated. In Section 3, a
continuous-time distributed proximal algorithm is developed and
the convergence analysis is provided. A simulation example is
given in Section 4 and some conclusions are drawn in Section 5.
2. Preliminaries and formulation
Notation. let Rn×n and Rn be the sets of n × n real matrices and
n dimension real vectors. Let In be an n × n identity matrix and 1n
be an n-dimensional vector with all entries being 1. ∥ · ∥ denotes
the Euclidean norm of a vector or a matrix. ⊗ represents the
operator of Kronecker product and λmin (A) denotes the smallest
eigenvalues of matrix A. Let col(x1 , . . . , xn ) be a column stack of
the vector xi , i = 1, 2, . . . , n, and diag(x1 , x2 , . . . , xn ) denotes the
diagonal matrix with diagonal entries being x1 to xn . Ker(M) and
Im(M) denote the kernel space and image space of matrix M,
respectively.
2.1. Graph theory
The communication network of n agents is described by an
undirected graph G = (V , E ), where V = {1, 2, . . . , n} and
E ⊆ V × V are the node set and the edge set, respectively. An
undirected edge (i, j) implies that the information flow between
nodes i and j is bidirectional. Ni = {j |(j, i) ∈ E } denotes the
neighbor set of agent i ∈ V . A weighted adjacency matrix of G is
denoted as A = [aij ] ∈ Rn×n , where aij > 0 if (j, i) ∈ E , and aij = 0
otherwise. Moreover, aij = aji holds for ∀(i, j) ∈ E of undirected
graph G . Here we assume that aii = 0 for ∀i ∈ V , which implies
that no loop exists. L = [Lij ] ∈ Rn×n is the Laplacian matrix of G ,
∑N
where Lii =
i=1 aij and Lij = −aij for i ̸ = j. For an undirected
graph G , L is a symmetric and positive semi-definite matrix, and
its eigenvalues are ordered as 0 ≤ λ2 (L) ≤ . . . ≤ λn (L).

(i) Unlike the resource allocation problems of Bai et al. (2017),
Deng et al. (2019), Li et al. (2018), Yi et al. (2016) where
only local constraints (e.g. local convex set or inequality constraint) were studied, the problem considered here
is more general and includes them as special cases. To
deal with this problem, we developed a new continuoustime distributed resource allocation algorithm by using the
operator splitting and primal–dual method.
(ii) In contrast to the distributed algorithms that require the
differentiability (Le et al., 2017; Yi et al., 2016) or strong
convexity (Wei et al., 2019, 2022; Zhu et al., 2021a, 2021b)
of the local cost functions, the proposed distributed algorithm can be applied to general local cost functions that are
only convex and not necessarily smooth. In addition, unlike
the optimization algorithms in Ge et al. (2022), Wei et al.
(2019, 2022), the proposed algorithm is fully distributed, in
which only one-hop neighbors’ information is involved and
the selection criterion of gain parameter is independent of
global information.

2.2. Proximal operators
Let f (x) be a lower semi-continuous convex function for x ∈
Rm . The proximal operator of f (x) at the point θ ∈ Rm is denoted
as
Proxf (θ ) = argminx∈Rm (f (x) +

1
2

∥x − θ∥2 ).

Let ∂ f be the subdifferential of f . If f (x) is convex, then ∂ f (x) is
monotone, i.e., (g1 − g2 )T (x1 − x2 ) ≥ 0, where g1 ∈ ∂ f (x1 ), g2 ∈
∂ f (x2 ). From the definition of proximal operator, one has that
θ − Proxf (θ ) ∈ ∂ f (Proxf (θ )). The proximal operator Proxf (·)
is nonexpansive, i.e., ∥Proxf (θ1 ) − Proxf (θ2 )∥ ≤ ∥θ1 − θ2 ∥ for
∀θ1 , θ2 ∈ Rm . A function f (x) is proximable if its proximal
operator has a closed or semi-closed form solution (Ryu & Yin,
2017; Wei et al., 2022).
2

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

2.3. Problem formulation

N
∑

Consider a group of n agents to cooperatively solve the following resource allocation problem.

i=1

min

F (x) =

N
∑

fi (xi )

hi (x∗i ) ≤ 0, µ∗0 ≥ 0, (µ∗0 )T

N
∑

hi (x∗i ) = 0.

(1b)

bi ,

(1c)

hi (xi ) ≤ 0,

(1d)

3.1. Distributed algorithm design

N
∑
i=1

i=1

∑N

To deal with the issue that the summation of fi1 (xi ) and fi2 (xi )
may not be proximable, inspired by the operator splitting technique in Davis and Yin (2017), we introduce an auxiliary variable
zi∗ ∈ Rni for each agent i ∈ V by splitting (2a) as

is the
where x = col(x1 , . . . , xN ) ∈ R , q =
i=1 ni , xi ∈ R
local decision variable of agent i ∈ V , fi (xi ) : Rni → R is the local
objective function composed of three components, in which fi0 (xi )
is a differentiable function, fi1 (xi ) and fi2 (xi ) are two nonsmooth
∑N
∑N
∑N
functions.
i=1 hi (xi ) ≤ 0
i=1 bi and h(x) =
i=1 Bi xi =
are the coupled equality and inequality constraints, where Bi ∈
Rp×ni , bi ∈ Rp are the local data, and hi (xi ) : Rni → Rm is the
local constraint function only known by agent i. Thereby, each
agent i ∈ V only has access to a partial information of the coupled
equality and inequality constraints.
q

ni

− ∇ fi0 (x∗i ) + BTi λ∗0
− ∇ hTi (x∗i )µ∗0 + γi zi∗ ∈ ∂ fi1 (x∗i ),

(3a)

− γi zi ∈ ∂ fi (xi ),
2

∗

∗

(3b)

where γi is a positive constant.
According to the definition of proximal operator, the above
Eq. (3) is transformed into

Remark 2.1. The optimization problem (1) is formulated in a
general form and includes many existing ones (e.g. Chen et al.
(2020), Yi et al. (2016), Zhu et al. (2021a)). For instance, if fi1 (xi ) is
an indicator function of convex set Ωi , then fi1 (xi ) is equivalent to
set constraint xi ∈ Ωi . In addition, if fi2 (x) = γ ∥xi ∥1 , the problem
(1) can represent the typical LASSO problem with set constraints.

x∗i = Proxf 1
i

(

x∗i − ∇ fi0 (x∗i )

+ BTi λ∗0 − ∇ hTi (x∗i )µ∗0

(4a)

)

+γi zi∗ ,

x∗i = Proxf 2 (x∗i − γi zi∗ ).

(4b)

i

It follows from (3) and Eq. (4) that −γi zi∗ is introduced by using
the operator splitting and its role is to estimate a subgradient of
fi2 (x∗i ).
To deal with the coupled equality and inequality constraints
in (1c) and (1d), similar to the results in Yi et al. (2016), we
first obtain the dual problem of (1) with respect to Lagrange
multipliers λc and µc

Assumption 2.1. The undirected graph G is connected.
Assumption 2.2. fi0 is a differentiable convex function. fi1 (xi ) and
fi2 (xi ) are nonsmooth lower semi-continuous convex functions
and are both proximable. hi (xi ) is differentiable and convex.
Assumption
2.3 (Slater’s
∑N
∑N Condition). There exists a point x̂i satisfying
B
x̂
=
i
i
i=1
i=1 bi and hi (x̂i ) < 0.

max

λc ∈Rp ,µc ∈Rm
+

Remark 2.2. Assumption 2.2 implies that each local function
fi (xi ) in the problem (1) is only convex but not strongly convex.
Assumption 2.3 is used to guarantee the existence of solution
to problem (1), which is widely common in the existing literatures (e.g.,Yi et al., 2016; Yun et al., 2019). Moreover, the
proximable condition in Assumption 2.2 can be satisfied for many
typical nonsmooth functions (e.g., indicator function of convex
set, l1 -norm function).

∑N

max
λ∈RNp ,µ∈RNm
+

Lemma 2.1.
Suppose that Assumptions 2.1–2.3 hold. x∗i is an optimal solution
to problem (1) if and only if there exist Lagrangian multipliers λ∗0 ∈
Rp , µ∗0 ∈ Rm such that the following conditions hold
0 ∈ ∇ fi0 (x∗i ) + ∂ fi1 (x∗i ) + ∂ fi2 (x∗i ),

−

λ0 + ∇
∗

hTi (x∗i )

g(λc , µc ) =

max

N
∑

λc ∈Rp ,µc ∈Rm
+

gi (λc , µc ),

(5)

i=1

where gi (λc , µc ) = infxi ∈Rni {fi (xi ) + λTc (bi − Bi xi ) + µTc hi (xi )}. It
follows from (5) that all the functions gi , i ∈ V share common
Lagrange multipliers λc and µc . To ensure that each function gi
owns local multipliers λi and µi , the consensus constraints (L ⊗
Ip )λ = 0 and (L1/2 ⊗ Ip )µ = 0, where L is the Laplacian matrix of
the graph G and L1/2 is its square root, λ = col(λ1 , . . . , λN ) ∈ RNp
and µ = col(µ1 , . . . , µN ) ∈ RNm
+ . Under Assumption 2.1, one has
that Ker(L) = Ker(L1/2 ) = Im(1N ) holds. Then, the dual problem
problem (5) can be rewritten as

The Lagrangian function of problem (1) is L1 (x, λ, µ) =
i=1
∑N
∑n
T
p
fi (xi ) +λTc
i=1 (bi − Bi xi ) +µc
i=1 hi (xi ), where λc ∈ R and µc ∈
Rm
+ are the Lagrange multipliers. Based on the Karush–Kuhn–
Tucker (KKT) conditions, the following lemma can be obtained.

BTi

(2c)

i=1

In this section, we first develop a continuous-time distributed
proximal algorithm such that the resource allocation problem can
be solved, and then the convergence analysis is provided.

Bi xi =
N
∑

(2b)

3. Main results

i=1

h(x) =

bi ,

i=1

(1a)

where fi (xi ) = fi0 (xi ) + fi1 (xi ) + fi2 (xi ),
s.t.

N
∑
i=1

N
∑

i=1

N
∑

Bi x∗i =

s.t.

g(λ, µ) =

N
∑

gi (λi , µi )

(6a)

i=1

(L ⊗ Ip )λ = 0, (L1/2 ⊗ Im )µ = 0

(6b)

The augmented Lagrangian duality of (6) with respect to multipliers y = col(y1 , . . . , yN ) ∈ RNp and s = col(s1 , . . . , sN ) ∈ RNm is
given by

(2a)

min

µ0 ,
∗

max

y∈RNp ,s∈RNm λ∈RNp ,µ∈RNm
+

3

L2 (λ, µ, y, z)

(7)

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

∑N

(iii) The dynamics of (µi , si ) is derived by using the projectionbased gradient decent ascent method for solving the max–
min problem (7), and its role is to guarantee the satisfaction
of the coupled inequality constraint (1d) and achieve consensus of µi for all i ∈ V . Unlike the dynamics ṡi =
∑
j∈Ni aij (µ̃i − µ̃j ) in Ge et al. (2022) that requires
∑ two-hop
neighbors’ information, the dynamics ṡi =
j∈Ni aij (µi −
µj ) in algorithm (1) only involves one-hop neighbors’ information.

1 T
T
where L2 (λ, µ, y, z) =
i=1 gi (λi , µi ) − y (L ⊗ Ip )λ− 2 λ (L ⊗ Ip )λ−
1 T
T 1/2
s (L ⊗ Im )µ − 2 µ (L ⊗ Im )µ.
According to the operator splitting in Eq. (4) and the augmented Lagrangian function L2 (λ, µ, y, z), we develop the following distributed proximal algorithm to solve the nonsmooth
resource allocation problem (1)

ẋi = Proxf 1 (xi − ∇ fi0 (xi ) + BTi λi

(8a)

i

−∇

hTi (xi )

µ̃i + γi zi + (1 + γi )żi ) − xi ,

żi = Proxf 2 (xi − γi zi ) − xi ,

(8b)

i

λ̇i = −(Bi xi + Bi ẋi − bi ) −

∑

aij (λi − λj )

Remark 3.1. Note from algorithm (8) that each agent only
has access to local gradient information ∇ fi0 (xi ), proximal operators Proxf 1 (·) and Proxf 2 (·), and local information Bi , bi , hi (xi )
i
i
in the coupled constraints (1c) and (1d). Each agent only needs
to receive three auxiliary variables (λj , yj , sj ) from its neighbors.
Moreover, unlike most algorithms that require the common gain
parameters of all the agents, the gain parameter γi of the algorithm (8) can be selected independently for each agent without
involving the global information of the network. Therefore, the
algorithm (8) is fully distributed and can be easily applied in the
large-scale networks.

(8c)

j∈Ni

−

∑

aij (yi − yj ),

j∈Ni

ẏi =

∑

aij (λi − λj ),

(8d)

j∈Ni

1

µ̇i = − (µi − µ̃i ),
∑2
ṡi =
aij (µi − µj ),

(8e)
(8f)

j∈Ni

3.2. Convergence analysis

where (xi , zi , λi , yi , µi , si ) ∈ Rni × Rni × Rp × Rp ×∑
Rm × Rm are
the local states of agent i ∈ V , µ̃i = (µi + hi (xi ) − j∈N aij (µi −
i
µj ) − si )+ , where (ζ )+ denotes a projection of any vector ζ ∈ Rm
+
on the set Rm
+ , i.e., (ζ ) = max(0, ζ ), and 0 < γi < 1. The initial
value si (0) is chosen to be zero, µi (0) ∈ Rm
+ and the other variables
can be initialized arbitrarily. The main steps of algorithm (8) are
shown as follows.

Let x = col(x1 , . . . , xN ), z = col(z1 , . . . , zN ), Proxf i (ξ ) =
col(Proxf i (ξ1 ), . . . , Proxf i (ξN )) for any vector ξ = col(ξ1 , . . . , ξN )
N
1
and i = 1, 2, ∇ h(x) = diag(∇ h1 (x1 ), . . . , ∇ hN (xN )), B =
diag(B1 , . . . , BN ), b = col(b1 , . . . , bN ), Γ = diag(γ1 In1 , . . . , γN InN ),
h(x) = col(h1 (x1 ), . . . , hN (xN )) and µ̃ = col(µ̃1 , . . ., µ̃N ) = (µ +
h(x) − (L ⊗ Im )µ − s)+ . Then, algorithm (8) can be written in the
following compact format

Algorithm 1 Continuous-time distributed proximal algorithm (1).

⎧
ẋ = Proxf 1 (x − ∇ f 0 (x) + BT λ − ∇ h(x)T µ̃
⎪
⎪
⎪
⎪
+ Γ z + (Iq + Γ )ż) − x
⎪
⎪
⎪
⎪
⎪
ż = Proxf 2 (x − Γ z) − x,
⎪
⎪
⎨
λ̇ = −(Bx + Bẋ − b)
⎪
− (L ⊗ Ip )λ − (L ⊗ Ip )y,
⎪
⎪
⎪
⎪
ẏ
=
(L ⊗ Ip )λ,
⎪
⎪
⎪
⎪
µ̇
=
− 21 µ + 12 (µ + h(x) − (L ⊗ Im )µ − s)+ ,
⎪
⎪
⎩
ṡ = (L ⊗ Im )µ.

1: Initiation:

2: xi (0), zi (0) ∈ R i , λi (0), yi (0) ∈ R , µi (0) ∈ R+ , si (0) = 0 for
n

p

m

∀i ∈ V ,
3: Update flows:
4: for Each agent i ∈ V (in parallel) do

Perform the update zi in (8b);
6: Perform the update xi in (8a);
7: Receive the neighbors’ variables λj , yj , µj , j ∈ Ni ;
8: Perform the updates λi , yi , µi , si in (8c)-(8f);
9: Broadcast the updated variables λi , yi , µi to j ∈ Ni .
10: end for
5:

(9)

Lemma 3.1.
Suppose that Assumptions 2.1–2.3 hold,
(x∗ , z ∗ , λ∗ , y∗ , µ∗ , s∗ ) is an equilibrium point of (9) if and only if
x∗ is an optimal solution to problem (1).

The algorithm design of (8) can be understood as follows.

Proof. Sufficiency: Let (x∗ , z ∗ , λ∗ , y∗ , µ∗ , s∗ ) be an equilibrium
point of (9), and then it follows that

(i) The role of xi is to estimate an optimal solution of problem
(1) and zi is to estimate a subgradient of − γ1 fi2 (xi ). With
i
the aid of zi , we obtain that the dynamics of xi is actually
equivalent to the gradient descent method that minimizes
the Lagrangian function L(x, λ, µ) with respect to xi . In
addition, the derivative feedback (1 + γi )żi is added in (8a)
to offset the cross term −x̄T ż brought by the derivative
of the Lyapunov function with respect to zi . Because the
cross term −x̄T ż can be neutralized, the strong convexity
of fi (xi ) required in Wei et al. (2022), Zhu et al. (2021a) can
be removed in our proposed algorithm (8).
(ii) The dynamics of (λi , yi ) is actually obtained by using the
gradient decent ascent method to solve the max–min problem (7), and its objective is to guarantee the satisfaction
of the coupled equality constraint (1c) and achieve the
consensus of λi . Here, we add the derivative feedback term
Bẋi to offset the cross term ẋT Bλ̄ brought by the derivative
of the Lyapunov function.

Proxf 1 (x∗ − ∇ f 0 (x∗ ) + BT λ∗

(10a)

− ∇ h (x )µ + Γ z ) = x ,
T

∗

∗

Proxf 2 (x − Γ z ∗ ) = x∗ ,
∗

∗

∗

(10b)

(Bx − b) + (L ⊗ Ip )y = 0,

(10c)

(L ⊗ Ip )λ∗ = 0,

(10d)

(µ∗ + h(x∗ ) − s∗ )+ = µ∗ ,

(10e)

(L ⊗ Im )µ = 0.

(10f)

∗

∗

∗

Under Assumption 2.1, it follows from (10d) and (10f) that λ∗ =
1N ⊗λ∗0 and µ∗ = 1N ⊗µ∗0 for some vectors λ∗0 ∈ Rp and µ∗0 ∈ Rm .
By left multiplying 1TN ⊗ Ip of (10c), one has that (2b) holds. Under
the initial condition si (0) = 0, ∀i ∈ V , it follows from the last
equation of (9) that (1TN ⊗ Im )s(t) = 0 for all t ≥ 0 and therefore
4

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

t, one has that s(t) − s(0) = (L1/2 ⊗ Im )ω(t) − (L1/2 ⊗ Im )ω(0).
Under the initial condition s(0) = (L1/2 ⊗ In )ω(0), we have that
s(t) = (L1/2 ⊗ Im )ω(t) holds for all t ≥ 0. Thus, (14) is obtained by
replacing s = (L1/2 ⊗ Im )ω in (9).
■

one further can derive that (1TN ⊗ Im )s∗ = 0. From (10e) and (10f),
one has that

µ∗ ≥ 0, h(x∗ ) − s∗ ≤ 0, (µ∗ )T (h(x∗ ) − s∗ ) = 0.

(11)

Substituting µ = 1N ⊗ µ0 into (11), we have that µ0 ≥ 0 and
∑N
T
∗
∗
(µ∗0 )T
i=1 hi (xi ) = 0 since (1N ⊗ Im )s = 0. By left
∑multiplying
N
T
1N ⊗ Im of the second inequality in (11), one has that i=1 hi (xi ) ≤
∗
0. This implies that there exists a vector µ0 such that (2c) can be
satisfied.
In addition, from (10a) and (10b), according to the properties
of proximal operator, we have that
∗

∗

− ∇ f (x∗ ) + BT (1N ⊗ λ∗0 )

∗

Remark 3.2. The initial value s(0) is set to be zero and it then
follows that (1TN ⊗ In )s(0) = 0. Under Assumption 2.1, one has
that Ker(1TN ⊗ Im ) = Im(L1/2 ⊗ Im ) and there exists ω(0) such that
s(0) = (L1/2 ⊗ Im )ω(0). This implies that the initial condition of
Proposition 3.1 is satisfied by setting s(0) to be zero.
Since algorithm (14) contains the terms (L1/2 ⊗ Im )ω and
(L1/2 ⊗ Im )µ involving the global information, it makes algorithm
(14) difficult to be implemented. Here, we should illustrate that
algorithm (14) is introduced to facilitate the convergence analysis
of (9). We next show the convergence analysis of (14). According
to the equivalency of algorithms (9) and (14), we also obtain the
convergence result of algorithm (9).
Let (x∗ , z ∗ , λ∗ , y∗ , µ∗ , ω∗ ) be an equilibrium point of (9). Define the error variables x̄ = x − x∗ , z̄ = z − z ∗ , ∑
λ̄ = λ − λ∗ , ȳ = y −
N
1
2
y∗ , µ̄ = µ − µ∗ , ω̄ = ω − ω∗ . Let φ (x, µ, s) =
i=1 fi (xi ) + 2 ∥µ̃∥ ,
1/2
+
where µ̃ = (µ + h(x) − (L ⊗ Im )µ − (L
⊗ Im )ω) . Define
the variable ξ = col(x̄, z̄ , λ̄, ȳ, µ̄, ω̄), and choose the Lyapunov
function V (ξ ) = V1 (x̄, z̄) + V2 (x̄, µ̄, ω̄) + V3 (λ̄, ȳ) + V4 (µ̄, ω̄), where

(12a)

− ∇ h(x∗ )T µ∗ + Γ z ∗ ∈ ∂ f 1 (x∗ ),
− Γ z ∗ ∈ ∂ f 2 (x∗ ).

(12b)

From (12a) and (12b), we can obtain that 0 ∈ −∇ f 0 (x∗ ) −
∂ f 1 (x∗ ) −∂ f 2 (x∗ ) + BT (1N ⊗λ∗0 ) −∇ h(x∗ )T µ∗ and one further derive
that 0 ∈ ∇ fi0 (x∗i ) + ∂ fi1 (x∗i ) + ∂ fi2 (x∗i ) − BTi λ∗0 + ∇ hTi (x∗i )µ∗0 . This
implies that (2a) holds. Based on Lemma 2.1, we conclude that x∗
is an optimal solution to problem (1).
Necessary: If x∗ is an optimal solution of problem (1), by using
Lemma 2.1, it follows that there exist λ∗0 , µ∗0 such that (2) can be
satisfied. Setting λ∗ = 1N ⊗ λ∗0 and µ∗ = 1N ⊗ µ∗0 , one has that
(10d) and (10f) hold. Note from (2b) that (1TN ⊗ Ip )(Bx∗ − b) = 0.
Since Ker(1TN ) = Im(L), there exists a constant vector y∗ ∈ RNp
such that (Bx∗ − b) + (L ⊗ Ip )y∗ = 0, i.e., (10c) holds. According
to (2c), we obtain that µ∗ ≥ 0 and (1TN ⊗ Im )h(x∗ ) ≤ 0. From
Ker(1TN ) = Im(L), one has that there exists a constant vector
χ ∗ ∈ RNm such that h(x∗ ) − (L ⊗ Im )χ ∗ ≤ 0. Also, we can verify
that (µ∗ )T (h(x∗ ) − (L ⊗ Im )χ ∗ ) = 0. By setting s∗ = (L ⊗ Im )χ ∗ , one
has that (10e) holds. In addition, it follows from (2a) that
0 ∈ −∇ f 0 (x∗ ) − ∂ f 1 (x∗ ) − ∂ f 2 (x∗ )

1 T
1
x̄ x̄ + z̄ T Γ z̄ − x̄T Γ z̄ ,
2
2
V2 (x̄, µ̄, ω̄) = φ (x, µ, s) − φ (x∗ , µ∗ , ω∗ )

V1 (x̄, z̄) =

− (x − x∗ )T (∇ f 0 (x∗ ) + ∇ hT (x∗ )µ∗ ) − (µ − µ∗ )T µ∗
+ (µ − µ∗ )T (L ⊗ Im )µ∗ + (ω − ω∗ )T (L1/2 ⊗ Im )µ∗ ,
1
1
∥λ̄∥2 + ∥ȳ∥2 ,
2
2
V4 (µ̄, ω̄) = V41 (µ̄) + V42 (µ̄, ω̄),
V3 (λ̄, ȳ) =

V41 (µ̄) =

(13)

2

∥µ̄∥2 , V42 (µ̄, ω̄) =

5
2

µ̄T (L ⊗ Im )µ̄

1

+ BT (1N ⊗ λ∗0 ) − ∇ h(x∗ )T µ∗ .

+ 2µ̄T (L1/2 ⊗ Im )ω̄ + ∥ω̄∥2 .
2

Then, there exists a constant vector −Γ z ∗ ∈ ∂ f 2 (x∗ ) such that
−∇ f 0 (x∗ ) + Γ z ∗ + BT (1N ⊗ λ∗0 ) − ∇ hT (x∗ )µ∗ ∈ ∂ f 1 (x∗ ) holds.
This implies that (10a) and (10b) can be satisfied. As a result, we
obtain that (x∗ , z ∗ , λ∗ , y∗ , µ∗ , s∗ ) is an equivalent point of (8). ■

Theorem 3.1. Suppose that Assumptions 2.1–2.3 hold and the gain
parameter satisfies 0 < γi < 1. Then, the proposed distributed
proximal algorithm (8) guarantees that the state x(t) converges to
an optimal solution of problem (1).

To facilitate the convergence analysis of the algorithm (9), an
equivalent algorithm is given as

⎧
ẋ = Proxf 1 (x − ∇ f 0 (x) + BT λ − ∇ h(x)T µ̃
⎪
⎪
⎪
⎪
+ Γ z + (Iq + Γ )ż) − x
⎪
⎪
⎪
⎪
⎪
ż
=
Proxf 2 (x − Γ z) − x,
⎪
⎪
⎪
⎪
λ̇
=
−(Bx + Bẋ − b)
⎨
− (L ⊗ Ip )λ − (L ⊗ Ip )y,
⎪
⎪
⎪ẏ = (L ⊗ Ip )λ,
⎪
⎪
⎪
⎪
µ̇ = − 21 µ + 12 (µ + h(x) − (L ⊗ Im )µ
⎪
⎪
⎪
⎪
− (L1/2 ⊗ Im )ω)+ ,
⎪
⎪
⎩
ω̇ = (L1/2 ⊗ Im )µ.

1

Proof. The proofs of Theorem 3.1 is divided into two steps.
Firstly, we prove that V (ξ ) is positive definite and radially unbounded. Secondly, we show that x(t) converges to an optimal
solution to problem (1).
Step 1: We show that V (ξ ) is positive definite and radially unbounded. From V1 (x̄, z̄), one has that

(14)

V1 (x̄, z̄) ≥ λmin (P1 )(∥x̄∥2 + ∥z̄ ∥2 ),

[
where P1 =

1/2Iq
−Γ /2

]
−Γ /2
is positive definite if 0 < γi < 1
Γ /2

for ∀i ∈ V .
Note that f 0 (x) is convex and µ̃ = (µ+ h(x) − (L ⊗ Im )µ− (L1/2 ⊗
Im )ω)+ is also convex with respect to (µ, x, ω), it then follows
that φ (x, µ, ω) is convex. In addition, one has that ∇x φ (x, µ, ω) =
∇ f 0 (x)+∇ hT (x)µ̃, ∇µ φ (x, µ, ω) = µ̃−(L⊗Im )µ̃ and ∇ω φ (x, µ, ω) =
−(L1/2 ⊗ Im )µ̃. According to the convex property of φ (x, µ, ω) and
µ̃∗ = µ∗ , it follows that

where ω ∈ RNm and the other variables are the same as those
in (9). Under Assumption 2.1, we know that L1/2 is still positive
semi-definite and symmetric.
Proposition 3.1. Algorithms (9) and (14) are equivalent if the initial
condition s(0) = (L1/2 ⊗ Im )ω(0) holds.

φ (x, µ, s) − φ (x∗ , µ∗ , s∗ )

1/2

Proof. Firstly, let s = (L
⊗ In )ω and it is obvious that (9) can
be derived from (14). Conversely, define ω̇ = (L1/2 ⊗ In )µ, and
it then follows that ṡ = (L1/2 ⊗ Im )ω̇. Integrating it from 0 to

− (x − x∗ )T (∇ f 0 (x∗ ) + ∇ hT (x∗ )µ∗ ) − (µ − µ∗ )T µ∗
+ (µ − µ∗ )T (L ⊗ Im )µ∗ + (ω − ω∗ )T (L1/2 ⊗ Im )µ∗ ≥ 0.
5

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

In addition, we note that

Note that (µ̃ − µ∗ + µ − µ∗ )T µ̇ = 12 (µ̃ + µ − 2µ∗ )T (µ̃ − µ) =
− 12 ∥µ − µ̃∥2 − (µ̃ − µ∗ )T (µ − µ̃). According to the property of
projection, one has that (µ + h(x) − (L ⊗ Im )µ − (L1/2 ⊗ Im )ω −
µ̃)T (µ̃ − µ∗ ) ≥ 0, which implies that

V4 (µ̄, ω̄) ≥ λmin (P2 )(∥µ̄∥2 + ∥ω̄∥2 ),
I + 5L 2L1/2
is positive definite. It then follows
where P2 =
2L1/2
I
2
that V (ξ ) ≥ λmin (P1 )(∥x̄∥ + ∥z̄ ∥2 ) + λmin (P2 )(∥µ̄∥2 + ∥ω̄∥2 ) +
1
∥λ̄∥2 + 12 ∥ȳ∥2 . Obviously, we know that V (ξ ) is positive definite
2
and radially unbounded. Also, we easily obtain that V (ξ ) = 0 if
and only if ξ = 0.
1
2

[

]

− (µ̃ − µ∗ )T (µ − µ̃)
≤ hT (x)(µ̃ − µ∗ ) − ωT (L1/2 ⊗ Im )(µ̃ − µ∗ )
− µT (L ⊗ Im )(µ̃ − µ∗ ).
It then follows that

Step 2: We show that x(t) converges to an optimal solution of
problem (1). Let Φ = ∇ f 0 (x) − BT λ + ∇ h(x)T µ̃ − γ z − (1 + γ )ż
and Φ ∗ = ∇ f 0 (x∗ ) − BT λ∗ + ∇ h(x∗ )T µ∗ − γ z ∗ . From (9) and (10),
one has that

− x̄T (∇ hT (x)µ̃ − ∇ hT (x∗ )µ∗ ) + (µ̃ − µ∗ + µ − µ∗ )T µ̇
1

= − ∥µ − µ̃∥2 + µ̃T (h(x) − ∇ h(x)x̄)

2
+ (µ∗ )T (−h(x) + ∇ hT (x∗ )x̄) − ωT (L1/2 ⊗ Im )(µ̃ − µ∗ )

ẋ + x = Proxf 1 (x − Φ ),

− µT (L ⊗ Im )(µ̃ − µ∗ )

x∗ = Proxf 1 (x∗ − Φ ∗ ),

1

ż + x = Proxf 2 (x − Γ z),

= − ∥µ − µ̃∥2 + µ̃T (h(x) − h(x∗ ) − ∇ h(x)x̄) + µ̃T h(x∗ )
2

x = Proxf 2 (x − Γ z ),
∗

∗

∗

+ (µ∗ )T (−h(x) + h(x∗ ) + ∇ hT (x∗ )x̄) − (µ∗ )T h(x∗ )
− ωT (L1/2 ⊗ Im )(µ̃ − µ∗ ) − µT (L ⊗ Im )(µ̃ − µ∗ )

From the definition of proximal operator, we have that

− Φ − ẋ ∈ ∂ f (ẋ + x), −Φ ∈ ∂ f (x ),
1

1

∗

∗

2

∗

1

(15a)

− Γ z − ż ∈ ∂ f (ż + x), −Γ z ∈ ∂ f (x ).
2

∗

≤ − ∥µ − µ̃∥2 + (µ̃ − µ∗ )T (h(x∗ ) − (L1/2 ⊗ Im )ω∗ )

2
+ (µ̃ − µ∗ )T (L1/2 ⊗ Im )ω∗ − ωT (L1/2 ⊗ Im )(µ̃ − µ∗ )

(15b)

By virtue of the convex property of fi1 (xi ) and fi2 (xi ), note from
(15) that
(−Φ − ẋ + Φ ∗ )T (ẋ + x̄) ≥ 0,

(16a)

(−Γ z̄ − ż)T (ż + x̄) ≥ 0.

(16b)

− µT (L ⊗ Im )(µ̃ − µ∗ )
1

≤ − ∥µ − µ̃∥2 − (µ̃ − µ∗ )T (L1/2 ⊗ Im )(ω − ω∗ )
2

− µT (L ⊗ Im )(µ̃ − µ∗ ),
where the first inequality is obtained by using the fact that h(x) −
h(x∗ )−∇ h(x)x̄ ≤ 0 and h(x)−h(x∗ )−∇ hT (x∗ )x̄ ≥ 0 according to the
convex property of h(x), and µ̃ ≥ 0, µ∗ ≥ 0. In addition, the last
inequality is derived by using (µ̃−µ∗ )T (h(x∗ ) − (L1/2 ⊗ Im )ω∗ ) ≤ 0.
Based on the above results, it follows that

With some calculations on (16), we obtain that
(x̄ + Φ − Φ ∗ )T ẋ ≤ −x̄T (Φ − Φ ∗ ) − ∥ẋ∥2 ,

(17a)

(x̄ + Γ z̄) ż ≤ −x̄ Γ z̄ − ∥ż ∥ .

(17b)

T

T

2

Taking the derivative of V1 (x̄, z̄), V2 (x̄, µ̄, s̄), V3 (λ̄, ȳ) and V4 (µ̄, s̄),
one has that

V̇1 + V̇2 + V̇41 ≤ −x̄T (∇ f 0 (x) − ∇ f 0 (x∗ )) − ∥ẋ∥2
1

+ λ̄T B(ẋ + x̄) − ∥ż ∥2 + ẋT (Iq + Γ )ż − ∥µ − µ̃∥2

V̇1 = (x̄ − Γ z̄)T ẋ + (z̄ − x̄)T Γ ż ,

2

∂ V2
∂ V2
∂ V2
ẋ +
µ̇ +
ω̇ = (∇ f 0 (x) − ∇ f 0 (x∗ )
∂x
∂µ
∂ω
+ ∇ hT (x)µ̃ − ∇ hT (x∗ )µ∗ )T ẋ + (µ̃ − µ∗ )T µ̇

− (µ̃ − µ∗ )T (L ⊗ Im )µ̇ − (µ̃ − µ∗ )T (L1/2 ⊗ Im )ω̇

V̇2 =

− (µ̃ − µ∗ )T (L1/2 ⊗ Im )ω̄ − (µ̃ − µ∗ )T (L ⊗ Im )µ,
According to the convexity property of f (x), one has that −x̄T
(∇ f (x) − ∇ f (x∗ )) ≤ 0. It then follows that

− (µ̃ − µ∗ )T (L ⊗ Im )µ̇ − (µ̃ − µ∗ )T (L ⊗ Im )ω̇,
V̇3 = λ̄T λ̇ + ȳT ẏ, V̇41 = µ̄T µ̇,

V̇1 + V̇2 + V̇41 ≤ −λmin (P3 )(∥ẋ∥2 + ∥ż ∥2 )

V̇42 = 5µ̄T (L ⊗ Im )µ̇ + 2µ̇T (L1/2 ⊗ Im )ω̄,

−

From (9), (10) and (17), V̇1 + V̇2 + V̇41 is calculated as
V̇1 + V̇2 + V̇41 = (x̄ − Γ z̄ + ∇ f 0 (x) − ∇ f 0 (x∗ )

+ ∇ hT (x)µ̃ − ∇ hT (x∗ )µ∗ )T ẋ + (z̄ − x̄)T Γ ż
− (µ̃ − µ ) (L

Υ0 = −(µ̃ − µ)T (L ⊗ Im )µ̇ − (µ − µ∗ )T (L ⊗ Im )µ̇

⊗ Im )ω̇

≤ −x̄T (Φ − Φ ∗ ) + λ̄T Bẋ + ż T (Iq + Γ )ẋ − ∥ẋ∥2 − x̄T Γ z̄

− (µ̃ − µ)T (L1/2 ⊗ Im )ω̇ − (µ − µ∗ )T (L1/2 ⊗ Im )ω̇

− ∥ż ∥ − x̄ (Iq + Γ )ż − (µ̃ − µ + µ − µ ) µ̇
2

T

2

(10), it follows that

+ (µ̃ − µ∗ + µ − µ∗ )T µ̇ − (µ̃ − µ∗ )T (L ⊗ Im )µ̇
1/2

∗

∗ T

− (µ̃ − µ)T (L1/2 ⊗ Im )ω̄ − (µ − µ∗ )T (L1/2 ⊗ Im )ω̄

− (µ̃ − µ∗ )T (L ⊗ Im )µ̇ − (µ̃ − µ∗ )T (L1/2 ⊗ Im )ω̇

− (µ̃ − µ)T (L ⊗ Im )µ − (µ − µ∗ )T (L ⊗ Im )µ

≤ −x̄T (∇ f 0 (x) − ∇ f 0 (x∗ )) − ∥ẋ∥2 + λ̄T B(ẋ + x̄)
− x̄T (∇ hT (x)µ̃ − ∇ hT (x∗ )µ∗ ) − ∥ż ∥2 + ż T (Iq + Γ )ẋ

= −2µ̇T (L ⊗ Im )µ̇ − µ̄T (L ⊗ Im )µ̇ − 2µ̇T (L ⊗ Im )µ

+ (µ̃ − µ∗ + µ − µ∗ )T µ̇ − (µ̃ − µ∗ )T (L ⊗ Im )µ̇

− µ̄T (L1/2 ⊗ Im )ω̇ − 2µ̇T (L1/2 ⊗ Im )ω̄ − ω̇T ω̄

1/2

− (µ̃ − µ ) (L
∗ T

(18)

∥µ − µ̃∥2 + λ̄T B(ẋ + x̄) + Υ0 ,
[
]
Iq
(Iq + Γ )/2
where P3 =
> 0 if 0 < γi < 1 for
(Iq + Γ )/2
Iq
∀i ∈ V , and Υ0 = −(µ̃ − µ∗ )T (L ⊗ Im )µ̇ − (µ̃ − µ∗ )T (L1/2 ⊗ Im )ω̇ −
(µ̃ − µ∗ )T (L1/2 ⊗ Im )ω̄ − (µ̃ − µ∗ )T (L ⊗ Im )µ. According to (9) and

+ 2µ̄T (L1/2 ⊗ Im )ω̇ + ω̄T ω̇.

∗ T

1

⊗ Im )ω̇.

− 2µ̇T (L ⊗ Im )µ − µ̄T (L1/2 ⊗ Im )ω̇
6

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

≤ −5µ̄T (L ⊗ Im )µ̇ − 2µ̄T (L1/2 ⊗ Im )ω̇

where fi (xi ) =

− 2µ̇T (L1/2 ⊗ Im )ω̄ − ω̇T ω̄

s.t
(19)

(22b)

B i xi =

N
∑

bi , h(x) =

i=1

N
∑

hi (xi ) ≤ 0,

(22c)

i=1

where the local objective function fi (xi ) of each agent i is composed of K + 1 functions fi0 (x), fi1 (x), . . . , fiK (x), in which fi0 (x) is a
j
differentiable function and fi (x) is a nonsmooth and proximable
function for ∀j = 1, . . . , K .
To solve the multiple nonsmooth resource allocation problem,
we develop the following distributed multi-proximal algorithm

V̇3 = −λ̄ (B(x̄ + ẋ) − (L ⊗ Ip )λ̄
T

− (L ⊗ Ip )ȳ) + ȳT (L ⊗ Ip )λ̄
(20)

By combining (18), (19) and (20), the derivative of V (ξ ) is simplified as
V̇ (ξ ) ≤ −λmin (P3 )(∥ẋ∥2 + ∥ż ∥2 )

N
∑
i=1

Similarly, by using (9) and (10), V̇3 is obtained as

≤ −λ̄T B(x̄ + ẋ) − λ̄T (L ⊗ Ip )λ̄.

fi (xi ),

j=0

Based on the derivative of V42 , we obtain that
V̇42 + Υ0 ≤ 0.

K
∑
j

ẋi = Proxf K
i

(

xi − ∇ fi0 (xi ) + BTi λi
K −1

(21)

− ∇ hTi (xi )µ̃i + γi

1

− ∥µ − µ̃∥2 − λ̄T (L ⊗ Ip )λ̄.
2

(23a)
K −1

∑ j

∑

j=1

j=1

zi + (1 + γi )

żi

)

−xi ,

j

j

żi = Proxf j (xi − γi zi ) − xi , j = 1, . . . , K − 1,

Together with the positiveness and unboundedness of V (ξ ), we
note from (21) that Ω1 = {ξ | V (ξ (t)) ≤ V (ξ (0))} is positively invariant. Then, we obtain that the equilibrium (x∗ , z ∗ , λ∗ , y∗ , µ∗ , ω∗ )
is Lyapunov stable and the trajectory ξ (t) is bounded for all t ≥ 0.
Define R = {(x̄, z̄ , λ̄, ȳ, µ̄, s̄) : V̇ (ξ ) = 0} ⊆ {(x̄, z̄, λ̄, ȳ, µ̄, ω̄) :
ẋ = 0, ż = 0, µ − µ̃ = 0, (L ⊗ Ip )λ̄ = 0}. Let M be the
largest invariant subset of R. Based on the invariance principle,
we have that (x̄(t), z̄(t), λ̄(t), ȳ(t), µ̄(t), ω̄(t)) → M when t →
∞, and the set M is positive invariant. This implies that if the
initial state (x̄(0), z̄(0), λ̄(0), ȳ(0), µ̄(0), ω̄(0)) belongs to the set
M, then (x̄(t), z̄(t), λ̄(t), ȳ(t), µ̄(t), ω̄(t)) ∈ M for all t ≥ 0. Then,
it follows that ẋ(t) ≡ 0 and ż(t) ≡ 0. From µ(t) − µ̃(t) ≡ 0 and
(L ⊗ Ip )λ̄(t) ≡ 0, one has that ẏ(t) ≡ 0 and µ̇(t) ≡ 0. In addition,
since x̄(t) ≡ x̄(0), ȳ(t) ≡ ȳ(0) and µ̄(t) ≡ µ̄(0), from the third and
sixth equations of (9), it follows that

(23b)

i

and the dynamics of the other states (λi , yi , µi , si ) are the same
as (8). In addition, to show the convergence of algorithm (23),
we modify the Lyapunov
V1 (x̄, z̄) of Theorem 3.1 as
∑K −function
1 1 j T
j
T
j
V1 (x̄, z̄) = 21 ∥x̄∥2 +
=
j=1 ( 2 (z̄ ) Γ z̄ − x̄ Γ z̄ ), where z̄
1
K −1
col(z̄ , . . . , z̄
). The convergence proof is almost similar to that
of Theorem 3.1 and is omitted here.
4. Numerical simulation
A simulation example is carried out to demonstrate the effectiveness of the proposed distributed proximal algorithm. The
algorithm (8) is applied to address the economic dispatch (ED)
problem in power systems, which is a typical class of resource
allocation problem. Consider a network of ten power generators.
Let Pi ∈ R be the output of the ith generator, and its cost function
fi (Pi ) is composed of a smooth function fi0 (Pi ) and two nonsmooth
functions fi1 (Pi ) and fi2 (Pi ). In particular, the smooth function
fi0 (Pi ) is described by fi0 (Pi ) = αi +βi Pi +ϖi Pi2 , where αi , βi , ϖi are
the ith generator’s parameters. The nonsmooth functions fi1 (Pi )
and fi2 (Pi ) are given by

λ̄˙ (t) = −Bx̄(0) − (L ⊗ Ip )ȳ(0) = C1 ,
ω̄˙ (t) = (L ⊗ Im )µ̄(0) = C2 ,
where C1 and C2 are constants. Then, if Ci ̸ = 0, i = 1, 2, one
has that λ̄(t) → ∞, s̄(t) → ∞ as t → ∞. This indicates
a contradiction to the boundedness of λ̄(t) and s̄(t). Thus, we
obtain that R ⊆ {(x̄, z̄ , λ̄, ȳ, µ̄, ω̄) | x̄˙ = 0, z̄˙ = 0, λ̄˙ = 0, ȳ˙ =
˙ = 0, ω̄˙ = 0}. This implies that any point in the set R is
0, µ̄
a Lyapunov stable equilibrium. Based on Theorem 4.4 of Khalil
(2002), it follows that the trajectory (x̄, z̄ , λ̄, ȳ, µ̄) converges to
a point in R. In addition, based on the results of Lemma 3.1,
we have that the state x(t) converges to an optimal solution to
problem (1).
■

fi1 (Pi ) =

if Pi ∈ Ωi 2
f (Pi ) = |Pi − ci |,
∞ if Pi ∈
/ Ωi , i

{

0

where Ωi = [Pimin , Pimax ] with Pimin , Pimax being the minimum and
maximum power outputs, and ci is a constant parameter. The
proximal operators of fi1 (xi ) and fi2 (xi ) are given by
Proxf 1 (Pi ) = PΩi (Pi ) = argminy∈Ωi |y − Pi |2 ,
i

Proxf 2 (Pi ) = [ψ1 , ψ2 ]T ,

Remark 3.3. From the proofs of Theorem 3.1, we note that the
derivative feedback items −(1 + γ )żi and −Bi ẋi in the dynamics
of xi and λi are very important to ensure the convergence of the
algorithm (1), which can be used to offset some cross terms of
V̇1 (x̄, z̄) + V̇2 (x̄, µ̄, ω̄) and V̇3 (λ̄, ȳ). In addition, by constructing
V4 (µ̄, ω̄), we obtain that V̇42 + Υ0 ≤ 0. This is a critical step
to ensure that the gain parameter of (8) can be selected without
requiring the global information of network.

i

where ψi for i = 1, 2 are

⎧
⎨Pi − 1 if Pi > ci + 1
ψi = ci
if ci − 1 ≤ Pi ≤ ci + 1
⎩
Pi + 1 if Pi < ci − 1.
∑10

In addition, we consider the coupled equality constraint i=1 Pi =
∑10
∑10
i Di and inequality constraint
i=1 C (Pi ) ≤ 0 in the ED
problem, where Di is the load demand of the ith generator and
1
C (Pi ) = 10
|Pi − πi,1 |2 − πi,2 .
We use a ring graph to describe the network topology G of
ten generators, i.e., 1 ⇔ 2 ⇔ 3 ⇔ . . . ⇔ 10 ⇔ 1, in
which aij = 1 for any edge (i, j) ∈ E . The parameter setting
of αi , βi , ϖi , ci , Di in these ten generators are given in Table 1.
The other parameters are chosen as Pimin = 0, Pimax = 40 −

Note that only two nonsmooth functions fi1 (xi ) and fi2 (xi ) are
considered in each local function of the optimization problem
(1). The distributed algorithm (8) can be also extended to solve
the multiple nonsmooth resource allocation problems, which is
formulated as
min F (x) =

N
∑

fi (xi )

(22a)

i=1

7

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309

Table 1
Parameters of ten generators.
Generator i

αi

βi

ϖi

ci

Di

1
2
3
4
5
6
7
8
9
10

10
15
18
19
10
26
11
20
23
14

3
7
8
9
10
5
4
6
2
4

2
4
1
1
2
1
0
0
0
0

20
11
19
10
17
18
20
22
15
10

10
20
20
15
12
14
20
10
22
13

Fig. 3. The trajectories of the variables λi , yi , µi and si for i ∈ V .

5. Conclusion
This paper develops a continuous-time distributed proximal
algorithm with Lipschitz continuity to solve the nonsmooth resource allocation problem with coupled equality and inequality
constraints. In contrast to the previous results, the proposed
distributed algorithm does not require the differentiability or
strong convexity of local cost function. Moreover, the proposed
algorithm is Lipschitz continuous and fully distributed, which can
be easily implemented in practical applications. Future work includes considering the distributed nonsmooth resource allocation
problem over general unbalanced directed graph.

Fig. 1. The decision variable Pi under the developed algorithm (8).

References
Aybat, N. S., Wang, Z., Lin, T., & Ma, S. (2017). Distributed linearized alternating direction method of multipliers for composite convex consensus
optimization. IEEE Transactions on Automatic Control, 63(1), 5–20.
Bai, L., Ye, M., Sun, C., & Hu, G. (2017). Distributed economic dispatch control
via saddle point dynamics and consensus algorithms. IEEE Transactions on
Control Systems Technology, 27(2), 898–905.
Boyd, S., Parikh, N., Chu, E., et al. (2011). Distributed optimization and statistical
learning via the alternating direction method of multipliers. Foundations and
Trends in Machine learning, 3(1), 1–122.
Buehrer, R. M., Wymeersch, H., & Vaghefi, R. M. (2018). Collaborative sensor
network localization: Algorithms and practical issues. Proceedings of the IEEE,
106(6), 1089–1114.
Chen, J., Sun, J., & Wang, G. (2022). From unmanned systems to autonomous
intelligent systems. Engineering, 12, 16–19.
Chen, G., & Yang, Q. (2019). Distributed constrained optimization for multi-agent
networks with nonsmooth objective functions. Systems & Control Letters, 124,
60–67.
Chen, G., Yang, Q., Song, Y., & Lewis, F. L. (2020). A distributed continuoustime algorithm for nonsmooth constrained optimization. IEEE Transactions
on Automatic Control, 65(11), 4914–4921.
Cherukuri, A., & Cortes, J. (2015). Distributed generator coordination for initialization and anytime optimization in economic dispatch. IEEE Transactions on
Control of Network Systems, 2(3), 226–237.
Cherukuri, A., & Cortes, J. (2016). Initialization-free distributed coordination
for economic dispatch under varying loads and generator commitment.
Automatica, 74, 183–193.
Davis, D., & Yin, W. (2017). A three-operator splitting scheme and its
optimization applications. Set-Valued and Variational Analysis, 25, 829–858.
Deng, Z., Liang, S., & Hong, Y. (2017). Distributed continuous-time algorithms for resource allocation problems over weight-balanced digraphs. IEEE
Transactions on Cybernetics, 48(11), 3116–3125.
Deng, Z., Nian, X., & Hu, C. (2019). Distributed algorithm design for nonsmooth resource allocation problems. IEEE Transactions on Cybernetics, 50(7),
3208–3217.
Ge, Y., Mei, X., Jiang, H., Qiu, J., & Yu, Z. (2022). A novel method for distributed optimization with globally coupled constraints based on multi-agent
systems. Neurocomputing, 487, 289–299.
Khalil, H. (2002). Nonlinear systems (3rd ed.). Upper Saddle River, NJ, USA:
Prentice-Hall.
Le, X., Chen, S., Yan, Z., & Xi, J. (2017). A neurodynamic approach to distributed optimization with globally coupled constraints. IEEE Transactions on
Cybernetics, 48(11), 3149–3158.

Fig. 2. The convergence results of coupled equality and inequality constraints.

i, πi,1 = 20, πi,2 = 20 for i ∈ V . According to the chosen
parameters in Table 1, we know that the cost function fi (Pi ), i =
7, . . . , 10 are only convex but not strongly convex. The initial
values xi (0), zi (0), λi (0), yi (0), µi (0), i ∈ V are randomly chosen in
the interval [0, 2] and the initial value si (0) is set to be zero. The
distributed proximal algorithm (8) is implemented with γi = 0.5
for i = 1, . . . , 5 and γi = 0.8 for i = 6, . . . , 10.
Fig. 1 shows that the outputs Pi of all the generators converge
to an optimal solution P ∗ = [4.559, 1.779, 6.618, 6.618, 2.809,
8.8118, 33, 32, 31, 30]T and
optimal value of global cost func∑the
10
tion is obtained as f (P) =
f
.934. Fig. 2 describes
i=1 i (Pi ) =
∑1215
10
that the coupled equality constraint
(P
i=1 i − Di ) converges to
∑10
zero and inequality constraint i=1 Ci (Pi ) converges to a negative
constant. This implies that the obtained optimal solution satisfies
these coupled constraints. Fig. 3 depicts the trajectories of the
variables λi , yi , µi and si , respectively. It can be seen from Fig. 3
that the dual variables λi , µi ultimately reach consensus, and µi
is always non-negative. The variables yi , si are both bounded and
finally converge to some constant vectors.
8

Y. Huang, Z. Meng, J. Sun et al.

Automatica 159 (2024) 111309
Ziyang Meng is currently an Associate professor with
the Department of Precision Instrument, Tsinghua University, China. He received his B.S. degree with honors
from Huazhong University of Science & Technology,
Wuhan, China, in 2006, and Ph.D. degree from Tsinghua
University, Beijing, China, in 2010. He was an exchange
Ph.D. student at Utah State University, Logan, USA from
2008 to 2009. Prior to joining Tsinghua University,
he held postdoc, researcher, and Humboldt research
fellow positions at, respectively, Shanghai Jiao Tong
University, Shanghai, China, KTH Royal Institute of
Technology, Stockholm, Sweden, and Technical University of Munich, Munich,
Germany from 2010–2015. His research interests include distributed control
and optimization, space science, and intelligent navigation technique. He serves
as associate editor for Systems & Control Letters and IET Control Theory &
Applications. He is a Senior Member of IEEE and Fellow of IET.

Li, H., Lu, Q., & Huang, T. (2018). Convergence analysis of a distributed optimization algorithm with a general unbalanced directed communication network.
IEEE Transactions on Network Science and Engineering, 6(3), 237–248.
Li, C., Yu, X., Huang, T., & He, X. (2017). Distributed optimal consensus over
resource allocation network and its application to dynamical economic
dispatch. IEEE Transactions on Neural Networks and Learning Systems, 29(6),
2407–2418.
Liang, S., & Yin, G. (2021). Distributed smooth convex optimization with coupled
constraints. IEEE Transactions on Automatic Control, 65(1), 347–353.
Liang, S., Zeng, X., & Hong, Y. (2017). Distributed nonsmooth optimization
with coupled inequality constraints via modified Lagrangian function. IEEE
Transactions on Automatic Control, 63(6), 1753–1759.
Liu, H., Zheng, W. X., & Yu, W. (2021). Continuous-time algorithm based on
finite-time consensus for distributed constrained convex optimization. IEEE
Transactions on Automatic Control, 67(5), 2552–2559.
Ryu, E. K., & Yin, W. (2017). Proximal-proximal-gradient method. arXiv preprint
arXiv:1708.06908.
Shi, W., Ling, Q., Wu, G., & Yin, W. (2015). A proximal gradient algorithm for
decentralized composite optimization. IEEE Transactions on Signal Processing,
63(22), 6013–6023.
Wang, Q., Chen, J., Zeng, X., & Xin, B. (2020). Distributed proximal-gradient
algorithms for nonsmooth convex optimization of second-order multiagent systems. International Journal of Robust and Nonlinear Control, 30(17),
7574–7592.
Wang, G., Giannakis, G. B., & Chen, J. (2019). Robust and scalable power system
state estimation via composite optimization. IEEE Transactions on Smart Grid,
10(6), 6137–6147.
Wang, Y., Zeng, X., Zhao, W., & Hong, Y. (2023). A zeroth-order algorithm for
distributed optimization with stochastic stripe observations. Science China.
Information Sciences, 66(9), Article 199202.
Wei, Y., Fang, H., Zeng, X., & Chen, J. (2019). A smooth double proximal primal–
dual algorithm for a class of distributed nonsmooth optimization problems.
IEEE Transactions on Automatic Control, 65(4), 1800–1806.
Wei, Y., Shang, C., Fang, H., Zeng, X., Dou, L., & Pardalos, P. (2022). Solving a class
of nonsmooth resource allocation problems with directed graphs through
distributed Lipschitz continuous multi-proximal algorithms. Automatica, 136,
Article 110071.
Yang, T., Yi, X., Wu, J., Yuan, Y., Wu, D., & Meng, Z. (2019). A survey of distributed
optimization. Annual Reviews in Control, 47, 278–305.
Yi, P., Hong, Y., & Liu, F. (2016). Initialization-free distributed algorithms for
optimal resource allocation with feasibility constraints and application to
economic dispatch of power systems. Automatica, 74, 259–269.
Yu, W., Liu, H., Zheng, W. X., & Zhu, Y. (2021). Distributed discrete-time
convex optimization with nonidentical local constraints over time-varying
unbalanced directed graphs. Automatica, 134, Article 109899.
Yun, H., Shim, H., & Ahn, H. S. (2019). Initialization-free privacy-guaranteed
distributed algorithm for economic dispatch problem. Automatica, 102,
86–93.
Zeng, X., Yi, P., Hong, Y., & Xie, L. (2018). Distributed continuous-time algorithms
for nonsmooth extended monotropic optimization problems. SIAM Journal on
Control and Optimization, 56(6), 3973–3993.
Zhu, Y., Wen, G., Yu, W., & Yu, X. (2021a). Continuous-time distributed
proximal gradient algorithms for nonsmooth resource allocation over general digraphs. IEEE Transactions on Network Science and Engineering, 8(2),
1733–1744.
Zhu, Y., Wen, G., Yu, W., & Yu, X. (2021b). Nonsmooth resource allocation of multiagent systems with disturbances: a proximal approach. IEEE Transactions on
Control of Network Systems, 8(3), 1454–1464.

Jian Sun received his B.S. degree from the Department
of Automation and Electric Engineering, Jilin Institute
of Technology, Changchun, China, in 2001, the M.Sc.
degree from the Changchun Institute of Optics, Fine
Mechanics and Physics, Chinese Academy of Sciences
(CAS), Changchun, China, in 2004, and the Ph.D. degree
from the Institute of Automation, CAS, Beijing, China,
in 2007. He was a Research Fellow with the Faculty
of Advanced Technology, University of Glamorgan, Pontypridd, U.K., from 2008 to 2009. He was a PostDoctoral
Research Fellow with the Beijing Institute of Technology, Beijing, from 2007 to 2010. In 2010, he joined the School of Automation,
Beijing Institute of Technology, where he has been a Professor since 2013.
His current research interests include networked control systems, time-delay
systems, and security of cyber–physical systems. Dr. Sun is an Editorial Board
Member of the IEEE Transactions on Systems, Man and Cybernetics: Systems,
the Journal of Systems Science & Complexity, and Acta. Automatica Sinica.

Gang Wang (Senior Member, IEEE) received a B.Eng.
degree in Automatic Control in 2011, and a Ph.D.
degree in Control Science and Engineering in 2018,
both from the Beijing Institute of Technology, Beijing,
China. He also received a Ph.D. degree in Electrical
and Computer Engineering from the University of Minnesota, Minneapolis, USA, in 2018, where he stayed
as a postdoctoral researcher until July 2020. Since
August 2020, he has been a professor with the School
of Automation at the Beijing Institute of Technology.
His research interests focus on the areas of signal
processing, control, and reinforcement learning with applications to cyber–
physical systems and multi-agent systems. He was the recipient of the Best Paper
Award from the Frontiers of Information Technology & Electronic Engineering
(FITEE) in 2021, the Excellent Doctoral Dissertation Award from the Chinese
Association of Automation in 2019, the Best Student Paper Award from the
2017 European Signal Processing Conference, and the Best Conference Paper
at the 2019 IEEE Power & Energy Society General Meeting. He is currently an
Associate Editor of Signal Processing and the IEEE Transactions on Signal and
Information Processing over Networks, the IEEE Open Journal of Control Systems,
as well as an Early Career Advisory Board member of the IEEE/CAA Journal of
Automatica Sinica.

Yi Huang received the Ph.D. degree in Applied Mathematics from Beihang University, Beijing, China, in 2019.
From 2019 to 2021, he was a Post-Doctoral Research
Fellow with the Department of Precision Instrument,
Tsinghua University, Beijing, China. He is currently an
Assistant Professor with the School of Automation, Beijing Institute of Technology, Beijing, Beijing, China. His
research interests include attitude coordinated control,
distributed control and optimization.

9

