DRHC synthesis for simultaneous tracking and formation of nonhomogeneous multi-agents with time-varying communication topology

Abstract

The synthesis approach of distributed receding horizon control is studied for the simultaneous tracking and formation problem of nonhomogeneous multi-agents. Different from the existing works, the communication topology between multi-agents is allowed to be time-varying in this article, which meets miscellaneous conditions in practice. To accommodate the time-varying communication topology, we refresh at each sampling instant the individual cost function for each agent, according to the real-time neighbourhood. Moreover, to guarantee the exponential stability of the overall closed-loop system, we design an auxiliary constraint and impose it in the individual optimization problem. The recursive feasibility of the auxiliary constraint can be guaranteed by updating the formation weighting scalars in real time. By solving the individual optimization problem (with respect to the input, state and auxiliary constraints) at each sampling instant, each agent can obtain its optimal control input sequence. The implementation of the first control input among the sequence for each agent can steer the overall multi-agent system to cooperatively achieve the desired tracking and formation objective. The effectiveness and practicability of our results are demonstrated through the illustrative examples.

Keywords

Distributed receding horizon control simultaneous tracking and formation time-varying communication topology recursive feasibility closed-loop stability linear matrix inequalities

Introduction

Due to the advantages in handling the physical constraints and optimizing the control performance, the receding horizon control (RHC), also called model predictive control, has become more and more popular in the control academia. The surveys of RHC can be found in the studies by Garcia et al.¹ and Mayne et al.² In RHC, the controller solves a constrained optimization problem (with respect to the cost function and the physical constraints associated with the future evolutions) at each sampling instant to obtain the optimal input sequence but implements only the first element among this sequence. However, the computation burden of RHC handling in the centralized form is usually tremendous, even intolerable, for the systems with high dimensional state variables (see the multi-agent systems and the large-scale industry systems). An efficient paradigm to tackle this issue is partitioning the overall system into a set of subsystems and applying a distributed RHC (DRHC) strategy for each subsystem.^3

–6

For the multi-agent systems, more and more attentions have been attracted on many important fields of research, such as the flocking problem,^7
–9 the consensus problem^10
–12 and the tracking and/or formation problem,^{5,6,13

–16} and many approaches have been proposed for the DRHC of these systems.^{3,5,6,13,17
–19} In this article, we focus on the simultaneous tracking and formation problem of these systems and the DRHC synthesis approach for this problem. In DRHC, the optimization problems of agents can be solved synchronously to save the time for all agents obtaining their control inputs^5,6; however, each agent has to use the assumed future evolutions of its neighbours (usually taken as the predicted evolutions solved at the last instant^5,6,13,18). As a result, the traditional stability paradigm of centralized strategy^2,20,21 is no longer suitable. To achieve the closed-loop stability guarantee in DRHC, some additional conditions or constraints are usually required for each agent, see for example, the sufficient stability condition derived in the study by Keviczky et al.,¹⁷ the additional contractive constraint employed in the study by Camponogara et al.³ and the compatibility constraints proposed in the literature.^13,18 Especially in the study by Wang and Ding,^5,6 the compatibility constraints are derived and imposed not only to guarantee the closed-loop stability but also to contribute to the collision avoidance.

However, all the above cited works study the DRHC for multi-agent systems based on the fixed communication topology, which fail when any communication link disconnects due to out of distance limitation, data packet dropout and so on. In practice, it is more reasonable to consider the time-varying communication topology due to miscellaneous reasons. Taking the time-varying topology into account, most of the researches focus on the consensus problem.^{10,12,14,22,23} As reported in the study by Ferrari-Trecate et al.,¹⁰ because of the receding horizon technique, the value of the consensus point depends not only on the initial conditions of the system but also on the sequence of agents’ states and on the communication topology along time. Thus, for the simultaneous tracking and formation problem, in which the objective to be achieved is pre-specified, the DRHC with respect to time-varying communication topology is still an issue to be solved.

Aiming at the simultaneous tracking and formation problem of nonhomogeneous multi-agent systems with time-varying communication topology, we contribute in this article to present a practicable DRHC synthesis approach by (i) online generating the real-time individual cost function for each agent, according to the pre-specified objective, the time-varying communication topology and the real-time updated formation weighing scalar; (ii) deriving an individual condition for each agent, collection of which can guarantee the closed-loop stability of the overall multi-agent system; and (iii) handling the derived condition as an auxiliary constraint and imposing it in the individual optimization problem with recursive feasibility guarantee. Following the imposed auxiliary constraint, the online generated cost function is substituted by its upper bound to be minimized, which endows the auxiliary constraint with the role of enhancing tracking and formation.

The following of this article is organized as follows. The various symbols used in this article are clarified in Table 1. In ‘System description and cost functions’ section, the systems of the multi-agent and reference object, the description of the simultaneous tracking and formation objective, the communication condition and the definition of neighbourhood are introduced. Then the individual cost function of each agent is online generated accordingly. In ‘Stability condition for the multi-agent system with primary distributed RHC’ section, the sufficient closed-loop stability condition of the overall multi-agent system is derived, with respect to the time-varying communication topology. ‘DRHC synthesis’ section then presents the DRHC synthesis approach by distributing the sufficient stability condition into individual ones, handling the individual condition as an auxiliary constraint with recursive feasibility, substituting the original cost function by its upper bound and providing the DRHC implementation algorithm. The properties of the proposed DRHC algorithm is analysed in ‘Property analysis’ section. Additionally, for easy application of the presented DRHC, ‘The optimization problem in LMIs form’ section provides a procedure to transform the individual optimization problem into linear matrix inequalities (LMIs) form, which completes the presentation of the DRHC approach. In ‘Illustrative examples’ section, the numerical and simulation examples are provided to demonstrate the effectiveness and practicability of our results. The ‘Conclusion’ section ends this article.

Table 1.

Notations.

Symbols	Meaning
$N_{a}$	Number of agents ( $N_{a} \geq 1$ )
ℕ_a	Index set of agents, ${1, 2, \dots, N_{a}}$
\|ℕ_a\|	Cardinality of the set ℕ_a
ℕ, ℕ⁻, ℕ⁺	${0, 1, \dots, N - 1}, {0, 1, \dots, N - 2}, {1, 2, \dots, N - 1}$
ℝⁿ	n-dimensional Euclidean space
$ℝ^{m \times n}$	m × n-dimensional real matrix set
I_n	n-order identity matrix
$diag {\dots}$	A block-diagonal matrix
⋆	Symmetric term with appropriate dimension in a block matrix
P > 0	Positive-definite matrix P
$\bar{λ} (P), \underline{λ} (P)$	Maximum, minimum eigenvalues of matrix P
$∥ x ∥, ∥ x ∥_{P}$	Two-norm of x, P-weight two-norm of x (i.e. $∥ x ∥ = \sqrt{x^{T} x}, ∥ x ∥_{P} = \sqrt{x^{T} P x}$ )
$x (t \| k)$	Future value of $x (k + t)$ , predicted at the instant k
$u^{*}, u^{\circ}$	Values of u corresponding to the optimal, feasible solutions of optimization
$[v; w]$	${[v^{T}, w^{T}]}^{T}$

System description and cost functions

Consider a multi-agent system with $N_{a}$ agents, which has the following discrete-time dynamic equations

x_{i} (k + 1) = A_{i} x_{i} (k) + B_{i} u_{i} (k), i \in ℕ_{a}, k \geq 0

where $x_{i} (k) = [p_{i} (k); v_{i} (k)] \in ℝ^{2 n_{i}} and u_{i} (k) \in ℝ^{m_{i}}$ denote the measurable state vector and control input of agent i, respectively; $p_{i} (k) \in ℝ^{n_{i}} and v_{i} (k) \in ℝ^{n_{i}}$ denote the position and velocity of agent i in n_i-D space, respectively. Usually, the state and input of each agent i are required to satisfy

x_{i} (k) \in X_{i} = {x_{i} \in ℝ^{2 n_{i}} | ∥ Ψ_{i, s} x_{i} ∥ \leq {\bar{ψ}}_{i, s}, s = 1, 2, \dots, c_{i}^{x}}, k \geq 0

u_{i} (k) \in U_{i} = {u_{i} \in ℝ^{m_{i}} | ∥ Φ_{i, s} u_{i} ∥ \leq {\bar{φ}}_{i, s}, s = 1, 2, \dots, c_{i}^{u}}, k \geq 0

where the bounded convex set $X_{i} (U_{i})$ denotes the admissible state (input) set, $Ψ_{i, s} \in ℝ^{n_{i, s}^{x} \times 2 n_{i}} with n_{i, s}^{x} rows and {\bar{ψ}}_{i, s} > 0 (Φ_{i, s} \in ℝ^{m_{i, s}^{u} \times m_{i}} with m_{i, s}^{u} rows and {\bar{φ}}_{i, s} > 0)$ are the state (input) constraint parameters and $c_{i}^{x} (c_{i}^{u})$ is the number of the state (input) constraints, of agent i. The reference agent to be tracked has the dynamics

x_{r} (k + 1) = A_{r} x_{r} (k) + B_{r} u_{r} (k), k \geq 0

where $x_{r} (k) = [p_{r} (k); v_{r} (k)] \in ℝ^{2 n_{r}} and u_{r} (k) \in ℝ^{m_{r}}$ are the reference state vector and control input, respectively. In this article, the dynamics parameters ${A_{i}, B_{i}} and {A_{r}, B_{r}}$ are allowed to be different, which accommodates the nonhomogeneous multi-agents.

The control objective of this article is to steer the multi-agents to simultaneously track the reference agent and achieve a pre-specified formation, without violating the required constraints (2) and (3). According to a simultaneous tracking and formation objective desired for the multi-agents, we denote $d_{i r}^{p} \in ℝ^{n_{i}}$ as the desired tracking position vector between agent i and reference agent, $d_{i j}^{p} \in ℝ^{n_{i}}$ as the desired formation position vector between the agents i and j, $d_{i r}^{x} = [d_{i r}^{p}; 0] \in ℝ^{2 n_{i}} and d_{i j}^{x} = [d_{i j}^{p}; 0] \in ℝ^{2 n_{i}}$ as the associated state vectors, $F = {(i, j) \in ℕ_{a} \times ℕ_{a} | d_{i j}^{p} is desired}$ as the desired formation pair index set and $ℕ_{i}^{0} = {j \in ℕ_{a} | (i, j) \in F}$ as the desired neighbour index set ( $j \in ℕ_{i}^{0} \Leftrightarrow i \in ℕ_{j}^{0}$ ). Because the multi-agents and the reference agent are allowed to be nonhomogeneous, even the dimensions of their states (inputs) can be not identical (i.e. $n_{i} \neq n_{j} \neq n_{r} and/or m_{i} \neq m_{j} \neq m_{r}$ for i ≠ j), we then define the coordinate transformation functions $F_{i r}^{x} (\cdot) : ℝ^{2 n_{r}} \to ℝ^{2 n_{i}}, F_{i r}^{u} (\cdot) : ℝ^{m_{r}} \to ℝ^{m_{i}} and F_{i j}^{x} (\cdot) : ℝ^{2 n_{j}} \to ℝ^{2 n_{i}}$ and formulate the control objective in a specific form as

