Formation Control of Multi-agent Nonlinear Systems using The State-Dependent Riccati Equation

This paper presents research on controlling a highly-populated multi-agent system within the framework of the state-dependent Riccati equation (SDRE), which is a nonlinear optimal closed-loop control policy. A multi-agent system consists of a leader and a population of followers, playing the role of agents. The agents in a swarm system show a collective behavior, with interactions achieving a common goal [1]. The previous reports on swarm robotics have been quite diverse in subject matter: technological aspects of swarm robotics [2,3]; arrangement, formation, and behavior [4–6]; and bio-inspiration [7,8], among others.

The application of swarm robotics has been reported in many works, including micro-robotics: the interaction and physical contact between the agents [9–11], nano-swarms for delivering medicine or performing an operation in human body [12], search and rescue [13–15], monitoring and data collection [16, 17], etc.

Multiple decision-making agents with auto–nomous operation, social ability interaction, reactivity perception, and pro-activeness exhibition, were motivated by complex inherently distributed systems to increase robustness. Multi-agent and swarm systems are similar in the sense that they require multiple agents that communicate and cooperate. Each agent in multi-agent systems can do some meaningful part of a task; however, in swarm systems, we would expect each agent to be too unaware of its environment to possibly even know what is going on around it. This emphasizes the existence of a large number of agents, and promotes scalability, soft collisions and physical contacts of the swarm in migration. This might not be the case for swarms in nature, i.e., in birds, which do not collide, and also do not receive a command from leader [18]. Here, in this current work, the modeling is not based on bio-inspiration, since we are far away from such sophisticated design in mother nature.

Multi-agent systems, or cooperative robots, can usually be fit into these two categories: systems that considered complex dynamics for each agent, though the number of agents is small [19–24], and systems with simple dynamics or only kinematics for each agent when the number of agents is large [25–32]. For example, Yang et al. presented a formation control strategy based on a stable control law to guide the agents towards a constraint region and simulated it for a total of 155 agents [32]. Collision avoidance and obstacle avoidance between the agents are also two important characteristics of the formation methods. Wang and Xin presented an optimal control approach for formation control taking into consideration obstacle avoidance for multiple unmanned aerial vehicle (UAV) systems [33]. Kumar et al. presented a velocity controller for a swarm of UAVs with obstacle avoidance capability [34]. Park and Yoo proposed a connectivity-based formation control with collisionavoidance properties for a swarm of unmanned surface vessels [35,36].

The application of the swarm regulation control is devoted to multi-copter UAV formation control for a large number of agents–1,089–to show the capabilities of the proposed approach. Swarm control using augmented methods such as graph theory and multilayer network design is a useful approach; however, in this work, the pure capacity of the SDRE is used to control a highly populated multi-agent system. This will simplify the design and limit the formulation to pure SDRE. The advantages of this point of view result in controlling large-scale systems with complex dynamics. A fully coupled six-degree-of-freedom (DoF) nonlinear model of the system is considered (12 state variables). The second case study is a leaderfollower system of 45-wheeled mobile robots (WMRs) with non-holonomic and holonomic constraints of wheels. A smaller number of robots was considered to highlight the effectiveness of obstacle and collision avoidance. SDRE has been used in both multi-agent systems [37–39], and leader-follower control [22,40 – 47]. The details of the number of agents using SDRE are reported on Table 1. The majority of these studies focused on the control and behavior of the leaderfollower system, first considering two agents and then sometimes increasing the number of agents to 5 or 10. The only large populated system using SDRE control included 1,024 satellites; obstacle avoidance and collision avoidance between agents were not considered [39]. The large distance between the agents and the pattern of the satellites did not necessitate collision/obstacle avoidance. In this work, a more mature version of multi-agent SDRE [39] is presented to add those characteristics, since the agents are working close together in both solved examples.

Table 1.

A detailed report on the population of agents in SDRE in previous, as literature compared to this work.

context	Ref.	No. of agents, type of agent
multi‐agentleader‐followerformation	[37][38][40][41][42][43][44][45][46][47][39]	5, single‐integrator2, mobile robot2, spacecraft5, single‐integrator2, spacecraft3, mobile robot2, aircraft3, spacecraft2, spacecraft2, satellite1024, satellite
consensus control	[22]	10, crane
cooperativemulti‐agents	[48][49][50][19][51][52]	2, manipulator2, manipulator3, spacecraft4, manipulator4, UAV3, mobile robot
multi‐agentsleader‐follower	this workthis workthis work	1,089, UAV45, mobile robot1,050, mobile robot

The swarm formulation in control methods could be solved by graph theory for the formation design of the whole swarm; in that case, another controller must control the agents, individually by checking the connection for the agents through multiple layers [53–55]. The design of the formation using graph theory, and adapting a controller for handling the multiagent system, might pose difficulty in implementation because a controller must include other layers of designs through external methods.

The reported number of agents was 16 satellites [53], 6 agents [54], and 90 agents [55]. The use of the Kronecker delta function for the selection of information from a particular agent limits the application of the swarm multi-agent formulation to a limited number. In the SDRE method particularly, the solution to the Riccati equation might become overly complex if the SDRE should be solved for the entire swarm. In this approach, the gain of the targeted agent would be collected from the overall matrix of the swarm. Again, this approach would limit the implementation of the formulation for a large number of agents. The reported SDRE agents, using consensus control and graph theory, were 10 [22] and 2 [38], respectively.

The focus of this work is large-scale multi-agent systems that need simplicity in their design. This simplicity does not imply a simple linear controller; on the contrary, the SDRE is a nonlinear sub-optimal controller with a relatively complex solution. However, the practical approach to handle the multi-agent system relies solely on the SDRE’s capabilities itself, and uses the weighting matrix of the states to handle distance constraints and obstacle avoidance. As a result, there would not be further complexity for practical implementation raised by the interaction of the SDRE with another layer or graph theory. This can be viewed as independent SDRE controllers working with neighbor components and their leader through a weighting matrix, which is inside the umbrella of the SDRE, not another system augmented by the multiagents or other complementary methods. This is the key to large-scale modeling and simulation–which is simplicity in the communication of the multi-agent system and reliance on the capacities of the SDRE itself.

The main contribution of this work is the presentation of a control structure for forming a multi-agent system with complex dynamics with obstacle and collision avoidance between the agents. Complex systems with a large number of agents are seldom reported; this work presents a suitable structure for that. This work adds collision and obstacle avoidance to SDRE multi-agent system control, motivated by the example set by [39], which reported solely multi-agent system SDRE without contact prevention. For very crowded systems, contact between the agents has been both expected and considered in the literature [30,32], but this work proposes a method that adds distance between the agents to reduce collisions during the control task.

The distance-based arrangement of agents, and avoiding collision between them, is basically the idea of multi-agent system control methods that use graph theory, potential field method, multi-layer designs, and the like; here in this work, the distribution, as well as obstacle, and collision avoidance, are structured within the SDRE formulation–specifically, the weighting matrix of states.

Section 2 presents the control structure of the SDRE and multi-agent system modeling. Section 3 describes the dynamics of the two case studies of this work, UAV, and WMR. The simulation results are presented in Section 4 and the concluding remarks are summarized in Section 5.