Tracking : lim_{k \to \infty} [x_{i} (k) - F_{i r}^{x} (x_{r} (k))] = d_{i r}^{x} and lim_{k \to \infty} [u_{i} (k) - F_{i r}^{u} (u_{r} (k))] = 0, \forall i \in ℕ_{a}

Formation : lim_{k \to \infty} [x_{i} (k) - F_{i j}^{x} (x_{j} (k))] = d_{i j}^{x}, \forall j \in ℕ_{i}^{0}

To guarantee the desired control objective being realizable for the multi-agents, the following assumption is required.

Assumption 1

The multi-agents and reference agent are compatible, and the simultaneous tracking and formation objective is achievable, that is, the coordinate transformation functions and the desired vectors satisfy

Compatibility: F_{i r}^{x} (x_{r} (k + 1)) + d_{i r}^{x} = A_{i} (F_{i r}^{x} (x_{r} (k)) + d_{i r}^{x}) + B_{i} F_{i r}^{u} (u_{r} (k)), \forall i \in ℕ_{a}, k \geq 0

Trackability: F_{i r}^{x} (x_{r} (k)) + d_{i r}^{x} \in X_{i} and F_{i r}^{u} (u_{r} (k)) \in U_{i}, \forall i \in ℕ_{a}, k \geq 0

Consistency: d_{i j}^{x} = d_{i r}^{x} - F_{i j}^{x} (d_{j r}^{x}) = d_{i l}^{x} - F_{i j}^{x} (d_{j l}^{x}), \forall (i, j), (i, l), (j, l) \in F

Additionally, the function $F_{i j}^{x} (\cdot)$ is Lipschitz, that is, there exists a constant ${\bar{F}}_{i j}^{x} > 0 such that ∥ F_{i j}^{x} (x) ∥ \leq {\bar{F}}_{i j}^{x} ∥ x ∥, \forall x \in ℝ^{n_{j}}, (i, j) \in F$ .

Due to miscellaneous practical conditions, such as the distance limitation of effective communication, communication delays, data packet dropouts and so on, the communication topology is usually time-varying. Thus, the neighbourhood between agents needs a real-time presentation, which is defined as follows.

Definition 1

At each instant k, the agent j is a real-time neighbour of agent i if it is a desired neighbour of agent i and agent i can receive data from agent j. Hence, the real-time neighbour index set of agent i can be formulated as $ℕ_{i} (k) = {j \in ℕ_{i}^{0} | agent i can receive data from agent j at the k th instant}$ .

Note that $j \in ℕ_{i} (k) ⇎ i \in ℕ_{j} (k)$ due to different conditions. According to the above control objective and time-varying neighbourhood, we generate the ideal real-time individual cost function for each agent i as

J_{i} (x_{i r}, x_{i ▪}, u_{i r}, k) = \underset{tracking cost}{\underset{︸}{\sum_{t \in ℕ} (α_{i} ∥ x_{i r} (t | k) ∥^{2} + ρ_{i} ∥ u_{i r} (t | k) ∥^{2})}} + \underset{formation cost}{\underset{︸}{\sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k)} β_{i} ∥ x_{i j} (t | k) ∥^{2}}} + \underset{terminal cost}{\underset{︸}{∥ x_{i r} (N | k) ∥_{P_{i} (k)}^{2}}}

where N is the predictive horizon; $x_{i r} (t | k) = x_{i} (t | k) - F_{i r}^{x} (x_{r} (k + t)) - d_{i r}^{x}, x_{i j} (t | k) = x_{i} (t | k) - F_{i j}^{x} (x_{j} (t | k)) - d_{i j}^{x} = x_{i r} (t | k) - F_{i j}^{x} (x_{j r} (t | k)) and u_{i r} (t | k) = u_{i} (t | k) - F_{i r}^{u} (u_{r} (k + t))$ denote the predictive deviations; $α_{i} > 0 and ρ_{i} > 0$ are the pre-specified tracking weighting scalars; $β_{i} > 0$ is the pre-specified formation weighting scalar; $P_{i} (k) > 0$ is the terminal weighting matrix determined in real time.

For brief expression, the cost function (10) can be reformulated as

J_{i} (x_{i r}, x_{i ▪}, u_{i r}, k) = \sum_{t \in ℕ} (∥ z_{i r} (t | k) ∥_{Q_{i} (k)}^{2} + ∥ u_{i r} (t | k) ∥_{R_{i}}^{2}) + ∥ x_{i r} (N | k) ∥_{P_{i} (k)}^{2}

where $z_{i r} = [x_{j_{1}, r}; \dots; x_{j_{| ℕ_{i}^{0} |}, r}; x_{i r}] \in ℝ^{2 n_{i}}$ is the local collected vector, with ${j_{1}, \dots, j_{| ℕ_{i}^{0} |}} = ℕ_{i}^{0} and n_{i} = \sum_{j \in ℕ_{i}^{0}} n_{j} + n_{i}; Q_{i} (k) > 0 and R_{i} > 0$ are the local weighting matrices calculated according to α_i, ρ_i and β_i as in proposition 1.

Proposition 1

$R_{i} = ρ_{i} I_{m_{i}} and Q_{i} (k) = G_{i}^{T} W_{i} (k) G_{i}$ with

G_{i} = [\begin{matrix} - F_{i, j_{1}}^{x} (\cdot) & 0 & \dots & 0 & I_{2 n_{i}} \\ 0 & - F_{i, j_{2}}^{x} (\cdot) & ⋱ & ⋮ & I_{2 n_{i}} \\ ⋮ & ⋱ & ⋱ & 0 & ⋮ \\ 0 & \dots & 0 & - F_{i, j_{| ℕ_{i}^{0} |}}^{x} (\cdot) & I_{2 n_{i}} \\ 0 & \dots & \dots & 0 & I_{2 n_{i}} \end{matrix}] \in ℝ^{2 (| ℕ_{i}^{0} | + 1) n_{i} \times 2 n_{i}}, {j_{1}, j_{2}, \dots, j_{| ℕ_{i}^{0} |}} = ℕ_{i}^{0}

and $W_{i} (k) = diag {w_{i,1} (k) I_{2 n_{i}}, \dots, w_{i, | ℕ_{i}^{0} |} (k) I_{2 n_{i}}, α_{i} I_{2 n_{i}}}, where w_{i, s} (k) = {\begin{matrix} β_{i}, & the s th element of ℕ_{i}^{0} belongs to ℕ_{i} (k); \\ 0, & otherwise \end{matrix} for s \in {1, 2, \dots, | ℕ_{i}^{0} |}$

According to equation (11), the brief expression of the global cost function can be formulated for the multi-agent system as

J_{Σ} (k) = \sum_{i \in ℕ_{a}} J_{i} (x_{i r}, x_{i ▪}, u_{i r}, k) = \sum_{t \in ℕ} (∥ z_{Σ r} (t | k) ∥_{Q (k)}^{2} + ∥ u_{Σ r} (t | k) ∥_{R}^{2}) + ∥ z_{Σ r} (N | k) ∥_{P (k)}^{2}

where $z_{Σ r} = [x_{1 r}; x_{2 r}; \dots; x_{N_{a}, r}] \in ℝ^{2 n} and u_{Σ r} = [u_{1 r}; u_{2 r}; \dots; u_{N_{a}, r}] \in ℝ^{m}$ are the global collected vectors, with $n = \sum_{i \in ℕ_{a}} n_{i} and m = \sum_{i \in ℕ_{a}} m_{i}; Q (k) > 0, R > 0 and P (k) > 0$ are the associated global weighting matrices determined as in proposition 2.

Proposition 2

$P (k) = diag {P_{1} (k), \dots, P_{N_{a}} (k)}, R = diag {R_{1}, \dots, R_{N_{a}}} and Q (k) = G^{T} W (k) G, where W (k) = diag {W_{1} (k), \dots, W_{N_{a}} (k)} with W_{i} (k)$ formulated in proposition 1, and $G = [G_{1}; \dots; G_{N_{a}}]$ with

G_{i} = [\begin{matrix} G_{i}^{1, 1} & \dots & G_{i}^{1, N_{a}} \\ ⋮ & ⋱ & ⋮ \\ G_{i}^{| ℕ_{i}^{0} | + 1, 1} & ⋮ & G_{i}^{| ℕ_{i}^{0} | + 1, N_{a}} \end{matrix}] \in ℝ^{2 (| ℕ_{i}^{0} | + 1) n_{i} \times 2 n}, G_{i}^{p q} = {\begin{array}{l} I_{2 n_{i}}, & q = i, p = 1, \dots, | ℕ_{i}^{0} |, | ℕ_{i}^{0} | + 1; \\ - F_{i, q}^{x} (\cdot), & q is the p th element of ℕ_{i}^{0}; \\ 0, & otherwise \end{array}

Due to the time-varying matrix $W_{i} (k)$ has a natural upper bound ${\bar{W}}_{i} = diag {β_{i} I_{2 n_{i}}, \dots, β_{i} I_{2 n_{i}}, α_{i} I_{2 n_{i}}} for any i \in ℕ_{a}$ , the upper bounds of the real-time local weighting matrix $Q_{i} (k)$ and global weighting matrix Q(k) can be easily determined as ${\bar{Q}}_{i} = G_{i}^{T} {\bar{W}}_{i} G_{i} and \bar{Q} = G^{T} diag {{\bar{W}}_{1}, \dots, {\bar{W}}_{N_{a}}} G$ , respectively.

Stability condition for the multi-agent system with primary distributed RHC

According to the traditional RHC synthesis approach in the study by Mayne et al.,² imposing $x_{i} (N | k) \in X_{i}^{f} (N | k) with X_{i}^{f} (N | k)$ being a terminal constraint set is a common practice for stability.^5,18 Generally, the terminal constraint set $X_{i}^{f} (N | k)$ is chosen to be admissible and positively invariant. Here, we intend to describe the admissible and positively invariant set for each agent i, by the following terminal constraints

∥ x_{i r} (N | k) ∥_{P_{i} (k)}^{2} \leq ε_{i} (k)

∥ x_{i r}^{K} (t + 1 | k) ∥_{P_{i} (k)}^{2} - ∥ x_{i r}^{K} (t | k) ∥_{P_{i} (k)}^{2} \leq - ∥ x_{i r}^{K} (t | k) ∥_{\bar{λ} (\bar{Q}) I_{2 n_{i}}}^{2} - ∥ u_{i r}^{K} (t | k) ∥_{R_{i}}^{2}, \forall t \geq N

x_{i}^{K} (t | k) \in X_{i}, \forall t \geq N

u_{i}^{K} (t | k) = K_{i} (k) x_{i r}^{K} (t | k) + F_{i r}^{u} (u_{r} (k + t)) \in U_{i}, \forall t \geq N

where $ε_{i} (k)$ is a positive scalar; $x_{i r}^{K} and u_{i r}^{K}$ are the deviations corresponding to $x_{i}^{K} and u_{i}^{K}$ , respectively; $x_{i}^{K} (t + 1 | k)$ is obtained along the dynamics (1) by implementing $u_{i}^{K} (t | k) = K_{i} (k) x_{i r}^{K} (t | k) + F_{i r}^{u} (u_{r} (k + t)), with x_{i r}^{K} (N | k) = x_{i r}^{*} (N | k); K_{i} (k)$ is the terminal state feedback law. Hence, $X_{i}^{f} (N | k) = {x_{i} (N | k) | (15) - (18)}$ is defined. For $x_{i} (N | k) \in X_{i}^{f} (N | k)$ , by summing (16) for all $i \in ℕ_{a}$ , with respect to equation (18), we have

∥ z_{Σ r}^{K} (t + 1 | k) ∥_{P (k)}^{2} - ∥ z_{Σ r}^{K} (t | k) ∥_{P (k)}^{2} \leq - ∥ z_{Σ r}^{K} (t | k) ∥_{\bar{Q} + K^{T} (k) R K (k)}^{2}, \forall t \geq N

where $z_{Σ r}^{K} = [x_{1 r}^{K}; x_{2 r}^{K}; \dots; x_{N_{a}, r}^{K}] and K (k) = diag {K_{1} (k), \dots, K_{N_{a}} (k)}$ .

Remark 1

In this article, the parameters $ε_{i} (k), P_{i} (k) and K_{i} (k)$ are solved online. Alternatively, one can fix ε_i, P_i and K_i offline⁵ for less online computational burden. Compared with the offline choice, online solution with relaxed terminal constraints can improve the control performance and enlarge the attraction region,^18,24 which prefers to be adopted in this article.

To implement an RHC strategy in the distributed way, each agent i has to utilize the future evolution $x_{j} (t | k), j \in ℕ_{i} (k)$ to minimize the ideal cost function (10). Because the future evolution is not determined before optimization, an assumed evolution ${\hat{x}}_{j} (t | k)$ has to be adopted as substitution. Generally, the assumed evolution is chosen as the optimal one solved at the last instant,^5,6,13,18 that is, each agent $i \in ℕ_{a} broadcasts {\hat{x}}_{i} (t | k) = x_{i}^{*} (t + 1 | k - 1), t \in ℕ$ as its future evolution before the current optimization. This manner brings to the actual cost function $J_{i} (x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)$ , which is constructed by substituting $x_{i j} (t | k)$ in equation (10) with ${\hat{x}}_{i j} (t | k) = x_{i} (t | k) - F_{i j}^{x} ({\hat{x}}_{j} (t | k)) + d_{i j}^{x}$ .

Following the above procedures, each agent $i \in ℕ_{a}$ can apply a primary DRHC strategy by solving the individual receding horizon optimization problem

min_{ε_{i} (k), P_{i} (k), K_{i} (k), u_{i} (t | k), t \in ℕ} J_{i} (x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)

subject to x_{i} (0 | k) = x_{i} (k)

x_{i} (t + 1 | k) = A_{i} x_{i} (t | k) + B_{i} u_{i} (t | k), t \in ℕ

x_{i} (t + 1 | k) \in X_{i}, t \in ℕ^{+}

u_{i} (t | k) \in U_{i}, t \in ℕ

x_{i} (N | k) \in X_{i}^{f} (N | k) = {x_{i} (N | k) | (15) - (18)}

and implementing the first element of the optimal control input, that is, $u_{i} (k) = u_{i}^{*} (0 | k)$ . From the property of traditional RHC synthesis approach, the recursive feasibility of the distributed optimization problem (20 to 25) can be naturally guaranteed, by choosing the current feasible solution as

ε_{i}^{\circ} (k + 1) = ε_{i}^{*} (k), P_{i}^{\circ} (k + 1) = P_{i}^{*} (k), K_{i}^{\circ} (k + 1) = K_{i}^{*} (k)

u_{i}^{\circ} (t | k + 1) = u_{i}^{*} (t + 1 | k) for t \in ℕ^{-}, u_{i}^{\circ} (N - 1 | k + 1) = u_{i}^{K} (N | k)

However, due to $x_{j} (t | k)$ is substituted by ${\hat{x}}_{j} (t | k)$ , and the time-varying neighbourhood is considered in the actual cost function, the guarantee for the stability of overall multi-agent system requires an additional condition. The required condition can be derived as follows.

According to the common practice of RHC synthesis, we can choose the Lyapunov function as the optimal cost function of the overall multi-agent system, that is, $V (k) = \sum_{i \in ℕ_{a}} J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ . Then, by applying the feasible solution (26 and 27) and utilizing the equation (19), we have

\begin{array}{l} V (k + 1) - V (k) \leq \sum_{i \in ℕ_{a}} J_{i} (x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k + 1) - \sum_{i \in ℕ_{a}} J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) \\ = \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} (α_{i} ∥ x_{i r}^{\circ} (t | k + 1) ∥^{2} + ρ_{i} ∥ u_{i r}^{\circ} (t | k + 1) ∥^{2} + \sum_{j \in ℕ_{i} (k + 1)} β_{i} ∥ {\hat{x}}_{i j}^{\circ} (t | k + 1) ∥) + ∥ z_{Σ r}^{\circ} (N | k + 1) ∥_{P^{\circ} (k + 1)}^{2} \\ - \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} (α_{i} ∥ x_{i r}^{*} (t | k) ∥^{2} + ρ_{i} ∥ u_{i r}^{*} (t | k) ∥^{2} + \sum_{j \in ℕ_{i} (k)} β_{i} ∥ {\hat{x}}_{i j}^{*} (t | k) ∥) - ∥ z_{Σ r}^{*} (N | k) ∥_{P^{*} (k)}^{2} \\ \leq \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} (\sum_{j \in ℕ_{i} (k + 1)} β_{i} ∥ x_{i j}^{*} (t | k) ∥^{2} - \sum_{j \in ℕ_{i} (k)} β_{i} ∥ {\hat{x}}_{i j}^{*} (t | k) ∥^{2}) - \sum_{i \in ℕ_{a}} (∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2} + ∥ u_{i r}^{*} (0 | k) ∥_{R_{i}}^{2}) \end{array}

By introducing a pre-specified convergence speed $δ \in [0, 1)$ , we intend to guarantee the exponential stability of the overall multi-agent system, by imposing $V (k + 1) - V (k) \leq - (1 - δ) \underline{λ} (Q (k)) ∥ z_{Σ r} (k) ∥^{2} \leq 0 for all k \geq 0$ . Following equation (28), it can be strictly guaranteed by

\sum_{i \in ℕ_{a}} [\sum_{t \in ℕ} (\sum_{j \in ℕ_{i} (k + 1)} β_{i} ∥ x_{i j}^{*} (t | k) ∥^{2} - \sum_{j \in ℕ_{i} (k)} β_{i} ∥ {\hat{x}}_{i j}^{*} (t | k) ∥^{2}) - δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}] \leq 0

Because the neighbour index set $ℕ_{i} (\cdot)$ is time-varying according to the real-time communication topology, one has to consider the extreme case $ℕ_{i} (k) = Ø$ . Hence, the satisfaction of the inequality (29) further requires

\sum_{i \in ℕ_{a}} (\sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k + 1)} β_{i} ∥ x_{i j}^{*} (t | k) ∥^{2} - δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}) \leq 0, k \geq 0

which serves as a sufficient exponential stability condition for the overall multi-agent system.

Remark 2

For the multi-agent system with fixed communication topology, that is, $ℕ_{i} (k) = ℕ_{i} (k + 1) = ℕ_{i}^{0}$ for all k ≥ 0, the inequality (29) also serve as a sufficient exponential stability condition (similar to the condition (38) in the study by Wang and Ding⁵). Then by properly distributing, an individual constraint (see the compatibility constraints in the literature^5,6,18) can be designed and imposed in the individual optimization problem to guarantee stability. However, the time-varying communication topology leads to the sufficient exponential stability condition (30) in this article, which invalidates the design of compatibility constraint. Note that the condition (30) is inevitably more conservative than (29), because the extreme case $ℕ_{i} (k) = Ø$ must be taken into account.

DRHC synthesis

Based on the sufficient exponential stability condition (30), we synthesize the DRHC in this section as follows.

Distribution of the sufficient stability condition

By applying the definition $x_{i j} = x_{i r} - F_{i j}^{x} (x_{j r})$ and the triangle inequality, we have

\begin{array}{l} \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k + 1)} β_{i} ∥ x_{i j}^{*} (t | k) ∥^{2} \leq \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k + 1)} 2 β_{i} (∥ x_{i r}^{*} (t | k) ∥^{2} + ∥ F_{i j}^{x} (x_{j r}^{*} (t | k)) ∥^{2}) \\ \leq \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k + 1)} 2 β_{i} (∥ x_{i r}^{*} (t | k) ∥^{2} + {\bar{F}}_{i j}^{x 2} ∥ x_{j r}^{*} (t | k) ∥^{2}) \\ = \sum_{i \in ℕ_{a}} \sum_{t \in ℕ} 2 (\sum_{j \in ℕ_{i} (k + 1)} β_{i} + \sum_{i \in ℕ_{j} (k + 1)} β_{j} {\bar{F}}_{ji}^{x 2}) ∥ x_{i r}^{*} (t | k) ∥^{2} \end{array}

which can be substituted into equation (30) to obtain the following distributed condition