2.

Control Structure: SDRE for Multi-agent Systems

2.1.

Conventional Controller: A Brief Preliminary Formulation

Consider N number of agents, randomly distributed in a predefined area, with continuous-time nonlinear systems and affine-in-control properties. The j-th agent’s dynamic is represented in state-space form as: 1 ${\dot{x}}_{j} (t) = f_{j} (x_{j} (t)) + g_{j} (x_{j} (t), u_{j} (t)),$ where x_j (t) ∈ ℝⁿ is the state vector and u_j (t) ∈ ℝ^m is the input vector of the j-th agent. f_j(x_j (t)) : ℝⁿ → ℝⁿ and g_j (x_j (t), u_j (t)) : ℝⁿ × ℝ^m → ℝⁿ in (1) are vectorvalued smooth functions that satisfy the local Lipschitz condition, with equilibrium point f_j(0) = 0. The nonlinear system (1) is transformed into the statedependent coefficient (SDC) parameterization: 2 ${\dot{x}}_{j} (t) = A_{j} (x_{j} (t)) x_{j} (t) + B_{j} (x_{j} (t)) u_{j} (t),$ where A_j(x_j(t)) : ℝⁿ → ℝ^n×n and B_j(x_j(t)) : ℝⁿ → ℝ^n×m are held.

The SDC parameterization (2) is referred to as apparent linearization; however, it preserves the nonlinearity of the system. For example, consider an over-actuated single-degree-of-freedom system in state-space form with n = 2 and m = 2, then the system vectors $\begin{array}{l} f (x) = [\begin{matrix} x_{2} \\ x_{1} x_{2}^{2} \cos x_{2} + x_{1} e^{x_{1}} \end{matrix}], \\ g (x, u) = [\begin{matrix} 0 \\ x_{2} u_{1} + x_{1} x_{2} u_{2} \end{matrix}], \end{array}$ can present a set of SDC matrices as $\begin{array}{l} f (x) & = \underset{A (x)}{\underset{︸}{[\begin{matrix} 0 & 1 \\ e^{x_{1}} & x_{1} x_{2} \cos x_{2} \end{matrix}]}} [\begin{array}{l} x_{1} \\ x_{2} \end{array}], \\ g (x, u) & = \underset{B (x)}{\underset{︸}{[\begin{matrix} 0 & 0 \\ x_{2} & x_{1} x_{2} \end{matrix}]}} [\begin{array}{l} u_{1} \\ u_{2} \end{array}] . \end{array}$

The result of multiplication of A(x)x and B(x)u will be the same original values of f(x) and g(x, u); hence, the SDRE method is a nonlinear optimal control that keeps the nonlinearities in the control structure. Another point on SDC matrices is that they are not unique representations, and can show numerous versions, such as: $\begin{array}{l} f (x) & = \underset{A (x)}{\underset{︸}{[\begin{matrix} 0 & 1 \\ e^{x_{1}} + x_{2}^{2} \cos x_{2} & 0 \end{matrix}]}} [\begin{array}{l} x_{1} \\ x_{2} \end{array}], \\ g (x, u) & = \underset{B (x, u)}{\underset{︸}{[\begin{matrix} 0 & 0 \\ x_{2} + x_{1} x_{2} \frac{u_{2}}{u_{1}} & 0 \end{matrix}]}} [\begin{array}{l} u_{1} \\ u_{2} \end{array}] . \end{array}$

B(x, u) in the second set of the SDC matrices is not ideal because it changes the problem into one with non-affine parameterization and singularity in the case of u₁ = 0. The limitation and proper selection of SDC matrices can be summarized in controllability assumption, as follows.

Assumption 1

The pair of {A_j(x_j(t)), B_j(x_j(t))} is a completely controllable parameterization of nonlinear affine-in-control system (1) for the states of the j-th agent x_j(t) ∈ ℝⁿ in t ∈ [0, ∞) [56].

The control objective is to minimize the cost function: 3 $J = \frac{1}{2} \int_{0}^{\infty} {x_{j}^{⊤} (t) Q_{j} (x_{j} (t)) x_{j} (t) + u_{j}^{⊤} (t) R_{j} (x_{j} (t)) u_{j} (t)} d t,$ where Q_j(x_j(t)) : ℝⁿ → ℝ^n×n and R_j(x_j(t)) : ℝⁿ → ℝ^m×m are weighting matrices of states and inputs of j-th agent, respectively. Q_j is symmetric positive semidefinite and R_j is symmetric positive definite.

Assumption 2

The pair of ${A_{j} (x_{j} (t)), Q_{j}^{1 / 2} (x_{j} (t))}$ is a completely observable parameterization of the non-linear affine-in-control system (1) for the states of j-th agent x_j(t) ∈ ℝⁿ in t ∈ [0, ∞), where $Q_{j}^{1 / 2}$ , in Eq. (3), is Cholesky decomposition of Q_j [57].

The control law of the SDRE is [56]: 4 $u_{j} (t) = - R_{j}^{- 1} (x_{j} (t)) B_{j}^{⊤} (x_{j} (t)) K_{j} (x_{j} (t)) x_{j} (t),$ where K_j(x_j(t)) is the positive-definite symmetric gain of the control law. The gain K_j(x_j(t)) is the solution to the state-dependent Riccati equation: 5 $A_{j}^{⊤} (x_{j} (t)) K_{j} (x_{j} (t)) + K_{j} (x_{j} (t)) A_{j} (x_{j} (t)) + Q_{j} (x_{j} (t)) - K_{j} (x_{j} (t)) B_{j} (x_{j} (t)) R_{j}^{- 1} (x_{j} (t)) B_{j}^{⊤} (x_{j} (t)) K_{j} (x_{j} (t)) = 0 .$

Without loss of generality, the input law (4) is transformed to: $u_{j} (t) = - R_{j}^{- 1} (x_{j} (t)) B_{j}^{⊤} (x_{j} (t)) K_{j} (x_{j} (t)) e_{j} (t),$ for controlling the system to a desired condition in which e_j(t) = x_j(t) – x_des,j(t) where x_des,j(t) is the desired state vector of j-th agent.

2.2.

Formation Control

The definition of formation control is introduced by a set of time-varying desired conditions for the agents to follow the leader of the group. The first agent, j = 1, is the leader and the rest of the agents are followers, j = 2, …, N. The leader is regulated (pointto-point motion) from an initial condition, x₁(0), to the desired condition, x_des,1(t_f), in which t_f is the final time of the control task. The schematic of the distribution of the agents, and the migration from the initial to the final condition, are illustrated in Fig. 1.

The rest of the agents follow the leader while they try to keep the predefined geometric formation during the task and sit at the exact predefined geometric formation at the end of the task: 6 $x_{des, j} (t) = x_{1} (t) - x_{f, 1} + x_{f, j}, j = 2, \dots, N,$ where x_f,j represents the desired coordinates (constant values) for the j-th agent in the predefined formation shape. Only the leader has a constant (timeinvariant) desired condition. Equation (6) shows that the time-varying desired conditions of the followers are fed by the online position of the leader and the constant desired formation. That means the followers try to keep the formation of the overall point-to-point trajectory of the leader from the initial to the final condition.

Remark 1