\sum_{t \in ℕ} 2 (\sum_{j \in ℕ_{i} (k + 1)} β_{i} + \sum_{i \in ℕ_{j} (k + 1)} β_{j} {\bar{F}}_{i j}^{x 2}) ∥ x_{i r}^{*} (t | k) ∥^{2} - δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2} \leq 0, \forall i \in ℕ_{a}

Note in condition (32) that, only the term $∥ x_{i r} (t | k) ∥^{2}$ is to be solved. Hence, we can assign the condition (32) to each agent as an individual constraint, that is

\sum_{t \in ℕ} ∥ x_{i r}^{*} (t | k) ∥^{2} \leq \frac{δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 (\sum_{j \in ℕ_{i} (k + 1)} β_{i} + \sum_{i \in ℕ_{j} (k + 1)} β_{j} {\bar{F}}_{i j}^{x 2})}

Design of an auxiliary constraint

Now, the exponential stability of the overall multi-agent system can be easily guaranteed, by imposing the individual constraint (33) in the optimization problem (20 to 25). However, the constraint (33) cannot be realized due to the undetermined sets $ℕ_{i} (k + 1)$ s, and its satisfaction can be prevented due to the fixed weighting scalars β_i s. Considering this issue, we enforce each agent i to update its formation weighting scalar according to

β_{i} (k + 1) = min_{j \in ℕ_{i} (k + 1)} {η_{i} (k), η_{j} (k), {\bar{β}}_{i}}

where $η_{i} (k)$ s are introduced positive decision variables, and ${\bar{β}}_{i}$ is specified as the originally fixed scalar β_i. According to the update equation (34), the following property holds

\sum_{j \in ℕ_{i} (k + 1)} β_{i} (k + 1) + \sum_{i \in ℕ_{j} (k + 1)} β_{j} (k + 1) {\bar{F}}_{j i}^{x 2} \leq \sum_{j \in ℕ_{i} (k + 1)} η_{i} (k) + \sum_{i \in ℕ_{j} (k + 1)} η_{i} (k) {\bar{F}}_{j i}^{x 2} \leq η_{i} (k) \sum_{j \in ℕ_{i}^{0}} (1 + {\bar{F}}_{j i}^{x 2})

Then, the individual constraint (33) is strictly guaranteed by an auxiliary constraint

\sum_{t \in ℕ} ∥ x_{i r} (t | k) ∥^{2} \leq \frac{δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 η_{i} (k) Δ_{i}}

where $Δ_{i} = \sum_{j \in ℕ_{i}^{0}} (1 + {\bar{F}}_{j i}^{x 2})$ . It deserves to note that, the constraint (36) can be easily realized by calculating the terms $∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}$ and Δ_i before optimization and can be easily satisfied by choosing a decision variable $η_{i} (k) > 0$ small enough.

Handling of the cost function

By imposing the auxiliary constraint (36) in the individual optimization problem (20 to 25), both the recursive feasibility and exponential stability can be guaranteed. However, because the auxiliary constraint (36) is more easily satisfiable for a smaller $η_{i} (k)$ , directly minimizing the cost function $J_{i} (x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)$ tends to obtain a smaller $η_{i} (k)$ . Then, following the update equation (34), a smaller $β_{i} (k + 1)$ will be obtained to weaken the formation effort. Hence, we prefer to solve a larger variable $η_{i} (k)$ in the optimization. Additionally, the terminal constraint set $X_{i}^{f} (N | k)$ prefers a smaller variable $ε_{i} (k)$ . To these ends, we apply the auxiliary constraint (36) and the terminal constraint (15) to determine an upper bound for $J_{i} (x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)$ as

{\bar{J}}_{i} (ε_{i}, η_{i}, x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k) = (1 + \frac{α_{i} δ}{2 η_{i} (k) Δ_{i}}) ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2} + \sum_{t \in ℕ} (ρ_{i} ∥ u_{i r} (t | k) ∥^{2} + \sum_{j \in ℕ_{i} (k)} β_{i} (k) ∥ {\hat{x}}_{i j} (t | k) ∥^{2}) + ε_{i} (k)

Note in equation (37) that, the upper bound ${\bar{J}}_{i} (ε_{i}, η_{i}, x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)$ is inversely proportional to the variable $η_{i} (k)$ and proportional to the variable $ε_{i} (k)$ . Hence, the maximum $η_{i} (k)$ and the minimum $ε_{i} (k)$ can be solved by minimizing equation (37).

Implementation of DRHC strategy

Following the above procedures, the actual individual optimization problem can be constructed for each agent $i \in ℕ_{a}$ as

min_{ε_{i} (k), P_{i} (k), K_{i} (k), η_{i} (k), u_{i} (t | k), t \in ℕ} (37) subject to (21) - (25) and (36)

based on which the associated DRHC algorithm is presented for implementation as follows.

Algorithm 1

Offline stage: For each agent $i \in ℕ_{a}$ , determine the set $ℕ_{i}^{0}$ according to the given tracking and formation vectors $d_{i r}^{p} and d_{i j}^{p}$ ; determine the functions $F_{i r}^{x} (\cdot), F_{i r}^{u} (\cdot) and F_{i j}^{x} (\cdot)$ s according to the dynamics of all agents and calculate the Lipschitz constants ${\bar{F}}_{i j}^{x}$ s; give the weighting scalars α_i, ρ_i and $β_{i} = {\bar{β}}_{i}$ and choose the convergence speed $δ \in [0, 1)$ and predictive horizon N ; calculate the scalars $\bar{λ} (\bar{Q}) and Δ_{i} = \sum_{j \in ℕ_{i}^{0}} (1 + {\bar{F}}_{i j}^{x 2})$ .

Online stage: At the initial instant k = 0, set $η_{i}^{*} (- 1) = {\bar{β}}_{i} and {\hat{x}}_{i} (t | 0) = x_{i} (0) for all i \in ℕ_{a} and t \in ℕ$ . At each instant k ≥ 0, each agent i performs as follows:

Broadcast $η_{i}^{*} (k - 1) and {\hat{x}}_{i} (t | k)$ , t ∈ ℕ and receive $η_{j}^{*} (k - 1) and {\hat{x}}_{j} (t | k)$ , t ∈ ℕ from its desired neighbours $j \in ℕ_{i}^{0}$ .

Determine the current set $ℕ_{i} (k)$ according to whether it received from $j \in ℕ_{i}^{0}$ , update the current scalar $β_{i} (k)$ according to equation (34) and generate the current cost according to equation (37).

Solve the optimization problem (38) to obtain the optimal solution ${ε_{i}^{*} (k), P_{i}^{*} (k), K_{i}^{*} (k), η_{i}^{*} (k), u_{i}^{*} (t | k), t \in ℕ}$ .

Implement the current control input $u_{i} (k) = u_{i}^{*} (0 | k)$ .

According to the procedure in ‘Handling of the cost function’ section, the optimization of equation (38) in step 1 will solve the minimum variable $η_{i} (k)$ . Then, the auxiliary constraint (36) can restrict the future tracking error $\sum_{t \in ℕ} ∥ x_{i r} (t | k) ∥^{2}$ to be small, that is, can enhance the tracking effort. Further according to the procedure in ‘Design of an auxiliary constraint’ section, updating according to equation (34) in step 2 will provide a larger formation weighting scalar $β_{i} (k + 1)$ . Then, the actually minimized cost (37) can punish the formation error term more heavily, that is, can enhance the formation effort. Hence, the auxiliary constraint (36), with the aid of update (34) and cost (37), plays an important role in the simultaneous tracking and formation enhancement. This is substantially different from the role of compatibility constraint imposed in the literature.^5,6,13,18

In DRHC strategy, the optimization time is generally assumed to be trivial for immediate implementation of control input; however, it is non-trivial in practice. For the real-time implementation with non-trivial optimization time (assumed to be smaller than t times sampling period), a predictive version of DRHC strategy (as referred in the literature^6,13) can be adopted. In a common predictive DRHC strategy, each agent solves in advance a predictive version of optimization problem at the instant k − t, with substituting the actual state $x_{i} (k)$ by the predictive state $x_{i}^{*} (t | k - t)$ . Because the real-time neighbour index set $ℕ_{i} (k)$ is considered for the time-varying communication topology in this article, one has to substitute it by $ℕ_{i} (k - t)$ in the predictive DRHC strategy. Note that, this substitution may lead to the sluggish formation of the multi-agents, which is the price for the real-time implementation in time-varying communication environment.

Property analysis

Since the above DRHC strategy is proposed according to the rationale of synthesis approach, the properties of a general synthesis approach (i.e. the recursive feasibility and stability) can be provided, which are analysed in the following.

Theorem 5.1

For each agent $i \in ℕ_{a}$ , if the individual optimization problem (38) is feasible at the initial instant k = 0, then by applying algorithm 1,

[Recursive Feasibility] the optimization problem (38) is feasible at any instant k > 0;

[Equivalent Optimality] the values of $J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) and {\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ are always identical, that is, for any k ≥ 0

∥ x_{i r}^{*} (N | k) ∥_{P_{i}^{*} (k)}^{2} = ε_{i}^{*} (k)

\sum_{t \in ℕ} ∥ x_{i r}^{*} (t | k) ∥^{2} = \frac{δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 η_{i}^{*} (k) Δ_{i}}

and the value of

J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)

is always minimal, that is

J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) \leq J_{i} (x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k), \forall k \geq 0

[Exponential Stability] the overall multi-agent system is exponentially stable, that is, the tracking and formation errors exponentially converge to zeros.

Proof

a) At the instant k, suppose that the optimization problem (38) is feasible, and the optimal solution is solved as ${ε_{i}^{*} (k), P_{i}^{*} (k), K_{i}^{*} (k), η_{i}^{*} (k), u_{i}^{*} (t | k), t \in ℕ}$ . By applying algorithm 1, two cases can be induced at the next instant k + 1.

Case 1

$∥ z_{i r} (k + 1) ∥_{Q_{i} (k + 1)}^{2} > 0$ , which implies that the simultaneous tracking and formation objective is not achieved. Then, by applying the control input sequence ${u_{i}^{\circ} (0 | k + 1), \dots, u_{i}^{\circ} (N - 1 | k + 1)} = {u_{i}^{*} (1 | k), \dots, u_{i}^{*} (N - 1 | k), u_{i}^{K} (N | k)}$ , the associated state sequence can be obtained as ${x_{i}^{\circ} (0 | k + 1), \dots, x_{i}^{\circ} (N | k + 1)} = {x_{i}^{*} (1 | k), \dots, x_{i}^{*} (N | k), x_{i}^{K} (N + 1 | k)}$ . Following the terminal constraints (17) and (18), we have $x_{i}^{\circ} (t | k + 1) \in X_{i} for all t \in ℕ^{+}, u_{i}^{\circ} (t | k + 1) \in U_{i}$ for all t ∈ ℕ and $x_{i}^{\circ} (N | k + 1) \in X_{i}^{f} (N | k)$ . Hence, the constraints (21) to (25) are all satisfied with $ε_{i}^{\circ} (k + 1) = ε_{i}^{*} (k), P_{i}^{\circ} (k + 1) = P_{i}^{*} (k) and K_{i}^{\circ} (k + 1) = K_{i}^{*} (k)$ . Additionally, by choosing a sufficient small $η_{i}^{\circ} (k + 1) > 0$ , the auxiliary constraint (36) can be easily satisfied. Thus, a feasible solution of the optimization problem (38) can be found at the instant k + 1.

Case 2

$∥ z_{i r} (k + 1) ∥_{Q_{i} (k + 1)}^{2} = 0$ , which implies that the simultaneous tracking and formation objective (with the current neighbours) is achieved, that is, $x_{i} (k + 1) = F_{i r}^{x} (x_{r} (k + 1)) + d_{i r}^{x} = F_{i j}^{x} (x_{j} (k + 1)) + d_{i j}^{x} for all j \in ℕ_{i} (k + 1)$ . Then, by applying the control input sequence $u_{i}^{\circ} (t | k + 1) = F_{i r}^{u} (u_{r} (k + t + 1))$ , t ∈ ℕ, the associated state sequence can be obtained as $x_{i}^{\circ} (t | k + 1) = F_{i r}^{x} (x_{r} (k + 1 + t)) + d_{i r}^{x}, t \in {1, \dots, N}$ , which leads to $x_{i r}^{\circ} (t | k + 1) = 0 for any t \in {1, \dots, N}$ . Following the trackability assumption (8), the constraints (21) to (24) are all satisfied. For any $ε_{i}^{\circ} (k + 1) > 0, P_{i}^{\circ} (k + 1) > 0 and K_{i}^{\circ} (k + 1)$ , the terminal constraint (25) is naturally satisfied. Additionally, the equality in auxiliary constraint (36) holds for any $η_{i}^{\circ} (k + 1) > 0$ . Thus, a feasible solution of the optimization problem (38) can be found at the instant k + 1.

In summary, the feasibility at the instant k guarantees the feasibility at the instant k + 1. Hence the conclusion (a) holds by induction.

b) Suppose that the equality (39) does not hold for the optimal solution ${ε_{i}^{*} (k), P_{i}^{*} (k), K_{i}^{*} (k), η_{i}^{*} (k), u_{i}^{*} (t | k), t \in ℕ}$ , then the inequality $∥ x_{i r}^{*} (N | k) ∥_{P_{i}^{*} (k)}^{2} < ε_{i}^{*} (k)$ holds. By choosing the variable $ε_{i}^{\circ} (k) = ∥ x_{i r}^{*} (N | k) ∥_{P_{i}^{*} (k)}^{2} < ε_{i}^{*} (k)$ , we have ${\bar{J}}_{i} (ε_{i}^{\circ}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) < {\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ , which contradicts with the minimality of ${\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ . Hence, the condition (39) holds. Similarly, the condition (40) holds. Following the equalities (39) and (40), we have $J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) = {\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ .

Suppose that the condition (41) does not hold for the optimal solution ${ε_{i}^{*} (k), P_{i}^{*} (k), K_{i}^{*} (k), η_{i}^{*} (k), u_{i}^{*} (t | k), t \in ℕ}$ , then there exist some feasible variables $P_{i}^{\circ} (k), K_{i}^{\circ} (k) and u_{i}^{\circ} (t | k)$ , t ∈ ℕ such that $J_{i} (x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k) < J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ . By choosing the variables $ε_{i}^{\circ} (k) = ∥ x_{i r}^{\circ} (N | k) ∥_{P_{i}^{\circ} (k)}^{2} and η_{i}^{\circ} (k) = \frac{δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 Δ_{i} Σ_{t \in ℕ} ∥ x_{i r}^{\circ} (t | k) ∥^{2}}$ , we have

{\bar{J}}_{i} (ε_{i}^{\circ}, η_{i}^{\circ}, x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k) = J_{i} (x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k) < J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) = {\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)

which contradicts with the minimality of ${\bar{J}}_{i} (ε_{i}^{*}, η_{i}^{*}, x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k)$ . Hence, the condition (41) holds.

c) Following the condition (41) and the procedures in ‘Stability condition for the multi-agent system with primary distributed RHC’, ‘Distribution of the sufficient stability condition’ and ‘Design of an auxiliary constraint’ sections, the recursive satisfaction of the auxiliary constraint (36) guarantees $V (k + 1) - V (k) \leq - (1 - δ) \underline{λ} (Q (k)) ∥ z_{Σ r} (k) ∥^{2} \leq 0$ for any k ≥ 0. Hence, the conclusion (c) holds.

These complete the proof. □

Now, following the proof of theorem 5.1, the main results of this article can be concluded as follows. For a nonhomogeneous multi-agent system (1) with time-varying communication topology, a reference agent (4) and a given control objective (5) and (6) satisfying assumption 1, if the optimization problem (38) of each agent $i \in ℕ_{a}$ is feasible at the initial instant k = 0, then the implementation of algorithm 1 for each agent can steer the multi-agents to achieve the simultaneous tracking and formation without violating the desired constraints (2) and (3).

It deserves to note that, both the recursive feasibility and stability are guaranteed independently of the real-time communication topology. Hence, the above results hold also for the real-time implementation of algorithm 1, in which the real-time set $ℕ_{i} (k)$ is substituted by $ℕ_{i} (k - t)$ . Indeed, with the implementation of the predictive version of algorithm 1, only the tracking and formation performance is affected by the time-varying communication topology.

The optimization problem in LMIs form

In this section, an LMIs form of the individual optimization problem (38) is provided via the following steps.

State constraint (23): The constraint (23) with $X_{i}$ in equation (2) can be expressed as $∥ Ψ_{i, s} x_{i} (t | k) ∥^{2} = x_{i}^{T} (t | k) Ψ_{i, s}^{T} Ψ_{i, s} x_{i} (t | k) \leq {\bar{ψ}}_{i, s}^{2} for all s = 1, \dots, c_{i}^{x} and t \in ℕ^{+}$ , which by applying Schur complement is equivalent to

[\begin{matrix} - {\bar{ψ}}_{i, s}^{2} & ⋆ \\ Ψ_{i, s} x_{i} (t | k) & - I_{n_{i, s}^{x}} \end{matrix}] \leq 0, s = 1, 2, \dots, c_{i}^{x}, t \in ℕ^{+}

Input constraint (24): Similarly to the above procedure, the constraint (24) with $U_{i}$ in equation (3) is equivalent to

[\begin{matrix} - {\bar{φ}}_{i, s}^{2} & ⋆ \\ Φ_{i, s} u_{i} (t | k) & - I_{m_{i, s}^{u}} \end{matrix}] \leq 0, s = 1, 2, \dots, c_{i}^{u}, t \in ℕ

Terminal constraint (25): For the satisfaction of the constraint (25), we can guarantee the terminal constraints (15) to (18) as follows (similar to the procedures in the literature^5,24).

Constraint (15): Define $X_{i} (k) = ε_{i} (k) P_{i}^{- 1} (k) and Y_{i} (k) = K_{i} (k) X_{i} (k)$ . The constraint (15) can be reformulated as $∥ x_{i r} (N | k) ∥_{X_{i}^{- 1} (k)}^{2} = x_{i r}^{T} (N | k) X_{i}^{- 1} (k) x_{i r} (N | k) \leq 1$ , which by applying Schur complement is equivalent to

[\begin{matrix} - 1 & ⋆ \\ x_{i r} (N | k) & - X_{i} (k) \end{matrix}] \leq 0

Constraint (16): By implementing the terminal control input $u_{i}^{K} (t | k) = K_{i} (k) x_{i r}^{K} (t | k) + F_{i r}^{u} (u_{r} (k + t))$ , the constraint (16) can be guaranteed by ${(A + B K_{i} (k))}^{T} P_{i} (k) (A + B K_{i} (k)) - P_{i} (k) + \bar{λ} (\bar{Q}) I_{2 n_{i}} + K_{i}^{T} (k) R_{i} K_{i} (k) \leq 0$ , which by applying Schur complement and congruent transformation is equivalent to

[\begin{matrix} - X_{i} (k) & ⋆ & ⋆ & ⋆ \\ A X_{i} (k) + B Y_{i} (k) & - X_{i} (k) & ⋆ & ⋆ \\ X_{i} (k) & 0 & - \frac{ε_{i} (k)}{\bar{λ} (\bar{Q})} I_{2 n_{i}} & ⋆ \\ Y_{i} (k) & 0 & 0 & - \frac{ε_{i} (k)}{ρ_{i}} I_{m_{i}} \end{matrix}] \leq 0

Note that, applying the LMIs (45) and (46) yields

∥ x_{i r}^{K} (t | k) ∥_{X_{i}^{- 1} (k)}^{2} \leq 1 for all t \geq N

Constraint (17): Using the triangle inequality, the constraint (17), that is, $∥ Ψ_{i, s} x_{i}^{K} (t | k) ∥ = ∥ Ψ_{i, s} [x_{i r}^{K} (t | k) + (F_{i r}^{x} (x_{r} (k + t)) + d_{i r}^{x})] ∥ \leq {\bar{ψ}}_{i, s}$ , is guaranteed by $∥ Ψ_{i, s} x_{i r}^{K} (t | k) ∥ \leq {\bar{ψ}}_{i, s} - {\bar{ψ}}_{r, s} with {\bar{ψ}}_{r, s} = {max}_{k \geq 0} ∥ Ψ_{i, s} (F_{i r}^{x} (x_{r} (k)) + d_{i r}^{x}) ∥$ . Following the property $∥ x_{i r}^{K} (t | k) ∥_{X_{i}^{- 1} (k)}^{2} \leq 1$ for all t ≥ N, we have $∥ Ψ_{i, s} x_{i r}^{K} (t | k) ∥ = ∥ Ψ_{i, s} X_{i}^{1 / 2} (k) X_{i}^{- 1 / 2} (k) x_{i r}^{K} (t | k) ∥ \leq ∥ Ψ_{i, s} X_{i}^{1 / 2} (k) ∥$ for all t ≥ N. Hence, the constraint (17) is guaranteed by $∥ Ψ_{i, s} X_{i}^{1 / 2} (k) ∥ \leq {\bar{ψ}}_{i, s} - {\bar{ψ}}_{r, s} for all s = 1, \dots, c_{i}^{x}$ , which by applying Schur complement is guaranteed by