It should be pointed out that x_f,j, is a state vector, and only the first two or three states are different from zero (those that represent the configuration space). The rest of the desired states could be set to zero. Common dynamics of multirotor UAVs usually have six degrees of freedom that provide 12 states for the system. For the UAVs, only three position states of Cartesian coordinate (x, y, z) must be considered in the formation. The orientation and velocity states might be independently controlled.

2.3.

Collision Avoidance and Obstacle Avoidance Between Agents

The weighting matrix for the states of the SDRE controller can include nonlinear functions of state variables in its components. This property provides obstacle avoidance in regulation control based on artificial potential field method [52,58,59]. The distance between an obstacle and an agent is: 7 $d_{o, j} (t) = {‖ x_{p, j} (t) - x_{o} ‖}_{2},$ for j = 1, …, N where x_p,j(t) is a part of state vector, which includes the position coordinates, x_p,j(t) ∈ x_j(t), and has x_o as the position of the obstacle.

The next step is to define the distance between two consecutive agents for collision avoidance. The distance between the j-th and (j – 1)-th agent is: 8 $d_{j} (t) = {‖ x_{p, j} (t) - x_{p, j - 1} (t) ‖}_{2},$ for j = 2, …, N, it should be noted that the leader does not need to include collision avoidance terms. The collision avoidance approach in Eq. (8) only avoids the impact between two consecutive agents. The topology of the agents is defined through the distributed pattern (6) where the main character is the time-varying position of the leader. Evidently, this pattern and leader information does not provide enough safety in the migration of the multi-agent system. Hence, the collision avoidance term was introduced (8) to add another safety condition for migration without impact. This feature creates communication between two consecutive agents for reading and sending the position data.

The index is assigned to a given agent randomly; as a result, two consecutive agents might not be neighbors at the initial condition and could even be far away from each other. This initialization might raise issues. The collision avoidance between neighbor agents (8) might not be effective and collision might happen between j-th agent and other members. So, it is not assumed that agents j and j – 1 are closed; however, starting the regulation control, they follow the leader and the pattern, and they merge to find their neighbor agents in a short period. To avoid impact and collision in a global scheme, the following terms will be considered as well.

To include the collision avoidance between the f-th agent and all previous members, the following condition is introduced: 9 $d_{s, j} (t) = {\begin{cases} d_{t o t, k, j} (t), & d_{t o t, k, j} (t) < r_{\min}, \\ 1, & d_{t o t, k, j} (t) \geq r_{\min}, \end{cases}$ for j = 3, …, N, and k = 1, …, j – 1 where d_tot,k,j(t) = ||x_p,j(t) – x_p,k(t)||₂, is the distance between the j-th agent and the k-th agent, and r_min is the minimum collision avoidance safety distance. The pseudocode of Eq. (9) is presented as

for j = 3, …, N,

for k = 1, …, j – 1,

d_tot,k,j(t) = ∥x_p,j(t) – x_p,k(t)||₂,

if d_tot,k,j(t) < r_min,

d_s,j(t) = d_tot,k,j(t),

else,

d_s,j(t) = 1,

end,

end.

Equation (9) is updated at each time step of simulation when all the agents are moving from the initial to the final condition. In cases in which the j-th agent is close to one agent, the term d_s,j(t) is updated in (9); otherwise d_s,j(t) = 1. The different distance terms (d_o,j(t), d_j(t), and d_s,j(t)) will be used in Section 2.4 to set the obstacle and collision avoidance strategy based on the modification of Q_j(x_j(t)) matrix.

2.4.

Definition of Weighting Matrix

The state-dependent weighting matrices of the SDRE are one of the advantages of this method over a linear quadratic regulator (LQR) controller. The weighting matrices of the LQR are constant and cannot include parameters related to trajectories such as–the introduced distance terms in Section 2.3– for obstacle avoidance, collision avoidance, or adaptive designs. The definitions of the weighting matrices for the leader and followers are different in order to gain different motion speeds. The leader must be slower and the followers must be faster to track the leader properly. The following definition of the weighting matrix will provide for this condition.

The weighting matrix of states for the leader only includes obstacle avoidance: 10 $Q_{m, m, 1} (x_{1} (t)) = 1 + \frac{1}{d_{0, 1} (t)},$ where Q_m,m,1 (x₁ (t)) is the m-th diagonal component of Q₁(x₁(t)) related to m-th position state; m could be a state variable of leader’s dynamics, i.e. the second state, generated by Q_2,2,1(x₁(t)). It should be noted that only a few diagonal components of Q₁(x_j(t)) are responsible for obstacle avoidance. For the example presented in Remark 1, only the first two diagonal components of the weighting matrix include (10), and the rest of the 10 components are tuned conventionally. The UAV position variables are indeed defined in X, Y, Z, though the obstacle avoidance terms are set only for in X, Y; the results were satisfactory. The choice is with the designer, and in some cases, it can also be assigned for all three axes in 3D platforms.

Remark 2

If we notice the term $1 + \frac{1}{d_{o, 1} (t)}$ in (10), $\frac{1}{d_{o, 1} (t)}$ might gain a very small value if the distance between the obstacle and the leader is far. The constant term “1” in Q_m,m,1(x₁(t)) makes sure that the velocity of the leader is properly defined when $\frac{1}{d_{o, 1} (t)} \to 0$ .

The followers j-th must avoid obstacles, as well as the collision between the other agents, through the matrix Q_j(x_j(t)): 11 $Q_{m, m, j} (x_{j} (t)) = \frac{1}{d_{o, j} (t)} + \frac{1}{d_{j} (t)} + \frac{1}{d_{s, j} (t)},$ for j = 2, …, N. If the distance terms become small, which shows there is a probability of collision between agents or obstacles and an agent, the diagonal component of Q gets bigger and that state deviates the trajectory of the system (in other words, it avoids the collision). Assigning all the position states for an obstacle, or between collision avoidance speeds up all the states simultaneously and the collision avoidance fails. More details could be revisited in [52]. The following is helpful to clarify this point. Imagine a mobile robot moves forward in a straight line and there is an obstacle on the way. The actuators are two motors on the wheels. If the obstacle avoidance mechanism commands both wheels to speed up at the same time, the robot just moves faster in a straight line. In order to move the robot to one side, only one wheel must rotate faster! Therefore, avoidance functions are introduced only in one or two components related to the Cartesian coordinates.

In this section, three different distance terms were introduced: d_o,j(t), d_j(t), and d_s,j(t). It must be noted that these distance terms will be applied in the weighting matrix of states, using the artificial potential field method. The distance terms increase the value of the weighting matrix, leading the agent to avoid the collision.

Each distance term plays its own role in the migration of the multi-agent system. d_o,j(t) avoids collision of the agent with the obstacle in the environment; d_j(t) avoids collision of the two neighbor agents; and in case of the proximity of any previous agents to current agent, d_s,j(t) will be activated to avoid the collision. In summary, the three components will contribute to the safety in case of collision of the agents with the obstacle and with other agents.

The definition of the rest of the components of the weighting matrix of states, as well as the inputs, could be done based on conventional tuning methods and rules [60].

3.

System Dynamics and SDC Parameterization

3.1.

Unmanned Aerial Vehicle

The dynamic equation of one quadrotor unmanned aerial vehicle (UAV) is presented in this section. The schematic view of the system and the configuration of the rotors are illustrated in Fig. 2. The agents are identical and share the same dynamics. The generalized coordinates of one system are ${x_{c, j} (t), y_{c, j} (t), z_{c, j} (t), ϕ_{j} (t), θ_{j} (t), ψ_{j} (t)}$ , where $ξ_{1, j} (t) = {[x_{c, j} (t), y_{c, j} (t), z_{c, j} (t)]}^{⊤} (m)$ defines the center-of-mass (CoM) of the j-th UAV in the Cartesian coordinate system, and $ξ_{2, j} (t) = {[ϕ_{j} (t), θ_{j} (t), ψ_{j} (t)]}^{⊤} (rad)$ are its Euler angles. The linear and angular velocities of the drone are presented υ_1,j(t) = [u_j(t), v_j(t), w_j(t)]^⊤(m/s) and υ_2,j(t) = [p_j(t), q_j(t), r_j(t)]^⊤(rad/s), respectively. The state vector of the system is defined as: 12 $x_{j} (t) = {[ξ_{1, j}^{⊤} (t), ξ_{2, j}^{⊤} (t), υ_{1, j}^{⊤} (t), υ_{2, j}^{⊤} (t)]}^{⊤} .$

The kinematics relations 13 $\begin{array}{l} {\dot{ξ}}_{1, j} (t) & = R_{Z Y X, j} (ξ_{2, j} (t)) υ_{1, j} (t), \\ {\dot{ξ}}_{2, j} (t) & = T_{j} (ξ_{2, j} (t)) υ_{2, j} (t), \end{array}$ are held for the quadrotor system. Considering that the quadrotor system flies based on small changes in the Euler angles, it is not supposed to perform aggressive and aerobatic maneuvers, and the hovering state is assumed for the flight condition.

Therefore, the derivative of state-vector (12) generates ${\dot{x}}_{j} (t) = {[{\dot{ξ}}_{1, j}^{⊤} (t), {\dot{ξ}}_{2, j}^{⊤} (t), {\dot{υ}}_{1, j}^{⊤} (t), {\dot{υ}}_{2, j}^{⊤} (t)]}^{⊤} .$

Assuming small changes in Euler angels during flight control, the derivative of the kinematics equation (13): $\begin{array}{l} {\ddot{ξ}}_{1, j} (t) & = \underset{\approx 0}{\underset{︸}{{\dot{R}}_{Z Y X, j} (ξ_{2, j} (t))}} υ_{1, j} (t) + \underset{\approx I}{\underset{︸}{R_{Z Y X, j} (ξ_{2, j} (t))}} {\dot{υ}}_{1, j} (t), \\ {\ddot{ξ}}_{2, j} (t) & = \underset{\approx 0}{\underset{︸}{{\dot{T}}_{j} (ξ_{2, j} (t))}} υ_{2, j} (t) + \underset{\approx I}{\underset{︸}{T_{j} (ξ_{2, j} (t))}} {\dot{υ}}_{2, j} (t), \end{array}$

shows that in a hovering state, ${\ddot{ξ}}_{1, j} (t) \approx {\dot{υ}}_{1, j} (t)$ and ${\ddot{ξ}}_{2, j} (t) \approx {\dot{υ}}_{2, j} (t)$ are valid. As a result, the state-space representation of one agent is [61]: 14 ${\dot{x}}_{j} = [\begin{matrix} R_{Z Y X, j} (ξ_{2, j}) υ_{1, j} \\ T_{j} (ξ_{2, j}) υ_{2, j} \\ \frac{1}{m_{j}} {R_{Z Y X, 3, j} (ξ_{2, j}) T_{B, j} - m_{j} g e_{3} - D_{f, j} {\dot{ξ}}_{1, j}} \\ J_{j}^{- 1} (ξ_{2, j}) {τ_{B, j} - C_{j} (ξ_{1, j}, ξ_{2, j}) {\dot{ξ}}_{2, j}} \end{matrix}],$ where m_j(kg) is the mass of the drone, T_B,j(t)(N) is the thrust force in Z_c direction (upward), g = 9.81(m/s²) is gravity constant, e₃ = [0, 0, 1]^⊤, D_f,j = diag(D_x, D_y, D_z)(kg/s) includes drag, aerodynamic effects, etc., and $τ_{B, j} (t) = {[τ_{ϕ, j} (t), τ_{θ, j} (t), τ_{ψ, j} (t)]}^{⊤}$ (Nm) is the input torque vector. Moreover, the details of $R_{Z Y X, j} (ξ_{2, j} (t)), T_{j} (ξ_{2, j} (t)), C_{j} (ξ_{1, j} (t), ξ_{2, j} (t))$ , and $J_{j} (ξ_{2, j} (t))$ can be found in Section III of Ref. [61], and $R_{Z Y X, 3, j} (ξ_{2, j} (t))$ is the third column of the rotation matrix $R_{Z Y X, j} (ξ_{2, j} (t))$ .

The quadrotor UAV is underactuated, and the control structure uses a cascade design, which controls the translation part of dynamics (14) in the first layer and the orientation in the second layer. The translation dynamics includes the states x_t,j(t) = [x_c,j(t), y_c,j(t), z_c,j(t), u_j(t), v_j(t), w_j(t)]^⊤, expressed in rows 1-3 and 7-9 of Eq. (14). The SDC parameterization considers full actuation for translation control: $\begin{array}{l} A_{t, j} (x_{t, j} (t)) & = [\begin{matrix} 0_{3 \times 3} & R_{Z Y X, j} (ξ_{2, j} (t)) \\ 0_{3 \times 3} & - \frac{1}{m_{j}} D_{f, j} R_{Z Y X, j} (ξ_{2, j} (t)) \end{matrix}], \\ B_{t, j} & = [\begin{matrix} 0_{3 \times 3} \\ \frac{1}{m_{j}} I_{3 \times 3} \end{matrix}] . \end{array}$

The component ${[R_{Z Y X, 3, j} (ξ_{2, j}) T_{B, j}]}_{3 \times 1}$ in Eq. (14) is a vector of thrust, and is distributed in the translation dynamics by the third component of the rotation matrix. Besides the SDRE design, there is a cascade controller as well, Eq. (17), to address the underactuation condition of the quadrotor.

In order to implement the cascade design and SDRE control, in the first step, the SDRE controller must be operated with the assumption of full actuation condition, which demands a $B_{t, j} = [\begin{matrix} 0_{3 \times 3} \\ \frac{1}{m_{j}} I_{3 \times 3} \end{matrix}]$ matrix of 6 × 3 dimension. If the third component of the rotation matrix enters the SDC parameterization, a scalar thrust term is obtained as the single input to the problem, which creates an obstacle in the design. This topic has been presented in previous implementations and exercised with different conditions, platforms, and controllers [61–63], but it was originally proposed as cascade control design; details can be seen in Ref. [64].

Assigning weighting matrices with proper dimensions and solving the SDRE for the translation part: $A_{t, j}^{⊤} (x_{t, j}) K_{t, j} (x_{t, j}) + K_{t, j} (x_{t, j}) A_{t, j} (x_{t, j}) - K_{t, j} (x_{t, j}) B_{t, j} R_{t, j}^{- 1} (x_{t, j}) B_{t, j}^{⊤} K_{t, j} (x_{t, j}) + Q_{t, j} (x_{t, j}) = 0,$ results in the sub-optimal gain K_t,j(x_t,j(t)) for the control law: 15 $u_{t, j} (t) = - R_{t, j}^{- 1} (x_{t, j} (t)) B_{t, j}^{⊤} K_{t, j} (x_{t, j} (t)) e_{t, j} (t),$ where e_t,j(t) = x_t,j(t) – x_t,des,j(t) in which x_t,des,j(t) includes the desired translation variables. The three components of the control law, in (15), are substituted in the thrust equation (16) [61]: 16 $T_{B, j} (t) = m_{j} {R_{Z Y X, 3, j}^{⊤} (ξ_{2, j} (t)) (u_{t, j} (t) + g e_{3})} .$

Note that gravity is compensated for in (16), and for this reason, the SDC matrix A_t,j(x_t,j(t)) excluded the gravity term. The second cascade control layer regulates the orientation of the drone to help the translation part, based on two constraints [62]: 17 $\begin{array}{l} θ_{des, j} (t) & = \tan^{- 1} \frac{u_{t, 1, j} (t) \cos (ψ_{des, j}) + u_{t, 2, j} (t) \sin (ψ_{des, j})}{u_{t, 3, j} (t) + g}, \\ ϕ_{des, j} (t) = & \sin^{- 1} \frac{u_{t, 1, j} (t) \sin (ψ_{des, j}) - u_{t, 2, j} (t) \cos (ψ_{des, j})}{\sqrt{u_{t, 1, j}^{2} (t) + u_{t, 2, j}^{2} (t) + {(u_{t, 3, j} (t) + g)}^{2}}}, \end{array}$ where ψ_des,j (rad) is the desired constant yaw of the system and i.e. u_t,2,j(t) is the second component of u_t,j(t). The state vector of orientation control is set as $x_{0, j} (t) = {[ϕ_{j} (t), θ_{j} (t), ψ_{j} (t), p_{j} (t), q_{j} (t), r_{j} (t)]}^{⊤}$ , and the relevant dynamics for orientation, rows 4–6 and 10–12 of (14), provides the SDC matrices: $\begin{array}{l} A_{0, j} (x_{0, j}) & = [\begin{matrix} 0_{3 \times 3} & T_{j} (ξ_{2, j}) \\ 0_{3 \times 3} & - J_{j}^{- 1} (ξ_{2, j}) C_{j} (ξ_{1, j}, ξ_{2, j}) \end{matrix}], \\ B_{o, j} (x_{0, j}) & = [\begin{matrix} 0_{3 \times 3} \\ J_{j}^{- 1} (ξ_{2, j}) \end{matrix}] . \end{array}$

The SDRE for the orientation control is: $A_{o, j}^{⊤} (x_{o, j}) K_{o, j} (x_{o, j}) + K_{o, j} (x_{o, j}) A_{o, j} (x_{o, j}) - K_{o, j} (x_{o, j}) B_{o, j} (x_{o, j}) R_{o, j}^{- 1} (x_{o, j}) B_{o, j}^{⊤} (x_{o, j}) K_{0, j} (x_{0, j}) + Q_{0, j} (x_{0, j}) = 0,$ which generates the sub-optimal gain K_o,j(x_o,j(t)) for the control law: $u_{0, j} (t) = - R_{o, j}^{- 1} (x_{0, j} (t)) B_{o, j}^{⊤} (x_{0, j} (t)) K_{o, j} (x_{0, j} (t)) e_{o, j} (t),$ where e_o,j(t) = x_o,j(t) – x_o,des,j(t) and x_o,des,j(t) includes the desired orientation variables. Note that ϕ_des,j(t) and θ_des,j(t) in x_o,des,j(t) are updated in each time-step of the simulation or experiment by (17).

The input thrust and torque vector to the system are generated by four rotating propellers: 18 $[\begin{array}{l} T_{B, j} (t) \\ τ_{ϕ, j} (t) \\ τ_{θ, j} (t) \\ τ_{ψ, j} (t) \end{array}] = M_{x, j} ω_{j}^{2} (t),$ where ω_j(t) ∈ ℝ⁴ is the angular velocity vector of the rotors and $M_{x, j} = [\begin{matrix} k & k & k & k \\ 0 & - l k & 0 & l k \\ - l k & 0 & l k & 0 \\ d & - d & d & - d \end{matrix}],$ is the mixer matrix of the quadrotor, which is defined based on the “plus” configuration of the rotors in Fig. 2, in which k(Ns²/rad²) is the lift constant or thrust factor, l(m) is the distance between motor and CoM of the quadrotor, and d(Nms²/rad²) is the drag constant.

Remark 3

To regulate the system (2) to the desired condition, the equilibrium point of the system should be zero. Many systems, such as robotic manipulators, multirotor UAVs, etc., work under the effect of gravity, where f_j(0) ≠ 0. For those systems, the gravity compensation should be considered as shifting the equilibrium to zero [65]. Here in Section 3.1, the gravity will be compensated for through the cascade control design of the UAVs, which performs translation control in the first layer and orientation control in the second one [62].

3.2.

Wheeled Mobile Robot

Consider j-th wheeled mobile robot of a multiagent system with differential wheels working on {X, Y} planar Cartesian coordinates. The coordinates, and global and local axes, are shown in Fig. 3. The center of mass of the robot is defined by {x_c,j(t), Y_c,j(t)} (m), and the rotation of the robot around Z axis, measured from X axis (counter-clockwise indicates positive direction), is called ϕ_j(t)(rad). The system is subjected to two constraints showing that the robot cannot move alongside the axis of the wheels and that the wheels only have a pure rolling motion (rolling without slipping).

The first auxiliary relation is found by combining the constraints on the left and right wheels, ${\dot{ϕ}}_{j} (t) = \frac{r}{2 b} ({\dot{θ}}_{r, j} (t) - {\dot{θ}}_{l, j} (t))$ . This integration generates a holonomic constraint [66]: 19 $ϕ_{j} (t) = \frac{r}{2 b} (θ_{r, j} (t) - θ_{l, j} (t)) .$

The non-holonomic constraint, derived from the rolling condition of the wheels, is obtained [67]: 20 ${\dot{x}}_{c, j} (t) \cos ϕ_{j} (t) + {\dot{y}}_{c, j} (t) \sin ϕ_{j} (t) = \frac{r}{2} ({\dot{θ}}_{r, j} (t) + {\dot{θ}}_{l, j} (t)) .$

Equation (19) is the holonomic constraint of the system and equations $({\dot{y}}_{c, j} (t) \cos ϕ_{j} (t) - {\dot{x}}_{c, j} (t) \sin ϕ_{j} (t) = 0)$ and (20) are non-holonomic constraints, which could be represented as $A_{j} (q_{j} (t)) {\dot{q}}_{j} (t) = 0$ , where 21 $A_{j} (q_{j} (t)) = [\begin{matrix} - \sin ϕ_{j} (t) & \cos ϕ_{j} (t) & 0 & 0 \\ - \cos ϕ_{j} (t) & - \sin ϕ_{j} (t) & \frac{r}{2} & \frac{r}{2} \end{matrix}],$ in which generalized coordinate vector of the j-th system is: 22 $q_{j} (t) = {[x_{c, j} (t), y_{c, j} (t), θ_{r, j} (t), θ_{l, j} (t)]}^{⊤} .$

The dynamics equation of the mobile robot is the common form ofa second-order differential equation: 23 $M_{j} (q_{j} (t)) {\overset{..}{q}}_{j} (t) + c_{j} (q_{j} (t), {\dot{q}}_{j} (t)) = E τ_{j} (t) - A_{j}^{⊤} (q_{j} (t)) λ_{j} (t),$ where M_j(q_j(t)) : ℝ⁴ → ℝ^4×4 is the inertia matrix; the Coriolis-centrifugal vector is $c_{j} (q_{j} (t), {\dot{q}}_{j} (t)) : ℝ^{4} \times ℝ^{4} \to ℝ^{4}; E = {[0_{2 \times 2}, I_{2 \times 2}]}^{⊤}, τ_{j} (t) = {[τ_{r, j} (t), τ_{l, j} (t)]}^{⊤}$ is the input torque vector to the wheels; A_j(q_j(t)) : ℝ⁴ → ℝ^2×4 is constraint matrix (21); and λ_j(t) = [λ_1,j(t),λ_2,j(t)]^⊤ is the vector of the Lagrange multipliers. The details of the inertia matrix and Coriolis vector can be found in Ref. [67]. One way to omitting unknown Lagrange multipliers from Eq. (23) is to find the null space of A_j(q_j(t)): $S_{j} (q_{j} (t)) = [\begin{matrix} \frac{r}{2} \cos (ϕ_{j} (t)) & \frac{r}{2} \cos (ϕ_{j} (t)) \\ \frac{r}{2} \sin (ϕ_{j} (t)) & \frac{r}{2} \sin (ϕ_{j} (t)) \\ 1 & 0 \\ 0 & 1 \end{matrix}]$ and pre-multiply the transpose of S_j(q_j(t)) to (23) to obtain [66]: 24 $S_{j}^{⊤} (q_{j} (t)) M_{j} (q_{j} (t)) {\overset{..}{q}}_{j} (t) + S_{j}^{⊤} (q_{j} (t)) c_{j} (q_{j} (t), {\dot{q}}_{j} (t)) = τ_{j} (t) .$

The angular velocities of the wheels are arranged in a vector $v_{j} (t) = {[{\dot{θ}}_{r, j} (t), {\dot{θ}}_{l, j} (t)]}^{⊤}; {\dot{q}}_{j} (t) = S_{j} (q_{j} (t)) v_{j} (t)$ is introduced; the second time derivative of the generalized coordinates results in: 25 ${\overset{..}{q}}_{j} (t) = {\dot{S}}_{j} (q_{j} (t)) v_{j} (t) + S_{j} (q_{j} (t)) {\dot{v}}_{j} (t) .$

The state-vector of the system is chosen as $x_{j} = {[q_{j}^{⊤} (t), v_{j}^{⊤} (t)]}^{⊤}$ . Substitution of ${\overset{..}{q}}_{j} (t)$ from (24) into (25) will result in a state-space representation of the system: 26 ${\dot{x}}_{j} = [\begin{matrix} S_{j} v_{j} \\ - {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} [S_{j}^{⊤} M_{j} {\dot{S}}_{j} v_{j} + S_{j}^{⊤} c_{j}] \end{matrix}] + [\begin{matrix} 0_{4 \times 2} \\ {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} \end{matrix}] τ_{j} .$

In the control scheme of WMR, the controlled output of the system is a point P_w(t), ahead of the CoM of the mobile robot; w_s is the distance between the CoM of the robot; P_w(t) (see Fig. 3). The point is related to the state variables of the WMR through a kinematics relation. To control this point, the outputand statedependent Riccati equation (OSDRE) approach is needed [57]. The output of the j-th WMR is defined as $y_{j} (t) = [\begin{array}{l} x_{c, j} (t) + w_{s} \cos ϕ_{j} (t) \\ y_{c, j} (t) + w_{s} \sin ϕ_{j} (t) \end{array}],$ and the first two components of S_j(q_j(t))v_j(t) results in 27 $\begin{array}{l} {\dot{x}}_{c, j} (t) & = \frac{r}{2} ({\dot{θ}}_{r, j} (t) + {\dot{θ}}_{l, j} (t)) \cos ϕ_{j} (t), \\ {\dot{y}}_{c, j} (t) & = \frac{r}{2} ({\dot{θ}}_{r, j} (t) + {\dot{θ}}_{l, j} (t)) \sin ϕ_{j} (t) . \end{array}$

Substituting (27) into ${\dot{y}}_{j} (t)$ , and another time derivative ${\overset{..}{y}}_{j} (t)$ results in: 28 ${\overset{..}{y}}_{j} (t) = α_{j} (x_{j} (t)) {\dot{v}}_{j} (t) + β_{j} (x_{j} (t)),$ where α_j(x_j(t)) : ℝ⁶ → ℝ^2×2 and β_j(x_j(t)) : ℝ⁶ → ℝ². α_j(x_j(t)) and β_j(x_j(t)) are a nonlinear matrix and vector, respectively, that have been generated during the computation of the derivative of the output. The symbolic representation of them can be revisited in [67]. α_j(x_j(t)) and β_j(x_j(t)) are kinematic transformations that change the dynamics of the mobile robot from a state to an output representation. Substituting ${\dot{v}}_{j} (t)$ from the lower set of (26) into (28) provides: 29 ${\overset{..}{y}}_{j} = α_{j} {- {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} [S_{j}^{⊤} M_{j} {\dot{S}}_{j} v_{j} + S_{j}^{⊤} c_{j}]} + β_{j} + α_{j} {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} τ_{j} .$

By defining a new state-vector for the output dynamics (29), $z_{j} (t) = {[y_{j}^{⊤} (t), {\dot{y}}_{j}^{⊤} (t)]}^{⊤}$ , one could find the state-space representation of the output dynamic: 30 ${\dot{z}}_{j} = [\begin{matrix} {\dot{y}}_{j} \\ α_{j} {- {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} [S_{j}^{⊤} M_{j} {\dot{S}}_{j} v_{j} + S_{j}^{⊤} c_{j}]} + β_{j}] \end{matrix} + [\begin{matrix} 0_{2 \times 2} \\ α_{j} {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} \end{matrix}] τ_{j.}$

The outputand state-dependent coefficient parameterization of system (30) is expressed as [57]: $\begin{matrix} A_{j} = [\begin{matrix} 0_{2 \times 2} & I_{2 \times 2} \\ 0_{2 \times 2} & [α_{j} {- {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} [S_{j}^{⊤} M_{j} {\dot{S}}_{j} v_{j} + S_{j}^{⊤} c_{j}} + β_{j}] {\dot{y}}_{j}^{†} \end{matrix}], \\ B_{j} = [\begin{matrix} 0_{2 \times 2} \\ α_{j} {(S_{j}^{⊤} M_{j} S_{j})}^{- 1} \end{matrix}], \end{matrix}$ where ${\dot{y}}_{j} (^{t) †}$ is pseudo-inverse of the velocity of the output ${\dot{y}}_{j} (t)$ . Finally, the control law is set in the form of: $u_{j} (t) = - R_{j}^{- 1} (x_{j} (t)) B_{j}^{⊤} (x_{j} (t)) K_{j} (x_{j} (t)) [z_{j} (t) - z_{des, j} (t)] .$

The necessary conditions of the output, controllability, and observability of the OSDRE, and condition on ${\dot{y}}_{j} (^{t) †}$ to avoid a singularity issue, were reported in [57].

4.

Simulations

4.1.

Multi-agent systems of UAV

A multi-agent system of N = 1,089 quadrotor agents is considered for a simulation that represents a regulation case from an initial to a final condition in t_f = 20(s). The agents must form a square shape when they are launched and keep the shape during the regulation. As they approach the final condition, the square shape and arrangement of the agents get more precise, and the error in regulation also converges to zero. The pseudocode of the square pattern of agents for the simulation case of 1,089 drones. The following loop will define the local arrangement of the agents in a square shape during the regulation, presented in the pseudocode:

for k = 1, …,N_n,

for j = 1+(k – 1)N_n,…,N_n + (k – 1)N_n,

$x_{c, f, j} = \frac{- (L + D_{a})}{2} + j D_{a} - (k - 1) L$ ,

$y_{c, f, j} = \frac{- (L + D_{a})}{2} + k D_{a}$ ,

z_c,f,j = 0,

end,

where length of square side is denoted by L = 30(m), N_n = 33 is number of rows and columns, and

D_{a} = \frac{L}{N_{n}} = 0.9091 (m)

.

The leader is the first agent of the system and guides the multi-agent system in point-to-point regulation from an arbitrary start to an endpoint. The 1,089 agents are initially scattered randomly inside a 10 × 10(m) square at the height of 10(m). The initial orientation angle and velocities of the translation and orientation states of the agents are zero. The three components of the pseudocode (x_c,f,j, y_c,f,j, z_c,f,j) are the first set of vector x_f,j ∈ ℝ¹² in (6) for the quadrotor multi-agent system.

The position of the obstacle is (x_o, y_o, z_o) = (15, 0.25, 6)(m), with a safety radius of 0.5 for collision avoidance. The final position coordinates for the leader are: $x_{f, 1} = [\begin{matrix} x_{c, des, 1} \\ y_{c, des, 1} \\ z_{c, des, 1} \\ 0_{9 \times 1} \end{matrix}] = [\begin{matrix} 20 + x_{c, f, 1} \\ y_{c, f, 1} \\ z_{c, f, 1} \\ 0_{9 \times 1} \end{matrix}],$ which is constant though the value for the rest of the agents is time-varying, presented in Eq. (6). This is a result of following the leader.

The physical characteristics of the agents are similar in the following way. The distance between the motor and the CoM of the quadrotor is l = 0.225(m); the radius of the propeller is R = 0.075(m); the lift constant is k = 2.98 × 10⁻⁶(Ns²rad⁻²); the drag constant is d = 1.14 × 10⁻⁷(Nms²rad⁻²); additionally I_xx = 4.856×10⁻³(kgm²), I_yy = 4.856×10⁻³(kgm²), and I_zz = 8.801 × 10⁻³(kgm²) are moments of inertia around x, y, z axes. The mass of the drone is m = 0.468(kg), the aerodynamics matrix is D_f = diag(0.25, 0.25, 0.25)(kg/s), the minimum and maximum rotor velocities are also ω_min,max = 496.48, 744.73(rad/s).

The motion of the multi-agent system is represented in Fig. 4. Because the agents are very close to each other at the beginning, and the collision avoidance scheme in tuning increases the Q matrix drastically, we see a rapid divergence at the commencement of the simulation. The beginning of the regulation is like an “explosion” of multiple agents escaping from collisions with each other. The Cartesian-coordinate errors of the agents in 3D positioning are illustrated in Fig. 5 for all agents, followed by a zoomed view that shows the leader first converging towards zero, then the followers catching up with the leader.

The weighting matrices of the leader are selected as: $\begin{array}{l} R_{t, 1} = 10 \times I_{3 \times 3}, R_{0, 1} = I_{3 \times 3}, \\ Q_{t, 1} (x) = 0.5 \times diag (1 + \frac{1}{d_{0, 1} (x)}, 1 + \frac{1}{d_{0, 1} (x)}, 1, 2, 2, 5), \\ Q_{0, 1} = diag {(2_{1 \times 3}, 0 . 5_{1 \times 3})}_{6 \times 6} . \end{array}$

The weighting matrices of the followers are: $\begin{array}{l} R_{t, j} & = I_{3 \times 3}, R_{0, j} = I_{3 \times 3}, \\ Q_{t, j} (x) & = d i a g (1 + Q_{D} (x), 1 + Q_{D} (x), 1, 2, 2, 5) \\ Q_{o, j} & = d i a g {(2_{1 \times 3}, 0 . 5_{1 \times 3})}_{6 \times 6}, \end{array}$ where $Q_{D} (x (t)) = \frac{1}{d_{j} (x (t))} + \frac{1}{d_{o, j} (x (t))} + \frac{1}{d_{s, j} (x (t))}$ , with d_s,j(x(t)) set by (9).

Obstacle avoidance by embedding an artificial potential field into the SDRE has been reported on and analyzed in the literature. It has been stated that this method does not guarantee absolute obstacle avoidance, though it reduces the chance of impact significantly [52]. Figure 6 shows that the minimum distance between some agents and the obstacle was 30(cm). The collision avoidance between the agents is reported in Fig. 7. A group of drones in Fig. 7 merged to almost a 30(cm) distance; those are the ones at the beginning of rows, with the previous agents located at the end of previous rows. For the rest of the agents, the collision avoidance term, d_j(t), shows a collision when it is lower than 0.48(m). The reason for the collision between the agents is the accumulation of 1, 089 drones in a 10 × 10(m) at the initial condition of the simulation. After the commencement of the simulation, the collision avoidance scheme guided the drones and increased the distance between them after 8 seconds.

Checking the collision between two consecutive agents does not necessarily show all the collisions. Other agents might also have collisions that are demonstrated in Fig. 8.

Contrary to d_j(t), the parameter d_s,j(t) checks the collision of j-th agent with all previous agents. The behavior of d_s,j(t) is similar to d_j(t); at the beginning, the number of collisions is higher, and moving towards the final pattern, the distance between the agents gets bigger, starting from first time-step and ending 30,000-th time-steps, discretized into 20 seconds of total simulation time (samples are shown in the axes instead of time as a result of the mesh presentation). The collision of the agents is obvious since their arrangement is quite ambitious, leaving only 42(cm) between two agents in the final square shape. It can be observed in Fig. 7 when the distance for some drones is less than 0.48(m).

The simulation of 1, 089 agents, presented in Fig. 4 might not reveal the complexity of the multi-agent system dynamics, therefore, the data of the simulation for the 5th agent is presented. The position and velocity states of the 5th drone are presented in Figs. 9 and 10, respectively. The equation of the mixer matrix (18) presents the input angular velocities of the rotors, limited to the lower and upper bounds, Fig. 11.

4.2.

WMR

A group N = 45, agents was chosen for a regulation simulation over t_f = 15(s) with circular final arrangements of the agents. The number of agents was selected smaller in comparison with Section 4.1 to highlight the role of obstacle avoidance in the results. The following algorithm defines the arrangement of the agents in a circular pattern:

for j = 1, …, N,

r_p = 1,

c_o = r_p,

$x_{c, f, j} = x_{shift} + r_{p} \cos \frac{(j - 1) 2 π}{N_{c}},$ ,

$y_{c, f, j} = y_{shift} + r_{p} \sin \frac{(j - 1) 2 π}{N_{c}},$ ,

for k = 1, …, N_r,

$if j > \frac{k (k + 1) N_{c}}{2},$

c_o = c_o + r_p,

$x_{c, f, j} = x_{shift} + r_{p} \cos \frac{(j - N_{c} - 1) 2 π r_{p}}{N_{c} c_{o}},$

$y_{c, f, j} = y_{shift} + r_{p} \sin \frac{(j - N_{c} - 1) 2 π r_{p}}{N_{c} c_{o}}$ ,

end,

where (x_shift, y_shift) = (40, 20)(m) defines the center of the circular pattern; N_c = 3 is the number of agents in the first circle; N_r = 5 is number of rings;

N = \frac{N_{r} (N_{r} + 1) N_{c}}{2} = 45

is total number of agents; and r_p = 1(m) is the distance between the rings.

The position of the obstacle is set as (x_o, y_o) = (25, 10)(m) with a safety radius of x_min = 0.5(m). The agents are randomly scattered in a 10 × 10(m) square at the origin. These random scattered agents are almost 40(m) away from the desired position of the pattern.

The physical characteristics of the mobile robots were chosen as follows: the radius of the wheels is denoted by r = 0.08(m); the length of the wheels’ axis by b = 0.145(m); mass of base is m_b = 6(kg); mass of wheels is m_w = 0.32(kg); moment of inertia components are I_zz,b = 0.06363,I_zz,w = 0.00006, I_yy,w = 0.00081(kgm²); and the distance of the look-ahead control term is selected as w_s = 0.1(m).

The weighting matrices of the leader are set: $\begin{array}{l} R_{1} & = I_{2 \times 2}, \\ Q_{1} (x (t)) & = diag (0.5, 0.5 + \frac{1}{\sqrt{{(d_{0, 1} (x (t)) - 2)}^{2}}}, 2, 2) \end{array} .$

The weighting matrices of the followers are selected as: $R_{j} = I_{2 \times 2}, Q_{j} (x (t)) = diag (0.5, 0.5 + Q_{D} (x (t)), 2, 2),$ where $Q_{D} (x (t)) = \frac{1}{d_{j} (x (t))} + \frac{1}{d_{s, j} (x (t))} + \frac{1}{\sqrt{{(d_{0, 1} (x (t)) - 2)}^{2}}},$ in which $d_{j} = \sqrt{{(x_{c, j} - x_{c, j - 1})}^{2} + {(y_{c, j} - y_{c, j - 1})}^{2}}$ .

The trajectories of the leader and follower agents are illustrated in Fig. 12. The leader and the followers avoided the obstacle, which are represented in the trajectories. The error convergence of the multi-agent system is presented in Fig. 13. The obstacle and collision avoidance terms are plotted in Figs. 14 and 15, respectively. The input signals and state variables of one agent are presented in Figs. 16 and 17, respectively. The error of the leader and 5th follower is found at 122.6741(mm) and 411.7172(mm) respectively. It is important to notice that the error is a function of the look-ahead control parameter which was set w_s = 0.1(m). This means that the robot is regulating the look-ahead point towards the desired condition, not the CoM of the robot. If w_s is reduced, the error changes. The overshoot and the behavior of the regulation also change as the look-ahead parameter varies. A comparative table is presented to clarify this effect in Table 2. Another simulation was done for 1, 050 WMRs to show the capacity of the method for a largepopulation multi-agent system. For the sake of brevity, only the migration of the systems from initial to final condition, following the leader and the pattern, is illustrated in Fig. 18.

Table 2.

The error of the system under different w_s values.

w_s(mm)	error leader (mm)	error 5th follower (mm)
20	144.3912	911.3158
50	76.4762	345.8130
80	104.3379	388.0787
100	124.3148	427.7242
200	224.7471	630.6206
300	324.6205	827.9986
400	424.4897	1025.9861
500	524.4070	1223.8401

4.3.

Tuning of Weighting Matrices

There are some general rules for the selection of weighting matrices of the SDRE. Q penalizes the states and R penalizes the inputs. Since the inverse of R is the first term of the control law, a decrease in R increases the input signal, which results in more energy consumption. Increase in Q also increases the precision of the state regulation. In this work, the weighting matrix of states has obstacle/collision avoidance terms for multi-agent control design. The definition of the weighting matrices for the leader and followers must be different. If the leader is too fast, the followers might not follow the leader with proper accuracy.

Hence, the diagonal components of Q relative to the position coordinates of the leader are smaller than the ones for the followers for both simulations of the UAVs and WMRs. Another important point for tuning the robots is the contrast between the weighting matrices’ value and the distance terms! If large values are selected for Q, the changes due to distance terms vanish, and collision/obstacle avoidance terms are not effective. These points were considered i n the tuning of the matrices in the simulation sections. Considering the aforementioned point, it must be said that several trials and errors are needed to tune the controller.

5.

Conclusion

This work presents an SDRE control design (without augmentation of other techniques) for the formation regulation of a multi-agent system of robotic systems that takes into account the complex dynamics of the agents. The command to the multi-agent system is given to the leader, and the followers pursue the leader in accordance with the predefined pattern of the multi-agent system. Obstacle and collision avoidance among the agents was designed through the modification of the weighting matrices of states and the artificial potential field method. The application of the proposed method was devoted to the control of unmanned multi-copter systems and differential wheeled mobile robots. A highly populated multi-agent system of 1,089 agents was simulated for the UAV, and another simulation was presented for the WMRs with 45, and then 1,050 agents which were successfully regulated in the control task, preserving the expected formation shapes.

Język:: Angielski

Częstotliwość wydawania:: 4 razy w roku
Dziedziny czasopisma:: Informatyka, Sztuczna inteligencja, Inżynieria, Elektrotechnika, Inżynieria sterowania, metrologia i testowanie, Inżynieria mechaniczna, Podstawy inżynierii mechanicznej

Kanał RSS czasopisma

Formation Control of Multi-agent Nonlinear Systems using The State-Dependent Riccati Equation

Saeed Rafee Nekoo

Data publikacji: 31 mar 2025

Zakres stron: 17 - 32

Otrzymano: 10 sty 2024

Przyjęty: 23 lip 2024

DOI: https://doi.org/10.14313/jamris-2025-003

Słowa kluczoweMulti-agent, Formation, Control, SDRE, Robotics

© 2025 Saeed Rafee Nekoo, published by Sciendo

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.

Słowa kluczowe
Multi-agent, Formation, Control, SDRE, Robotics