[\begin{matrix} - X_{i} (k) & ⋆ \\ Ψ_{i, s} X_{i} (k) & - {({\bar{ψ}}_{i, s} - {\bar{ψ}}_{r, s})}^{2} I_{n_{i, s}^{x}} \end{matrix}] \leq 0, s = 1, 2, \dots, c_{i}^{x}

Constraint (18): Similar to the above procedure, the constraint (18) is guaranteed by

[\begin{matrix} - X_{i} (k) & ⋆ \\ Φ_{i, s} Y_{i} (k) & - {({\bar{φ}}_{i, s} - {\bar{φ}}_{r, s})}^{2} I_{m_{i, s}^{u}} \end{matrix}] \leq 0, s = 1, 2, \dots, c_{i}^{u}

where

{\bar{φ}}_{r, s} = {max}_{k \geq 0} ∥ Φ_{i, s} F_{i r}^{u} (u_{r} (k)) ∥

Auxiliary constraint (36): Define $θ_{i} (k) = η_{i}^{- 1} (k) and x_{i r}^{ℕ^{+}} = [x_{i r} (1 | k); \dots; x_{i r} (N - 1 | k)] \in ℝ^{2 (N - 1) n_{i}}$ . The constraint (36) by applying Schur complement is equivalent to

[\begin{matrix} - \frac{δ θ_{i} (k) ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 Δ_{i}} & ⋆ \\ x_{i r}^{ℕ^{+}} & - I_{2 (N - 1) n_{i}} \end{matrix}] \leq 0

Cost function (37): Define $u_{i r}^{ℕ} = [u_{i r} (0 | k); \dots; u_{i r} (N - 1 | k)] \in ℝ^{N m_{i}}, {\hat{x}}_{i j}^{ℕ^{+}} = [{\hat{x}}_{i j} (1 | k); \dots; {\hat{x}}_{i j} (N - 1 | k)] \in ℝ^{2 (N - 1) n_{i}} and {\hat{x}}_{i, ℕ_{i} (k)}^{ℕ^{+}} = [{\hat{x}}_{i j}^{ℕ^{+}}; \dots; {\hat{x}}_{i l}^{ℕ^{+}}] \in ℝ^{2 (N - 1) | ℕ_{i} (k) | n_{i}} with {j, \dots, l} = ℕ_{i} (k)$ . Let $γ_{i} (k) \geq \sum_{t \in ℕ} ρ_{i} ∥ u_{i r} (t | k) ∥^{2} + \sum_{t \in ℕ} \sum_{j \in ℕ_{i} (k)} β_{i} (k) ∥ {\hat{x}}_{i j} (t | k) ∥^{2} = ρ_{i} ∥ u_{i r}^{ℕ} ∥^{2} + β_{i} (k) ∥ {\hat{x}}_{i, ℕ_{i} (k)}^{ℕ^{+}} ∥^{2}$ , which by applying Schur complement is equivalent to

[\begin{matrix} - γ_{i} (k) & ⋆ & ⋆ \\ u_{i r}^{ℕ} & - ρ_{i}^{- 1} I_{N m_{i}} & ⋆ \\ {\hat{x}}_{i, ℕ_{i} (k)}^{ℕ^{+}} & 0 & - β_{i}^{- 1} (k) I_{2 (N - 1) | ℕ_{i} (k) | n_{i}} \end{matrix}] \leq 0

The cost

{\bar{J}}_{i} (ε_{i}, η_{i}, x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)

can be bounded from above by

{\bar{\bar{J}}}_{i} (θ_{i}, γ_{i}, ε_{i}, k) = (1 + \frac{α_{i} δ θ_{i} (k)}{2 Δ_{i}}) ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2} + γ_{i} (k) + ε_{i} (k)

Optimization problem (38): By incorporating the constraints (43) to (50) and the cost ${\bar{\bar{J}}}_{i} (θ_{i}, γ_{i}, ε_{i}, k)$ , the LMIs optimization problem is formulated as

min_{ε_{i} (k), X_{i} (k), Y_{i} (k), θ_{i} (k), γ_{i} (k), u_{i} (t | k), t \in ℕ} {\bar{\bar{J}}}_{i} (θ_{i}, γ_{i}, ε_{i}, k) subject to (43) - (50)

By applying algorithm 1 to each agent $i \in ℕ_{a}$ by solving the LMIs optimization problem (51), all the properties shown in ‘Property analysis’ section are still guaranteed, which can be similarly verified.

Illustrative examples

To demonstrate the effectiveness and practicability of our result, we provide a comparison and a verification in the following numerical and simulation examples, respectively.

Numerical example

Consider a very simple multi-agent system with $N_{a} = 2$ homogeneous agents ( $i = 1, 2$ ) in 1D space. The coordinate transform functions are $F_{12}^{x} (x) = F_{21}^{x} (x) = I_{2} x$ , and the Lipschitz constants are calculated as ${\bar{F}}_{12}^{x} = {\bar{F}}_{21}^{x} = 1$ . The communication topology is assumed to be time-varying due to the distance between agents, that is, each agent can receive data from the other only for $∥ d_{12} (k) ∥ = ∥ d_{21} (k) ∥ = ∥ p_{1} (k) - p_{2} (k) ∥ \leq \bar{d}, where \bar{d} = 5$ denotes the distance limitation of effective communication for each agent. Let us test the monotonous decrease of Lyapunov function, with the fixed formation weighting scalar (as in the literature^5,13) and the real-time updated formation weighting scalars (as in this article).

Give $ℕ_{1}^{0} = {2}, ℕ_{2}^{0} = {1}, α_{1} = α_{2} = ρ_{1} = ρ_{2} = 1 and N = 2$ and take the convergence speed as $δ = 0.9$ . For comparison, we choose the fixed formation weighting scalars as $β_{1} = β_{2} = 2$ and set the upper bounds ${\bar{β}}_{1} = {\bar{β}}_{2} = 2$ for the real-time updated formation weighting scalars $β_{1} (k) and β_{2} (k)$ .

Suppose the following conditions occur at the instant k

x_{1 r}^{*} (0 | k) = [4; 2], x_{1 r}^{*} (1 | k) = [2; 1], x_{1 r}^{*} (2 | k) = [0; 0], u_{1 r}^{*} (0 | k) = - 1, u_{1 r}^{*} (1 | k) = - 1

x_{2 r}^{*} (0 | k) = - [4; 2], x_{2 r}^{*} (1 | k) = - [2; 1], x_{2 r}^{*} (2 | k) = [0; 0], u_{2 r}^{*} (0 | k) = 1, u_{2 r}^{*} (1 | k) = 1

Then, the optimal costs of all agents are calculated as

\begin{matrix} J_{1} (x_{1 r}^{*}, {\hat{x}}_{12}^{*}, u_{1 r}^{*}, k) = 1 \times (∥ [4; 2] ∥^{2} + ∥ [2; 1] ∥^{2}) + 1 \times (∥ - 1 ∥^{2} + ∥ - 1 ∥^{2}) + ∥ [0; 0] ∥_{P_{1}^{*} (k)}^{2} = 27 \\ J_{2} (x_{2 r}^{*}, {\hat{x}}_{21}^{*}, u_{2 r}^{*}, k) = 1 \times (∥ - [4; 2] ∥^{2} + ∥ - [2; 1] ∥^{2}) + 1 \times (∥ 1 ∥^{2} + ∥ 1 ∥^{2}) + ∥ [0; 0] ∥_{P_{2}^{*} (k)}^{2} = 27 \end{matrix}

It deserves to note that, because $∥ d_{12} (k) ∥ = ∥ d_{21} (k) ∥ = ∥ p_{1 r}^{*} (0 | k) - p_{2 r}^{*} (0 | k) ∥ = ∥ 4 - (- 4) ∥ = 8 > \bar{d} induces ℕ_{1} (k) = ℕ_{2} (k) = Ø$ , the terms ${\hat{x}}_{12} (0 | k), {\hat{x}}_{12} (1 | k), {\hat{x}}_{21} (0 | k) and {\hat{x}}_{21} (1 | k)$ are not incorporated in the real-time costs.

At the instant k + 1, because $∥ d_{12} (k + 1) ∥ = ∥ d_{21} (k + 1) ∥ = ∥ p_{1 r}^{*} (1 | k) - p_{2 r}^{*} (1 | k) ∥ = ∥ 2 - (- 2) ∥ = 4 < \bar{d} induces ℕ_{1} (k + 1) = ℕ_{1}^{0} and ℕ_{2} (k + 1) = ℕ_{2}^{0}$ , the feasible state and input deviation sequences are obtained, according to case 1 in proof (a) of theorem 5.1 as

\begin{array}{l} x_{1 r}^{\circ} (0 | k + 1) = x_{1 r}^{*} (1 | k) = [2; 1], x_{1 r}^{\circ} (1 | k + 1) = x_{1 r}^{*} (2 | k) = [0; 0], x_{1 r}^{\circ} (2 | k + 1) = x_{1 r}^{K} (3 | k) = [0; 0] \\ {\hat{x}}_{12}^{\circ} (0 | k + 1) = x_{1}^{\circ} (0 | k + 1) - {\hat{x}}_{2} (0 | k + 1) = x_{1 r}^{*} (1 | k) - x_{2 r}^{*} (1 | k) = [4; 2], u_{1 r}^{\circ} (0 | k + 1) = u_{1 r}^{*} (1 | k) = - 1 \\ {\hat{x}}_{12}^{\circ} (1 | k + 1) = x_{1}^{\circ} (1 | k + 1) - {\hat{x}}_{2} (1 | k + 1) = x_{1 r}^{*} (2 | k) - x_{2 r}^{*} (2 | k) = [0; 0], u_{1 r}^{\circ} (1 | k + 1)) = K_{1}^{*} (k) x_{1 r}^{*} (2 | k) = 0 \end{array}

\begin{matrix} x_{2 r}^{\circ} (0 | k + 1) = x_{2 r}^{*} (1 | k) = - [2; 1], x_{2 r}^{\circ} (1 | k + 1) = x_{2 r}^{*} (2 | k) = [0; 0], x_{2 r}^{\circ} (2 | k + 1) = x_{2 r}^{K} (3 | k) = [0; 0] \\ {\hat{x}}_{21}^{\circ} (0 | k + 1) = x_{2}^{\circ} (0 | k + 1) - {\hat{x}}_{1} (0 | k + 1) = x_{2 r}^{*} (1 | k) - x_{1 r}^{*} (1 | k) = - [4; 2], u_{2 r}^{\circ} (0 | k + 1) = u_{2 r}^{*} (1 | k) = 1 \\ {\hat{x}}_{21}^{\circ} (1 | k + 1) = x_{2}^{\circ} (1 | k + 1) - {\hat{x}}_{1} (1 | k + 1) = x_{2 r}^{*} (2 | k) - x_{1 r}^{*} (2 | k) = [0; 0], u_{2 r}^{\circ} (1 | k + 1) = K_{2}^{*} (k) x_{2 r}^{*} (2 | k) = 0 \end{matrix}

With the fixed formation weighting scalars $β_{1} = β_{2} = 2$ , the associated feasible costs of all agents are calculated as

\begin{matrix} J_{1} (x_{1 r}^{\circ}, {\hat{x}}_{12}^{\circ}, u_{1 r}^{\circ}, k + 1) = 1 \times (∥ [2; 1] ∥^{2} + ∥ [0; 0] ∥^{2}) + 1 \times (∥ - 1 ∥^{2} + ∥ 0 ∥^{2}) + 2 \times (∥ [4; 2] ∥^{2} + ∥ [0; 0] ∥^{2}) + ∥ [0; 0] ∥_{P_{1}^{*} (k)}^{2} = 46 \\ J_{2} (x_{2 r}^{\circ}, {\hat{x}}_{21}^{\circ}, u_{2 r}^{\circ}, k + 1) = 1 \times (∥ - [2; 1] ∥^{2} + ∥ [0; 0] ∥^{2}) + 1 \times (∥ 1 ∥^{2} + ∥ 0 ∥^{2}) + 2 \times (∥ - [4; 2] ∥^{2} + ∥ [0; 0] ∥^{2}) + ∥ [0; 0] ∥_{P_{2}^{*} (k)}^{2} = 46 \end{matrix}

which yield $V (k + 1) - V (k) \leq \sum_{i = 1}^{2} J_{i} (x_{i r}^{\circ}, {\hat{x}}_{i ▪}^{\circ}, u_{i r}^{\circ}, k + 1) - \sum_{i = 1}^{2} J_{i} (x_{i r}^{*}, {\hat{x}}_{i ▪}^{*}, u_{i r}^{*}, k) \leq (46 + 46) - (27 + 27) = 38$ . Hence, the guarantee for the monotonous decrease of the Lyapunov function may be lost.

By applying algorithm 1, the conditions (52) and (53) yield (according to the condition (40))

η_{i}^{*} (k) = \frac{δ ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}}{2 Δ_{i} \sum_{t = 1}^{N - 1} ∥ x_{i r}^{*} (t | k) ∥^{2}} = \frac{0.9 \times 20}{2 \times 2 \times 5} = 0.9, i = 1, 2,

and $β_{1} (k + 1) = β_{2} (k + 1) = 0.9$ (according to the update equation (34)). Then, the feasible costs of all agents at the instant k + 1 are calculated as

\begin{matrix} J_{1} (x_{1 r}^{\circ}, {\hat{x}}_{12}^{\circ}, u_{1 r}^{\circ}, k + 1) = 1 \times (∥ [2; 1] ∥^{2} + ∥ [0; 0] ∥^{2}) + 1 \times (∥ - 1 ∥^{2} + ∥ 0 ∥^{2}) + 0.9 \times (∥ [4; 2] ∥^{2} + ∥ [0; 0] ∥^{2}) + ∥ [0; 0] ∥_{P_{1}^{*} (k)}^{2} = 24 \\ J_{2} (x_{2 r}^{\circ}, {\hat{x}}_{21}^{\circ}, u_{2 r}^{\circ}, k + 1) = 1 \times (∥ - [2; 1] ∥^{2} + ∥ [0; 0] ∥^{2}) + 1 \times (∥ 1 ∥^{2} + ∥ 0 ∥^{2}) + 0.9 \times (∥ - [4; 2] ∥^{2} + ∥ [0; 0] ∥^{2}) + ∥ [0; 0] ∥_{P_{2}^{*} (k)}^{2} = 24 \end{matrix}

These lead to $V (k + 1) - V (k) \leq (24 + 24) - (27 + 27) = - 6$ , that is, the monotonous decrease of the Lyapunov function is strictly guaranteed.

Simulation example

Consider a simultaneous tracking and formation problem of nonhomogeneous multi-agents with $N_{a} = 6$ vehicles in 3D space. The vehicles are composed of three 2D vehicles and three 3D vehicles, whose dynamics are formulated as in equation (1) with the parameters $A_{i} = [\begin{matrix} I_{2} & I_{2} \\ 0 & I_{2} \end{matrix}] and B_{i} = [\begin{matrix} 0.5 I_{2} \\ I_{2} \end{matrix}] for i = 1, 2, 3 and A_{i} = [\begin{matrix} I_{3} & I_{3} \\ 0 & I_{3} \end{matrix}] and B_{i} = [\begin{matrix} 0.5 I_{3} \\ I_{3} \end{matrix}] for i = 4, 5, 6$ . These dynamics models are obtained by discretizing the associated continuous-time models with sampling period 1 s. The initial states of these vehicles are given as $x_{1} (0) = [0; - 7.5; 0; 0], x_{2} (0) = [0; 4.5; 0; 0], x_{3} (0) = [0; - 1.5; 0; 0], x_{4} (0) = [0; - 4.5; 0; 0; 0; 0], x_{5} (0) = [0; 7.5; 0; 0; 0; 0] and x_{6} = [0; 1.5; 0; 0; 0; 0]$ .

The velocity and acceleration (i.e. input) constraints are desired for each vehicle $i \in ℕ_{a} as ∥ v_{i} (k) ∥ \leq {\bar{v}}_{i} and ∥ u_{i} (k) ∥ \leq {\bar{u}}_{i}$ , with the upper bounds ${\bar{v}}_{i} = \sqrt{2} and {\bar{u}}_{i} = 0.2$ . Additionally, the vertical acceleration constraint of 3D vehicles are desired as $∥ u_{i,3} (k) ∥ \leq {\bar{u}}_{v} for i = 4, 5, 6$ , with the upper bounds ${\bar{u}}_{v} = 0.05$ . By reformulating the above constraints according to equations (2) and (3), we can obtain the constraint parameters as $Ψ_{i} = {\begin{array}{l} [0, I_{2}] \in ℝ^{2 \times 4}, & i = 1, 2, 3 \\ [0, I_{3}] \in ℝ^{3 \times 6}, & i = 4, 5, 6 \end{array}, Φ_{i,1} = {\begin{array}{l} I_{2}, & i = 1, 2, 3 \\ I_{3}, & i = 4, 5, 6 \end{array}, {\bar{ψ}}_{i} = \sqrt{2} and {\bar{φ}}_{i,1} = 0.2 for i \in ℕ_{a}, Φ_{i,2} = [0, 0, 1] and {\bar{φ}}_{i,2} = 0.05 for i = 4, 5, 6$ .

The reference agent is a virtual vehicle with the same dynamic equation with the 2D vehicles, whose initial state and control input are given as $x_{r} (0) = [0; 0; 0; 0] and u_{r} (k) = [u_{r,1} (k); u_{r,2} (k)] = [{\bar{u}}_{r} cos (0.1 k); {\bar{u}}_{r} cos (0.25 k)] with {\bar{u}}_{r} = 0.1$ . The desired tracking and formation vectors are given as $d_{1 r}^{p} = [2; 0], d_{2 r}^{p} = [- 2; - 2], d_{3 r}^{p} = [- 2; 2], d_{4 r}^{p} = [- 2; 0; 4], d_{5 r}^{p} = [4; - 2; 4], d_{6 r}^{p} = [3; 2; 4], d_{12}^{p} = [4; 2], d_{13}^{p} = [4; - 2], d_{15}^{p} = [0; 2], d_{16}^{p} = [0; - 2], d_{21}^{p} = [- 4; - 2], d_{23}^{p} = [0; - 4], d_{24}^{p} = [0; - 2], d_{31}^{p} = [- 4; 2], d_{32}^{p} = [0; 4], d_{34}^{p} = [0; 2], d_{42}^{p} = [0; 2; 4], d_{43}^{p} = [0; - 2; 4], d_{45}^{p} = [- 4; 2; 0], d_{46}^{p} = [- 4; - 2; 0], d_{51}^{p} = [0; - 2; 4], d_{54}^{p} = [4; - 2; 0], d_{56}^{p} = [0; - 4; 0], d_{61}^{p} = [0; 2; 4], d_{64}^{p} = [4; 2; 0] and d_{65}^{p} = [0; 4; 0]$ .

In this simulation, the scalars $α_{i} = 1, ρ_{i} = 1 and {\bar{β}}_{i} = 1$ are given for each vehicle $i \in ℕ_{a}$ ; the scalars δ = 0.99 and N = 5 are chosen for all vehicles. According to the dynamics of these vehicles, the coordinate transformation functions can be determined as $F_{i r}^{x} (x) = {\begin{array}{l} I_{4} x, & i = 1, 2, 3 \\ diag {[I_{2}; 0], [I_{2}; 0]} x, & i = 4, 5, 6 \end{array}, F_{i r}^{u} (u) = {\begin{array}{l} I_{2} u, & i = 1, 2, 3 \\ [I_{2}; 0] u, & i = 4, 5, 6 \end{array},and F_{i j}^{x} (x) = {\begin{array}{l} diag {[I_{2},0], [I_{2},0]} x, & i = 1, 2, 3 and j = 4, 5, 6 \\ diag {[I_{2}; 0], [I_{2}; 0]} x, & i = 4, 5, 6 and j = 1, 2, 3 \\ I_{4} x, & i, j = 1, 2, 3; \\ I_{6} x, & i, j = 4, 5, 6, \end{array}$ ,

following which the Lipschitz constants can be calculated as ${\bar{F}}_{i j}^{x} = 1 for all (i, j) \in F$ . The parameters $\bar{λ} (\bar{Q}) \approx 13, Δ_{i} = 8 for i = 1, 4 and Δ_{i} = 6 for i = 2, 3, 5, 6$ can be calculated. The communication topology is time-varying due to the distance limitation of effective communication, which is set as $\bar{d} = 5$ for all vehicles. Note that, the communications between some desired pairs of vehicles fail at the initial instant, since the distances between these desired pairs are greater than $\bar{d}$ . Table 2 shows the real-time neighbourhood of vehicles over the simulation horizon.

Table 2.

Real-time neighbour index sets of all vehicles.

Instant k	0	1	2	3	4	5	6	7	8	9	10	11	12–30
$ℕ_{1} (k)$	Ø	Ø	Ø	Ø	Ø	Ø	Ø	{2}	{2}	{2,5}	{2,5}	{2,5,6}	{2,3,5,6}
$ℕ_{2} (k)$	Ø	Ø	Ø	{3}	{3}	{3,4}	{3,4}	{1,3,4}	{1,3,4}	{1,3,4}	{1,3,4}	{1,3,4}	{1,3,4}
$ℕ_{3} (k)$	{4}	{4}	{4}	{2,4}	{2,4}	{2,4}	{2,4}	{2,4}	{2,4}	{2,4}	{2,4}	{2,4}	{1,2,4}
$ℕ_{4} (k)$	{3}	{3}	{3}	{3}	{3}	{2,3}	{2,3}	{2,3,5,6}	{2,3,5,6}	{2,3,5,6}	{2,3,5,6}	{2,3,5,6}	{2,3,5,6}
$ℕ_{5} (k)$	Ø	Ø	Ø	{6}	{6}	{6}	{6}	{4,6}	{4,6}	{1,4,6}	{1,4,6}	{1,4,6}	{1,4,6}
$ℕ_{6} (k)$	Ø	Ø	Ø	{5}	{5}	{5}	{5}	{4,5}	{4,5}	{4,5}	{4,5}	{1,4,5}	{1,4,5}

By applying algorithm 1 by solving the optimization problem (51), we simulate this example through MATLAB on a PC with a dual-core 3.2 GHz Intel i5 CPU and a 3.47 GB RAM. Over the simulation horizon, all the optimal solutions are solved and implemented at each instant for each vehicle, that is, the recursive feasibility of the optimization problem (51) is demonstrated. Because the optimization problem (51) is an equivalent LMIs version of the optimization problem (38), the conclusion (a) of theorem 5.1 is verified. According to the solved solutions, the left and right terms of equation (39) are calculated and shown in Figure 1. It can be easily seen that the values of both sides in equation (39) are always identical for each vehicle; hence, the equality (39) is verified. Similarly, Figure 2 verifies the equality (40). These further verify the conclusion (b) of theorem 5.1 following the proof by contradiction. With the implementation of the solved solutions, the values of each vehicle’s real-time cost and the global cost are calculated and depicted in Figure 4. The monotonous decrease of the global cost verifies the conclusion (c) of theorem 5.1. More detailed simulation results are provided with analysis in the following.

Figure 1.

The left and right terms of equation (39), where ${Lt}_{i} (k) = ∥ x_{i r}^{*} (N | k) ∥_{P_{i}^{*} (k)}^{2} and {Rt}_{i} (k) = ε_{i}^{*} (k)$ .

Figure 2.

The left and right terms of equation (40), where ${Lt}_{i} (k) = \sum_{t \in ℕ} ∥ x_{i r}^{*} (t | k) ∥^{2} and {Rt}_{i} (k) = \frac{δ}{2 η_{i}^{*} (k) Δ_{i}} ∥ z_{i r} (k) ∥_{Q_{i} (k)}^{2}$ .

In Figure 3, the simulation trajectories of all vehicles and reference object are depicted, where the initial and final positions of all vehicles and reference object are marked with indices, and the desired formation is denoted by black dashed lines. It can be easily seen that all the 2D and 3D vehicles successfully track the reference object and achieve the desired formation, although the communication topology is time-varying. Additionally, the velocities and accelerations of all vehicles converge to those of the reference object (see Figures 5 to 7), which means that the vehicles will keep the desired tracking and formation in future. In Figures 5 to 7, the two-norms of all vehicles’ velocities do not exceed ${\bar{v}}_{i} = \sqrt{2}$ , the two-norms of all vehicles’ inputs (i.e. accelerations) do not exceed ${\bar{u}}_{i} = 0.2$ and the vertical accelerations of 3D vehicles $i = 4, 5, 6$ do not exceed the range $[- 0.05, 0.05]$ . These mean that all the desired state and input constraints are not violated over the simulation horizon. Thus, the effectiveness and practicability of the proposed DRHC strategy are illustrated.

Figure 3.

The trajectories of vehicles and reference object.

Figure 4.

The values of the cost functions.

Figure 5.

The two-norms of velocities.

Figure 6.

The two-norms of accelerations.

Figure 7.

The vertical accelerations of 3D vehicles.

To further verify the purpose discussed in ‘Implementation of DRHC strategy’ section (see the first paragraph after algorithm 1), we depict the update trends of the formation weighting scalars $β_{i} (k)$ s in Figure 8. It can be easily found that minimizing ${\bar{J}}_{i} (ε_{i}, η_{i}, x_{i r}, {\hat{x}}_{i ▪}, u_{i r}, k)$ tends to solve large $η_{i} (k)$ s and then to update large formation weighting scalars $β_{i} (k)$ s. For the real-time implementation, we show the optimization times $T c_{i} (k)$ of each vehicle $i \in ℕ_{a}$ in Figure 9. Note that, all the optimization times are smaller than the sampling period 1 s; hence, a real-time implementation can be simulated by solving in advance a predictive version of the optimization problem (51) at each instant k − 1 (more details refer to the second paragraph after algorithm 1). Because the simulation results are similar to those provided above, the detailed graphs of the simulation trajectories, costs, velocities and accelerations are omitted here.

Figure 8.

The real-time formation weighting scalars.

Figure 9.

The optimization times of all vehicles.

Conclusion

For the simultaneous tracking and formation problem of nonhomogeneous multi-agents with time-varying communication topology, we study the synthesis of DRHC in this article. By considering the control objective and real-time communication topology, we generate the individual cost function for each agent online, design an auxiliary constraint in distributed form and impose it in the optimization problem and propose an update equation for the formation weighting parameter. Accordingly, the individual optimization problem and associated DRHC algorithm are presented. By applying the presented DRHC algorithm, the recursive feasibility of the optimization problem and the exponential stability (i.e. the exponential convergence of tracking and formation) of the overall multi-agent system are guaranteed. Due to the handling of the effects of the time-varying communication topology on the recursive feasibility and stability, our result is preferred for miscellaneous practical applications.

Indeed, the practical application of nonhomogeneous multi-agents involves more factors, such as the collision avoidance, obstacle avoidance, formation switching, disturbance handling and so on, which will arouse a series of theoretical issues to be considered in our future researches.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Nature Science Foundation of China under Grants 61403414, 61571458 and 41601436, the Natural Science Basic Research Plan in Shaanxi Province of China under Grants 2016JQ6070 and 2016JM6050, the China Postdoctoral Science Foundation under Grant 2016M603042, and the Aero Science Fund under Grant 20160196005.

References

Garcia

Prett

Morari

. Model predictive control: theory and practice – a survey. Automatica 1989; 25(3): 335–348.

Mayne

Rawlings

Rao

. Constrained model predictive control: stability and optimality. Automatica 2000; 36(6): 789–814.

Camponogara

Jia

Krogh

. Distributed model predictive control. IEEE Control Syst Magaz 2002; 22(1): 44–52.

Scattolini

. Architectures for distributed and hierarchical model predictive control - a review. J Proc Control 2009; 19(5): 723–731.

Wang

Ding

. A synthesis approach of distributed model predictive control for homogeneous multi-agent system with collision avoidance. Int J Control 2014; 87(1): 52–63.

Wang

Ding

. Distributed RHC for tracking and formation of nonholonomic multi-vehicle systems. IEEE Trans Autom Control 2014; 59(6): 1439–1453.

Tanner

Jadbabaie

Pappas

. Flocking in fixed and switching networks. IEEE Trans Autom Control 2007; 52(5): 863–868.

Wang

Chen

. A connectivity-preserving flocking algorithm for multi-agent systems based only on position measurements. Int J Control 2009; 82(7): 1334–1343.

Yang

Zhang

Jiang

. Flocking of multi-agents with time delay. Int J Syst Sci 2012; 43(11): 2125–2134.

10.

Ferrari-Trecate

Galbusera

Marciandi

MPE

. Model predictive control schemes for consensus in multi-agent systems with single- and double-integrator dynamics. IEEE Trans Autom Control 2009; 54(11): 2560–2572.

11.

Shang

. Finite-time consensus for multi-agent systems with fixed topologies. Int J Syst Sci 2012; 43(3): 499–506.

12.

Cheng

Zhang

Fan

. Distributed consensus of multi-agent systems with input constraints: a model predictive control approach. IEEE Trans Circ Syst I 2015; 62(3): 825–834.

13.

Dunbar

Murray

. Distributed receding horizon control for multi-vehicle formation stabilization. Automatica 2006; 42(4): 549–558.

14.

Jiang

Zhou

. Consensus of multi-agent linear dynamic systems via impulsive control protocols. Int J Syst Sci 2011; 42(6): 967–976.

15.

Zhao

Zou

. A finite-time approach to formation control of multiple mobile robots with terminal sliding mode. Int J Syst Sci 2012; 43(11): 1998–2014.

16.

Lee

Kim

Myung

. Cooperative coevolutionary algorithm-based model predictive control guaranteeing stability of multirobot formation. IEEE Trans Control Syst Technol 2015; 23(1): 37–51.

17.

Keviczky

Borrelli

Balas

. Decentralized receding horizon control for large scale dynamically decoupled systems. Automatica 2006; 42(12): 2105–2115.

18.

Ding

Xie

Cai

. Distributed model predictive control for constrained linear systems. Int J Robust Nonlin Control 2010; 20(11): 1285–1298.

19.

Wang

Ong

. Distributed model predictive control of dynamically decoupled systems with coupled cost. Automatica 2010; 46(12): 2053–2058.

20.

Primbs

Nevistic

. Feasibility and stability of constrained finite receding horizon control. Automatica 2000; 36(7): 965–971.

21.

Primbs

Nevistic

. A new approach to stability analysis for constrained finite receding horizon control without end constraints. IEEE Trans Autom Control 2000; 45(8): 1507–1512.

22.

Olfati-Saber

Murray

. Consensus problems in networks of agents with switching topology and time-delays. IEEE Trans Autom Control 2004; 49(9): 1520–1533.

23.

Ren

Beard

. Consensus seeking in multiagent systems under dynamically changing interaction topologies. IEEE Trans Autom Control 2005; 50(5): 655–661.

24.

Kothare

Balakrishnan

Morari

. Robust constrained model predictive control using linear matrix inequalities. Automatica 1996; 32(10): 1361–1379.