\section{Introduction}
\IEEEPARstart{I}{ntelligent} mobile robots are widely used to replace humans in performing complex, monotonous, and dangerous tasks due to their compactness, flexibility, and ability to carry a variety of sensors.  Nowadays, with the increasing complexity of operational tasks in the fields of warehousing, logistics, agriculture, subsea exploration, etc., there are apparent productivity and efficiency bottlenecks to operate with a single robot. For example, a single robot has limited perception capability and operating range. Multiple robots cooperating in the same space have advantages in efficiency, robustness, and task completion rate.

Motion planning technology is critical to realize the autonomous collaborative operation of multiple robots.  The general description of multi-robot cooperative motion planning is: in the same space,  multiple robots need to start from different initial points to reach the corresponding target points.  The motion planner is required to plan a set of collision-free and dynamic-feasible trajectories with minimum time or energy consumption. Achieving such task with various interaction patterns is quite challenging, especially in dense environments where there exist multiple active agents under non-communication and local observable conditions. Most of the time, we have to consider dynamic constraints of different types of robots, including velocity, acceleration, curvature, jerk, snap, etc., so that trajectories are truly executable.

\subsection{Related Works}
Multi-robot motion planning approaches can be divided into centralized methods and Decentralized methods.  Centralized methods contain optimization-based trajectory generation methods such as \cite{mellinger2012mixed,augugliaro2012generation}, and heuristic search-based planning methods such as \cite{yu2016optimal,svancara2017new}. Centralized methods have the advantage of completeness or probabilistic completeness.  However, this type of method requires obtaining global state information and cannot handle local observable situations.  In addition, centralized methods are constantly suffering from scalability issues.  As the number of robots in the task scenario increases, the computation complexity would rise exponentially, which is pretty challenging for the central computational device.

\begin{figure}[!t]
	\centering
	\includegraphics[width=3.5 in]{Tot_MP_Arch.png}
	\caption{The architecture of the proposed MASAC-based hybrid motion planner. The CTDE RL paradigm-based front-end of this planner is responsible for searching the joint discrete waypoints. The minimal snap with constraints optimizer-based back-end of the planner is responsible for generating the final motion trajectory.}
	\label{Fig_TrajOptim}
\end{figure}

In contrast, decentralized methods are widely studied for their stronger robustness and scalability.  The classical sampling path planning algorithm-based multi-robot motion planner is one of the mainstream research trends. In \cite{desaraju2011decentralized},   Desaraju et al. present DMA-RRT. This planner combines the rapid-exploring random tree (RRT) with the distributed multi-robot cooperative planning paradigm and can generate multiple paths for different robots.  However, DMA-RRT requires that each robot can communicate with others, and DMA-RRT does not consider the dynamic constraints of each robot.  In \cite{le2019multi}, Duong Le et al. utilize multi-agent search methods to guide the sampling-based planning process of each robot to effectively multi-robot motion planning problems with kino-dynamic constraints. The common problem of the above methods is that they do not separate the global and local planners. The whole planning process still relies on the priori map.

The velocity obstacle (VO) based methods are another type of decentralized multi-robot motion planning method that has been widely researched.  VO-based methods are based on reactive mechanisms and have the advantages of real-time and high efficiency.  VO does not consider the interaction pattern of other active agents in the same space \cite{douthwaite2018comparative}, which can cause the oscillation problem of final trajectories.  Reciprocal velocity obstacle (RVO) and its variants introduce implicit speed selection mechanisms and adopt probabilistic approaches to deal with uncertain interaction issues \cite{snape2011hybrid}. This operation can effectively solve the oscillation problem. Optimal reciprocal collision avoidance (ORCA) and its variants further improve the efficiency of real-time multi-robot trajectory generation.  In each time step, the robot obtains the velocity state of other robots and constructs its own speed space without collision.  The intersecting area from different speed spaces of robots constitutes the final optimal speed space.  Ultimately, the optimal joint motion strategy can be derived by solving linear programming problems \cite{van2011reciprocal}.  In \cite{douthwaite2018comparative},  scholars have demonstrated the superiority of ORCA over other VO-based methods in ideal scenarios.  ORCA can only handle motion planning problems where each robot is isomorphic and does not have non-holonomic constraints. In \cite{wang2018reciprocal,huang2019generalized},    researchers optimize the limitation of ORCA so that it can deal with more general task scenarios.  It is worth noting that there are many prerequisites for the application of ORCA and its variants. First, the robot has the perfect perception ability and can obtain the position, shape, and motion strategy information of $n$ robots within a specific range of itself.

In recent years, with the development of machine learning,  learning-based model-free multi-robot motion planning methods have gradually become a research hotspot.  For example, in \cite{qureshi2020motion}, AH Qureshi et al. present motion planning networks (MPNet). MPNet is a neural motion planner based on the continuous learning paradigm and the expert supervision.  MPNet successfully constructs the mapping relationship between raw sensor streaming data and final bidirectional connected paths.  In \cite{riviere2020glas}. Riviere Benjamin et al. propose global-to-local safe autonomy synthesis (GLAS)  for multi-robot motion planning. GLAS integrates centralized methods with centralized methods. In GLAS, each robot acquires a distributed optimal policy through imitation learning. The expert experience of imitation learning comes from a priori centralized global optimal motion planner with resolution completeness.  The limitation of the above methods  is that the performance of the algorithm is depends too much on the quality of the priori expert demonstrations or the labeled dataset.

Reinforcement learning (RL) approaches have received more attention than data-driven learning methods. RL is model-free and can learn in a ``trial-and-error" manner. Through interacting with the environments, the RL agent iteratively learns the policy to maximize the cumulative return \cite{he2021review,wang2021research}.  This pattern allows RL to integrate the learning process with the decision-making process and achieve end-to-end mapping from raw state input to the action output.   In the single mobile robot navigation task, RL-based motion planners have proven to deal with dense and dynamic scenarios without priori maps \cite{wang2018learning, shi2019end, he2021asynchronous}. Likewise, there are many scholars dedicated to researching multi-robot cooperative motion planning problems. Many studies have demonstrated that centralized training with decentralized execution (CTDE) -based multi-agent reinforcement learning (MARL) algorithms are pretty suitable for handling multi-robot collaborative planning problems \cite{chen2017decentralized, long2018towards}.  Compared to centralized MARL,  CTDE-based MARL has better scalability.  Besides, the centralized training process in CTDE avoids credit assignment and dynamic environment issues in decentralized MARL methods. In \cite{lowe2017multi, iqbal2019actor,yu2021surprising }, CTDE-based MARL algorithms such as MADDPG, MAAC, MAPPO have achieved great results in cooperative navigation scenarios in multi-agent particle simulation environments developed by OpenAI.  However, these works do not consider the trajectory shape and kino-dynamic constraints of real robots.    The ACL laboratory of MIT has made remarkable achievements in the field of  RL-based multi-robot motion planning \cite{chen2017decentralized, chen2017socially, semnani2020multi, everett2021collision}. In \cite{chen2017decentralized}, they propose collision avoidance with deep RL (CADRL) to solve multi-robot collaborative motion planning problems.  In \cite{chen2017socially}, they present socially aware CADRL (SA-CADRL) to address pedestrian-robot interaction issues on the basis of CADRL. Later, they consider the stochastic behavior model and introduce a supervised learning stage to solve the problem of dynamic environment information encoding \cite{everett2021collision}. In \cite{semnani2020multi} ,  they redesign the reward function and integrate original GA3C-CADRL with the force-based motion planning method to solve the long-range navigation problem.  The above researches promote the development of RL-based multi-robot motion planning. However, these works utilize the agent-level information of all robots. The observation input of each neural network contains policy information of other robots.

To cope with the above limitations, we develop a hybrid motion planner suitable for non-communication and partial observable multi-robot motion planning tasks.  The specific architecture is shown in Fig. \ref{Fig_TrajOptim}. We combine the advantage of the CTDE-paradigm based MARL algorithm with the minimal snap-based trajectory optimization algorithm to establish an end-to-end mapping from local observations to the final executable, dynamic-feasible, smooth, and continuous trajectories of multiple robots.  In the front-end waypoints searching module,  we utilize the pre-trained MARL model to generate collaborative discrete waypoints. We extend soft actor-critic (SAC) to a MARL algorithm to handle the decentralized partially observable Markov decision problem (Dec-POMDP). In the back-end trajectory optimization module, we construct an optimization problem for solving the optimal coefficients of the continuous trajectory polynomials.  To sum up, our contributions are summarized as follows:

\begin{itemize}
	\item We propose a hybrid multi-robot motion planner that separates the front-end waypoints searching module and the back-end trajectories optimization module. This overall motion planner is model-free and does not rely on map priors.
	\item We develop an off-policy multi-agent soft actor-critic-based cooperative waypoints searching method with self-adjusting policy entropy. This algorithm follows the CTDE paradigm. We set up multiple multi-robot motion planning scenarios under local observable and non-communication conditions and various baselines to verify the performance of this method.
	\item We utilize a minimal snap trajectory optimization method with safety zone constraints to do the post-processing. This module performs shape adjustment and optimization for discrete paths and outputs smooth, continuous, and executable joint trajectories. Also, this separated module reduces the difficulty of offline MARL training and avoids the problem of non-convergence caused by introducing too many optimization goals in the RL framework.
\end{itemize}

The remainder of this paper is organized as follows. The details of the front-end waypoints searching module and the back-end trajectory optimization module are described in Section II and Section III, respectively. Section IV presents the experimental results. Section V concludes this paper.

\section{Front-end Waypoints Searching}
In this section,  we describe in detail the front-end waypoints search module of the hybrid motion planner proposed in this paper. The specific content includes an introduction of the multi-agent soft actor-critic (MASAC) MARL framework based on the CTDE paradigm and a specification of the configuration of the state space, the action space, and the reward function in the multi-robot collaborative waypoints searching task.

\subsection{Multi-agent Soft Actor Critic Framework}
\subsubsection{Soft Actor Critic}
Before the soft actor-critic (SAC) algorithm is proposed, mainstream single agent model-free RL algorithms have some limitations. For example, the sample complexity of high-dimensional tasks leads to low sampling efficiency; a large number of hyper-parameters leads to unstable algorithm performance and weak generalization ability. The off-policy SAC balances sample utilization and algorithm stability. In addition,  the stochastic policy selection and the policy entropy mechanism are integrated into SAC. This operation enables SAC to encourage policy exploration by maximizing policy entropy, thus assigning nearly equal probability to those near-optimal actions with similar action-state values, avoiding repeatedly choosing the same action to fall into the sub-optimality. Meanwhile, the maximizing reward item ensures that the update process of the algorithm does not deviate from the overall optimization direction. Therefore, compared to DDPG, TD3, and other deterministic and continuous control RL algorithms, SAC has stronger policy exploration ability, generalization ability, and robustness, and thus is widely used in the field of robot learning.

The core of SAC described above can be represented as follows:
\begin{equation}\label{equation_sac}
\pi_{\max }^{*}=\arg \max _{\pi} \sum_{t=0}^{T} E_{\left(s_{t}, a_{t}\right) \sim \rho_{z}}\left[r\left(s_{t},a_{t}\right)+\bm{\alpha} H\left(\pi\left(\cdot \mid s_{t}\right)\right)\right]
\end{equation}
where ${H(\pi(\cdot \mid s_t))=-\log \pi(\cdot \mid s_t)}$ is the policy entropy, which represents the degree of randomization of the policy. $\bm{\alpha}$ is the temperature coefficient, which represents whether the optimization objective of SAC is more inclined to maximize rewards or maximize policy entropy. ${\sum_{t=0}^{T} \mathbb{E}_{\left(\mathbf{s}_{t}, \mathbf{a}_{t}\right) \sim \rho_{\pi}}\left[r\left(\mathbf{s}_{t}, \mathbf{a}_{t}\right)\right]}$ is the basic maximization reward item in the objective function.

\subsubsection{Multi-agent Soft Actor Critic Training Paradigm}
The policy entropy maximization property and the hyperparameter insensitivity property of SAC determine that the agent naturally has certain robustness and generalization. This advantage makes SAC more suitable for robot learning with external disturbances and uncertain factors. Therefore, we propose a MARL training paradigm for the multi-robot waypoints searching method based on SAC.

\begin{figure}[!t]
	\centering
	\includegraphics[width=3.5 in]{Fig_training.pdf}
	\caption{The detailed structure of MASAC off-policy training paradigm.}
	\label{Fig_TrajOptim}
\end{figure}

MASAC follows the CTDE paradigm. The robots represented by each actor are independent and cannot communicate with others. At each timestep, no robot can obtain the current motion policy of others. The advantage of CTDE is that there is no need to consider the credit assignment problem among agents. The optimization objective is to maximize the total reward, which is the sum of the rewards obtained by each robot after interacting with the environment. The reward function of each robot are independent of each other and only need to be set separately according to their own task requirements. The detailed structure of the MASAC training paradigm is illustrated in Fig. \ref{Fig_TrajOptim}. In this paradigm, the multi-robot waypoints search task can be modeled as a Dec-POMDP.
\begin{equation}\label{DecPOMDP}
	G=\langle \bm{\hat{S}}, \bm{A}, \bm{P}, \bm{R}, \bm{\Omega}, \bm{I}\rangle
\end{equation}
where ${i \in I=\{1,2, \ldots, N\}}$ represents the index set of each robot.  ${\bm{\hat{S}}=\langle \bm{S},\bm{O} \rangle}$ includes the global state and the collection of local observations of each robot. ${a_{i} \in A, \boldsymbol{a} \in A^{N}}$ represents the joint action of all robots at timestep ${t}$. ${\bm{R}=\{R_1,R_2,...R_N\}}$ is a tensor containing the reward signal of each agent.
${P\left(s^{\prime} \mid s, \boldsymbol{a}\right)}$ is the transition function from the current state to the next state. ${o_{i} \in \bm{\Omega} \sim \bm{O}(s, i)}$ is the observation function.

The global critic consists of two sets of action-value networks, one of which is a global twin soft Q network, and the other is the global target Q network that performs soft update based on the parameters of the twin soft Q network. Also, each actor includes a target policy network. The total twin soft Q Network outputs a pair of total Q values. We tend to take a smaller Q to weaken the overestimation bias in the Q learning process. The utilization of target networks of the global critic and actors is intended to stabilize the training process.

The loss function of the global twin soft Q network is shown as follows:
\begin{equation}\label{Eq_LossFunQ}
\begin{aligned}
	L(\theta_i) &=E_{\mathcal{D}}\left[\left(r_i\left(\boldsymbol{s}_{t}, \boldsymbol{o}_{t}, \boldsymbol{a}_{t}\right)+\gamma \min _{j \in 1,2} Q_{\theta_{j}}^{\prime t a r}\left(\boldsymbol{s}_{t+1}, \boldsymbol{o}_{t+1}, \tilde{\boldsymbol{a}}_{t+1}\right)\right.\right.\\
	&\left.\left.-Q_{\theta}\left(\boldsymbol{s}_{t}, \boldsymbol{o}_{t}, \boldsymbol{a}_{t}\right)\right)^{2}\right]
\end{aligned}
\end{equation}
where
\begin{equation}\label{Eq_Qtar_prime}
\begin{aligned}
	Q_{\theta_{j}^{\prime}}^{\prime t a r}\left(\boldsymbol{s}_{t+1}, \boldsymbol{o}_{t+1}, \boldsymbol{a}_{t+1}\right) &=\mathrm{E}_{\pi^{i}}\left[Q_{\theta_{j}^{\prime}}^{\operatorname{tar}}\left(\boldsymbol{s}_{t+1}, \boldsymbol{o}_{t+1}, \tilde{\boldsymbol{a}}_{t+1}\right)\right. \\
	&\left.-\alpha_{t}^{i} \log \pi^{i}\left(\tilde{\boldsymbol{a}}_{t+1} \mid o_{t+1}^{i}\right)\right].
\end{aligned}
\end{equation}
The joint trajectory ${\left(\bm{r}_t,\bm{s}_t, \bm{o}_t, \bm{a}_t, \bm{o}_{t+1}, \bm{s}_{t+1} \right)}$ is sampled from the replay buffer ${\mathcal{D}}$ at every specific timestep. ${\tilde{\bm{a}}_{t+1}}$ is obtained by re-sampling from current policy networks $\pi$. ${\gamma}$ is the discounted factor. ${\min _{j \in 1,2} Q_{\theta_{j}}^{\prime t a r}}$ is the twin soft Q values. We take the minimal one as the target. Equation (\ref{Eq_Qtar_prime}) illustrates its calculation process. In this equation, the policy entropy item is derived from each soft reward function
\begin{equation}
r_i^{soft}\left(\bm{s}_{t}, \bm{o}_t,\bm{a}_{t}\right)=r_i\left(\bm{s}_{t}, \bm{o}_t,\bm{a}_{t}\right)+\gamma {\alpha}_t^i \mathbb{E}_{\pi^{i}} H\left(\pi_{tar}^i\left(\cdot \right)\right).
\end{equation}
where ${Q_{\theta_{j}^{\prime}}^{\operatorname{tar}}}$ is the target Q network with parameters ${\theta^{\prime}}$. ${{\alpha}_t^i}$ is the exclusive temperature coefficient of each actor, which is used to control the weight of the policy entropy.

The loss function of each policy network is shown as follows:
\begin{equation}\label{Eq_policyUpdate}
	\begin{aligned}
	L(\phi_i)&=\mathbb{E}_{\bm{o}_{t} \sim \mathcal{D}, \varepsilon \sim \mathcal{N}}\left[\alpha_t^i \log \pi_{\phi_i}\left(f_{\phi_i}\left(\varepsilon_{t} ; \bm{o}_{t}\right) \mid \bm{o}_{t}\right)\right. \\
	&\left.-Q_{\theta}\left(\bm{s}_t, \bm{o}_{t}, f_{\phi_i}\left(\varepsilon_{t} ; \bm{o}_{t}\right)\right)\right].
	\end{aligned}
\end{equation}
This function is derived from the soft policy iteration procedure. We also utilize the reparameterization trick here. For each sample of ${\pi_{\phi_i}(\cdot \mid \bm{o})}$, it is jointly determined by the current observation $\bm{o}$, parameters ${\phi}$, and the independent noise ${\varepsilon}$ that conforms to the standard normal distribution. We adopt ${f_{\phi_i}\left(\varepsilon_{t} ; \bm{o}_{t}\right)}$ to represent this squashing function. Its specific form is shown as follows:
\begin{equation}\label{Eq_SquashingGaussian}
	{f_{\phi_i}}(\bm{o}, \varepsilon )=\tanh \left(\mu_{\phi_i}(\bm{o})+\sigma_{\phi_i}(\bm{o}) \odot \varepsilon \right), \quad \varepsilon \sim \mathcal{N}(0, I)
\end{equation}
where ${\text{tanh}}$ activation function limits the final action to ${[-1,1]}$. $\mu_{\phi_i}$ and $\sigma_{\phi_i}$ are the outputs of the policy network. This trick allows us to optimize the policy parameters by computing the gradient of soft Q value directly.

In addition, in order to better control the tradeoff between the exploitation and the exploration, and improve the parameters insensitivity property of the algorithm, we adopt an adaptive learning approach to control the temperature coefficient $\alpha_i$ of each actor. The purpose is to make each agent in MASAC pay more attention to exploring to improve the policy diversity in the early stage, and more inclined to utilize effective strategies to improve the training stability in the later stage. The objective function of $\alpha_i$ at timestep $t$ is as follows:
\begin{equation}
	L(\alpha_i)=\mathbb{E}_{\bm{a}_{t} \sim \pi^i_{t}}\left[-\alpha_i \log \pi^i_{t}\left(\bm{a}_{t} \mid \tau_{t}\right)-\alpha_i \overline{\mathcal{H}}\right]
\end{equation}
$\overline{\mathcal{H}}$ is the target entropy,  usually set to the dimension of the action space ${-\text{dim}(\mathcal{A}_i)}$ of each actor.

\subsection{RL Training Framework Configuration of Waypoints Searching Task}
The previous sub-section focuses on the description of the MARL framework based on the CTDE training paradigm. In this sub-section, we combine the functional requirements of the front-end waypoints search module in the hybrid motion planner to introduce our configuration methods for the action space, the observation space, and the reward function.

The information contained in the continuous observation space of each mobile robot includes measurement data from onboard sensors and preset prior data. We hope that each mobile robot cannot directly obtain the action information of other robots, but learns to make inferences and predictions from limited local observation, and make its own optimal or near-optimal motion decisions. The specific configuration is as follows:
\begin{equation}
	\bm{o}_i = \left[i,\bm{v}_i,\bm{p}_i,\bm{\tilde{p}}_{\text{goal}},\bm{\tilde{p}}_{\text{others}},r_{\text{safe}}\right]
\end{equation}
where $i$ is the index number of current robot. $\bm{v}_i$ is the speed strategy of the robot $i$ at current timestep. $\bm{p}_i$ is the current global coordinate of the robot $i$ at current timestep. $\bm{\tilde{p}}_{\text{goal}}$ represents the relative coordinate of the goal with respect to the robot $i$. $\bm{\tilde{p}}_{\text{others}}$ represents the relative coordinate of the position of other robots with respect to the robot $i$. $r_{\text{safe}}$ denotes the radius of the rectangular-shaped safety zone of the robot $i$. Moreover, the global state configuration to input global soft Q network is as follows:
\begin{equation}
	\bm{s}_i = \left[ \bm{P}_{\text{robots}}, \bm{P}_{\text{goals}} \right]
\end{equation}
It contains global positions and corresponding target points of all robots at a specific timestep.

The continuous action space of each mobile robots is set to $\bm{a}_i=\left[a_x, a_y\right] $. $a_x$ and $a_y$ are the components of the resultant acceleration generated by the robot actuator on the x-axis and y-axis, respectively. They are all in the value range of $[-1,1]$. In the actual deployment, the action range can change according to the complexity of the external environment.

In this paper, we assume that each mobile robot is homogeneous and has the same task objective. We design the reward function of each mobile robot as follows:
\begin{equation}\label{Eq_reward}
	r_t\left({\bm{o}_i,\bm{s}_i}\right)=\left\{\begin{array}{ll}
		1 & \text { if } d_g<0.1 \\
		-0.35 & \text { if } \text{any}(\bm{d}_{\text{robot}})<0 \\
		-0.35 + 0.05d_{robot}&\text{ if } 0<\text{any}(\bm{d}_{\text{robot}})<r_{th}\\
		d^{t-1}_g-d^{t}_g & \text { otherwise }.
	\end{array}\right.
\end{equation}
where ${d_{g}=\left\|\bm{p}_{g}-\bm{p}\right\|_{2}}$ is the Euclidean distance between the current position of the robot $i$ and the target position. ${\bm{d}_{robot}}$ is the distance vector between the robot $i$ and other mobile robots. $r_{th}$ is the safety distance threshold between two mobile robots. It can be found that in addition to the final goal-reaching reward and collision penalty, we add many intermediate state reward signals to enable the robot to learn continuously and avoid some problems caused by the reward sparsity problem.

Finally, after setting up MASAC training framework according to the above description and running the centralized training process with multiple episodes, the pre-trained actor model can be directly uploaded to the corresponding robot to perform decentralized collaborative waypoints online searching tasks without the global perfect sensing assumption.

The details of training the MASAC-based front-end waypoints searching model are summarized in Algorithm \ref{Algortihm_Front_end}.

\begin{algorithm}[htb]
	\caption{MASAC Training Paradigm.}
	\label{Algortihm_Front_end}
	\begin{algorithmic}[1]
		\State \textbf{Initialize} each policy network parameters ${[\phi_i]}_{i=1}^n$.  Initialize centralized twin soft Q networks   ${[\theta_i^1,\theta_i^2]}_{i=1}^n$. Initialize temperature coefficients of policy entropy ${[\alpha_i]}_{i=1}^n$.
		\State \textbf{Initialize} replay buffers $\mathcal{D}$, independent noise $\varepsilon \sim \mathcal{N}(0,1)$.
		\State \textbf{Initialize} parameters of target networks ${\phi}^{\prime}_i \leftarrow \phi_i, {\theta}^{\prime}_i \leftarrow \theta_i$.
		\For {$eps = 1$ to $M$}
			\If {$eps<M_{init}$}
				\State Warm-up stage.
			\EndIf
		\For{$t = 1$ to $eps_T$}
			\State select action $\bm{a}^t$ according to equation (\ref{Eq_SquashingGaussian}).
			\State $\bm{s}^{t+1} \sim p\left(\bm{s}^{t+1} \mid \bm{s}^{t}, \bm{a}^{t}\right)$.
			\State add $( \bm{s}^t, \bm{o}^t, \bm{a}^t, \bm{r}^t,
			\bm{o}^{t+1}, \bm{s}^{t+1} )$ to $\mathcal{D}$.
			\If {$\text{update}$}
				\State $\theta_i \leftarrow \theta_i - \alpha^{Q}_i \nabla_{\boldsymbol{\theta}_{i}} L(\theta_i) $
				\State $\phi_i \leftarrow \phi_i - \alpha^{\pi}_i \nabla_{\boldsymbol{\phi}_{i}} L(\phi_i) $
				\State $ \alpha_i \leftarrow \alpha_i - \nabla_{\boldsymbol{\alpha}_{i}} L(\alpha_i)$
			\EndIf
			\If {$\text{soft update}$}
				\State $\theta_i^{\prime,j} \leftarrow \tau \theta_i^{j} - (1-\tau)\theta_i^{\prime,j}$, $ j \in [1,2] $
				\State $\phi_i^{\prime} \leftarrow \tau \phi_i - (1-\tau)\phi_i^{\prime}$
			\EndIf
		\EndFor
			\State \textbf{end for}
		\EndFor
		\State \textbf{end for}
	\end{algorithmic}
\end{algorithm}

\section{Back-end Trajectory Optimizing}

At present, in mainstream RL-based motion planning methods, mobile robots would directly track unsmooth folded segment trajectories formed by the discrete waypoints. However, in practical applications, the speed command and the acceleration of the mobile robot cannot change abruptly. This limitation makes mobile robots unable to follow the preset trajectory precisely, thus producing additional motion overshoot and resulting in safety risks such as collisions. Meanwhile, mobile robots need to accelerate and decelerate frequently  to track such polyline trajectories, which is extremely energy-consuming. If the trajectory generation and optimization are directly integrated into the RL architecture to realize the end-to-end mapping from state input to executable smooth trajectories, iterative reward function debugging processes have to be considered. Also, this method requires coupling kino-dynamic constraints and safety constraints of the robots, which undoubtedly increase the convergence difficulty of the algorithm.

In this section, we separate the trajectory optimization process from the RL-based end-to-end motion planner, and utilize quadratic programming methods to generate optimal trajectories with minimal snap. This cascaded hybrid motion planning approach facilitates us to introduce kino-dynamics and safety constraints to solve smooth and executable trajectories with continuous velocity, acceleration, and jerk. Also, the operation of separating the trajectory optimization process avoids the introduction of a multi-objective reward form in the training phase, reduces the convergence difficulty, and allows the robot to focus on learning to obtain collision-free discrete waypoints.  Most importantly, from the front-end waypoints search to the back-end trajectory optimization, this hybrid motion planner helps us avoid the complex modeling process.

\subsection{Trajectory Generation}
A piece of continuous trajectory formed by any two waypoints can be described by a segment of $n$-th order polynomial function:
\begin{equation}\label{Eq_traj}
	f_i(t)=\left[1, t, t^{2}, \ldots, t^{n}\right] \cdot \bm{p}_i
\end{equation}
where
\begin{equation}
	\bm{p}_i=\left[p_{0,i}, p_{1,i}, \ldots, p_{n,i}\right]^{T}
\end{equation}
${[p_{0,i},p_{1,i},...,p_{n,i}]}$ are the trajectory parameters, with the number of $n+1$. So, we can derive the speed, acceleration, jerk,  and snap representation of this two-point trajectory according to equation (\ref{Eq_traj}). The details are shown in equation (\ref{Eq_vajs}).
\begin{equation}\label{Eq_vajs}
	\begin{array}{l}
v_i(t)=f^{\prime}(t)=\left[0,1,2 t, 3 t^{2}, 4 t^{3}, \ldots, n t^{n-1}\right] \cdot \bm{p}_i \\
a_i(t)=f^{\prime \prime}(t)=\left[0,0,2,6 t, 12 t^{2}, \ldots, n(n-1) t^{n-2}\right] \cdot \bm{p}_i \\
\text{jerk}_i(t)=f^{(3)}(t)=\left[0,0,0,6,24 t, \ldots, \frac{n !}{(n-3 !)} t^{n-3}\right] \cdot \bm{p}_i \\
\text{snap}_i(t)=f^{(4)}(t)=\left[0,0,0,0,24, \ldots, \frac{n !}{(n-4 !)} t^{n-4}\right] \cdot \bm{p}_i
\end{array}
\end{equation}
Now, we assume there exists $M$ mobile robots. For any robot $m$, there are $K+1$ waypoints generated by the MARL-based front-end. Therefore, the whole motion trajectory of this robot composed of $K$-segment polynomials for this robot has the following representation:
\begin{equation}
	\begin{aligned}
	f(t)=\left\{\begin{array}{ll}
f_{1}(t) \doteq \sum_{i=0}^{n} \bm{p}_{1, i} t^{i}   &T_{0} \leq t \leq T_{1} \\
f_{2}(t)  \doteq \sum_{i=0}^{n} \bm{p}_{2, i} t^{i}  &T_{1} \leq t \leq T_{2} \\
\vdots  \\
f_{K}(t)  \doteq \sum_{i=0}^{n} \bm{p}_{K, i} t^{i}  &T_{K-1} \leq t \leq T_{K}
\end{array}\right.
	\end{aligned}
\end{equation}
where $T$ is the time node of each segment.

\subsection{Minimal Snap with Safety Zone Constraints}
Combined with the above derivation, our optimization objective is to select a set of optimal polynomial parameter combinations $P=[P_{R1},P_{R2},...,P_{RM}]$ (where $P_{Rm}=[\bm{p}_1^m,\bm{p}_2^m,...,\bm{p}_{K}^m]$) under various constraints to minimize the snap of the trajectory of each robot. In this way, the actuator thrust of each robot changes as smoothly as possible, thereby minimizing energy consumption.

For a single trajectory of the mobile robot $m$, the cost function of minimal snap can be written as:
\begin{equation}
\begin{aligned}
J(T) &=\min \int_{0}^{T}\left(f^{(4)}(t)\right)^{2} \mathrm{d} t=\min \sum_{i=1}^{K} \int_{T_{i-1}}^{T_{i}}\left(f^{(4)}(t)\right)^{2} \mathrm{d} t \\
&=\min \sum_{i=1}^{K} \bm{p}_i^{T} \bm{Q}_{i} \bm{p}_i
\end{aligned}
\end{equation}
where $f^{4}(t)$ can be obtained in equation (\ref{Eq_vajs}). $\bm{p}_i$ represents the polynomial parameters of each segment. $\bm{Q}_i$ is the Hessian matrix, the detail is as follows:
\begin{equation}
\bm{Q}_i=\left[\begin{array}{cc}
0_{4 \times 4} & 0_{4 \times(n-3)} \\
0_{(n-3) \times 4} & \begin{array}{cc}\frac{r(r-1)(r-2)(r-3)c(c-1)(c-2)(c-3)}{r+c-7}\times\\ \left(t_{i}^{(r+c-7)}-t_{i-1}^{(r+c-7)}\right)
\end{array}\end{array}\right]_{r\times c}
\end{equation}
where $r$ and $c$ respectively represent the number of rows and columns of $Q_i$.

In the task scenario described in this paper, the constraints of the trajectory optimization process of each mobile robot include several equality constraints and several inequality constraints.
\subsubsection{Equality Constraints}
First, we introduce the initial and terminal state constraints of the mobile robot, including  position, speed, and acceleration:
\begin{equation}
f^{(d)}\left(T_{0,T}\right)=s_{0,T}^{(d)}.
\end{equation}
We transform it into the standard input form of the QP optimizer:
\begin{equation}
	\bm{A}_{0,T}\bm{p}_{0,T} = \bm{s}_{0,T}
\end{equation}
where $\bm{A}_{0,T} = [\bm{A}_0,\bm{A}_T]^T$ is the mapping matrix between polynomial parameters and state, $\bm{s}=[x,v,a]^T$.

The other type is the continuity constraint. We ensure the continuity of the whole trajectory by constraining endpoint derivatives of the segment $i$ to be equal to initial derivatives of the segment $i+1$:
\begin{equation}
\left.\begin{array}{l}
f_{i}^{(d)}\left(T_{j}\right)=f_{i+1}^{(d)}\left(T_{j}\right) \\
\Rightarrow\left[\bm{A}^{i}-\bm{A}^{i+1}\right]\left[\begin{array}{c}
\bm{p}_{i} \\
\bm{p}_{i+1}
\end{array}\right]=0
\end{array}\right.
\end{equation}

\subsubsection{Inequality Constraints}
First, we need to set ${v_{min}, v_{max}, a_{min}, a_{max}}$ to constrain the motion speed and acceleration of each mobile robot. Moreover, we introduce a set of rectangular-shaped safety zones to constrain the middle points of each segment in the trajectory instead of fixing each middle point according to the classical minimal snap method. In detail, we perform multiple middle path point sampling processes in every intermediate segment. For each path point, we add two range inequality constraints on the x-axis and y-axis:
\begin{equation}
	\begin{array}{c}
{A_{T_j} \bm{p}_i \leq f_i\left(T_{j}\right)+d_{\text{safe}}} \\
{-A_{T_j} \bm{p}_i \leq-\left(f_i\left({T_j}\right)-d_{\text{safe}}\right)}
\end{array}
\end{equation}
This corridor-based soft constraints approach contains an implicit time allocation mechanism. Furthermore, it avoids the problem of overshoot when optimizing trajectory by the classical minimal snap method with solid constraints.

\begin{table}[!t]
	\centering
	\caption{Hyper-parameters configuration of MASAC training process}
	\setlength{\tabcolsep}{5mm}
	\label{Table_MASACTraining}
	\begin{tabular}{lc}
		\toprule[1pt]
		\textbf{Parameters}   & \textbf{Value}                     \\ \hline
		$\bm{\alpha}$(init)        & [0.01,0.01,0.01]                               \\
		target entropy        & [-dim(12),-dim(12),-dim(12)]                          \\
		learning rate(actor)  & 0.001                               \\
		learning rate(critic) & 0.001                              \\
		optimizer             & Adam                               \\
		tau                   & 0.005                              \\
		mini-batch size       & 2048                               \\
		episodes              & 5e4                                \\
		warm-up episodes      & 1e3                                \\
		memory length         & 1e6                                \\
		timesteps per episode & 50								  \\
		update frequency      & 50                                \\
		actor delay frequency & 1                                  \\
		$d_{\text{safe}}$(m)         & 0.3 							   \\
		\bottomrule[1pt]
	\end{tabular}
\end{table}

\begin{figure}[!t]
	\centering
	\includegraphics[width=3.0 in]{RewardCurves.pdf}
	\caption{The average reward curves for MASAC, MATD3, MADDPG, and MAAC on the 3-robot waypoints searching task scenario.}
	\label{Fig_RewardCurves}
\end{figure}

\begin{figure*}[!t]
	\centering
	\includegraphics[width=5.0 in]{AlphaDiscussion.png}
	\caption{(a) represents the self-learning process of average temperature coefficients of different agents with different initial values. (b) represents the average reward curves of MASAC methods with different initial $\alpha$ during the training process. }
	\label{Alpha_tuning}
\end{figure*}

\begin{figure*}[]
	\centering
	\includegraphics[width=6 in]{3robots_planning_results.png}
	\caption{The visualization inference results of different algorithms in the unknown 3-robot fully interaction planning scenario. The rightmost bar chart compares the minimal distance between robots during the planning process.}
	\label{Fig_Comparision3agents}
\end{figure*}

\begin{table*}[!t]
	\centering
	\caption{Performance metrics comparison of 3 robots collaborative path planning based on different algorithms}
	\label{Tab_all_metrics}
	\begin{tabular}{cccccc}
		\toprule[1pt]
		\textbf{Algorithms}                         & \multicolumn{5}{c}{\textbf{Different Performance Indicators}}                                                                         \\ \hline
		\multicolumn{1}{l}{\textbf{Motion Planner}} & \textbf{Avg Total Distance} & \textbf{Avg Total Time} & \textbf{Avg Total Reward} & \textbf{Collsion Numbers} &\textbf{Success Rate} \\ \hline
		ORCA                                        & 1.928                       & 114.538                 & -                         & 10(9049)                  & 73.627\%              \\
		MATD3                                      & 3.497                       & 124.013                 & -28.804                   & 25(9673)                  & 92.301\%              \\
		MAAC(CAS)                                   & 3.209                       & 123.679                 & -33.695                   & 14(9647)                  & 92.051\%              \\ \hline
		\textbf{MASAC-auto}                         & \textbf{2.969}              & \textbf{118.897}        & \textbf{-39.545}          & \textbf{2(9274)}          & \textbf{93.592\%}     \\ \bottomrule[1pt]
	\end{tabular}
\end{table*}

\begin{figure*}[!t]
	\centering
	\includegraphics[width=5.5 in]{6_12agents_MODIFY.png}
	\caption{6-robot and 12-robot planning scenarios. ORCA is still selected as an ideal comparison algorithm. The total planning distance length and average time consumption metrics of these two methods are compared in the rightmost bar charts.}
	\label{Fig_6agentsAnd12agents}
\end{figure*}

\section{Experiments}
In this section, we first present the training process and the collaborative planning effect of the MASAC-based front-end waypoints searching module, and analyze the comparison results between this algorithm with multiple model-free algorithms with the CTDE paradigm and a mainstream model-based algorithm. Then, we combine the pre-trained front-end waypoints module with the back-end trajectory optimization module, and demonstrate the trajectory planning effect of the final hybrid motion planner and make some analysis.

\subsection{Experiments of the Front-end Waypoints Searching Module}
We select the multi-robot collaborative navigation task with preset target points as the benchmark test scenario. As for the selection of the comparison baselines, we adopt state-of-the-art MATD3 (an improved algorithm on the basis of MADDPG) and MAAC algorithms which are all training under the CTDE paradigm. Furthermore, we set optimal reciprocal collision avoidance (ORCA) as the ideal baseline. ORCA is a commonly used multi-robot interaction algorithm based on the velocity obstacle (VO). ORCA is an ideal algorithm. It has several preconditions. First, all robots should be homogeneous. Besides, there exists a perfects communication assumption between robots, i.e., each robot can obtain the motion strategy of other robots at each timestep. This baseline facilitates us to compare the performance difference between our model-free and local observable method with the communicable method.

The centralized critic of the MASAC architecture of our waypoints searching module consists of two soft Q networks. Each soft Q network contains three fully connected hidden layers, and the number of units in each layer is $[1024,512,300]$. Each distributed policy network contains two fully connected hidden layers, and the number of units in each layer is $[500,128]$. During the training process, we introduce the warm-up training stage and reward scaling trick for stability purposes. Also, we assign an independent temperature coefficient to each actor and utilize the self-tuning approach to balance the exploration and exploitation. At the beginning of each training phase, we randomly initialize the positions of every robot and corresponding target point. The specific hyperparameters are configured as shown in Table \ref{Table_MASACTraining} (3-robot waypoints searching task).

We deploy all algorithms on a computer with AMD Ryzen7 5800H CPU and NVIDIA RTX 3060 GPU for training. We take the 3-robot waypoints search task as an example and list the average reward curves of different algorithms. The details are visualized in Fig. \ref{Fig_RewardCurves}. It can be found that MASAC with multiple self-learning temperature coefficients $\bm{\alpha}=[\alpha_1,\alpha_2,...,\alpha_n]$ can speed up the training process and help the robot obtain the maximum average reward without communication. Also, during the hyper-parameter-tuning process, we find that the initial value of the exploration temperature coefficient $\alpha$ has a certain impact on the final performance of MASAC. As shown in Fig. \ref{Alpha_tuning}, the subplot (a) shows that MASAC can indeed effectively adjust the exploration strength during the training process. Also, different initial $\alpha$ eventually converges to similar final states. Subplot (b) illustrates that when the initial value of each $\alpha_i$ is 0.01, MASAC has a faster convergence speed compared to the other three groups. However, it should not be ignored that a slight performance degradation occurs at the end of the training phase. A comprehensive comparison reveals that the degree of this decay disappears with increasing of the initial $\alpha_i$ of each agent.

The real-time inference test results are summarized in Fig. \ref{Fig_Comparision3agents}. We design an unknown scenario that enables each mobile robot could fully interact with others as an inference test benchmark. The results show that ORCA allows robots to interact with an approximate minimum safety threshold (0.301) under the condition of perfect communication. On the other hand, MATD3, MAAC, and MASAC are running under non-communication and partial observation conditions. The minimum distance between robots during the interaction process of MATD3 is 0.269. It can also be verified in the second sub-figure that the red robot collides with the green one at the 11-th timestep. The robots trained in MAAC architecture are relatively conservative. They tend to select suboptimal motion strategies to ensure a sufficient safety distance (0.564).

The other two MASAC-based waypoints searching methods obtain more intuitive collaborative planning results. Under the premise that other robot motion strategies cannot be obtained, each robot chooses to slow down or bypass in advance to avoid others. When we set the initial exploration temperature coefficient $\alpha_i = 0.01$, the robot behaves more conservatively (0.311). In contrast, the final result is closer to the ORCA when we adjust $\alpha_i$ to $0.1$ (0.330).  This is because a bigger exploration degree in the early stage of the training phase will increase the complexity and diversity of samples, which will sacrifice part of the training speed. However, the algorithm will therefore have better robust performance and inference generalization ability.

Further, we design an algorithm to randomly generate 100 task scenarios of 3-robot waypoints searching for evaluation. Also, we selected a variety of evaluation metrics, including the average distance sum of robots, the average time consumption of robots,  the average reward sum,  total collision numbers, and the search success rate within limited timesteps (Note that ``success" is recorded when all robots reach their target points.) All results are aggregated in Table \ref{Tab_all_metrics}. This result suggests that MASAC has better comprehensive waypoints searching performance and is closer to the performance of the ideal ORCA algorithm compared to MAAC(continuous action space version) and MATD3. Moreover, the MASAC-based searching method has the highest success rate with limited timesteps.

For further evaluating the inference ability of the proposed method, we set up more complex unknown fully interaction scenarios. The details are shown in Fig. \ref{Fig_6agentsAnd12agents}. The first waypoints searching scenario contains six robots. All robots interact fully from the same side to the other side. The second one is a 12-robot bilateral bi-directional interaction waypoints searching scenario.  It should be noted that in the 12-robot environments, we choose a smaller scenario scale to increase the difficulty of inference. We also utilize optimal ORCA(Global distance perception) as a comparison benchmark. The aggregated results in this figure suggest that the MASAC-based method has great scene generalization ability. In the 6-robot long-distance planning scenario, the proposed method approaches the performance of perfect sensing ORCA in the total distance length indicator, and has less timestep consumption. In the narrower 12-robot scenario, ORCA with the real-time policy sharing mechanism makes each robot go around a distance in advance. In contrast, our method achieves better results on the matrics of total distance length and average timestep consumption.

\subsection{Experiments of the Hybrid Motion Planner}

\begin{table}[!t]\label{Table_HybridMotionPlanner}
	\centering
	\setlength{\tabcolsep}{5mm}
	\caption{Hyper-parameters configuration of minimal snap trajectory optimization with safety zone constraints}
	\label{Table_TrajOptimProcess}
	\begin{tabular}{lc}
		\toprule[1pt]
		\textbf{Parameters}            & \multicolumn{1}{l}{\textbf{Value}} \\ \hline
		order $n$                         & 5                                  \\
		length of the safety zone $l$     & 0.1                                \\
		width of the safety zone $w$      & 0.1                                \\
		intermediate sampling interval $r$ & 0.05                               \\
		max iteration  $N$               & 1000                               \\
		action range $[a_{min},a_{max}]$   &
		[-1.0,1.0]                         \\
		\bottomrule[1pt]
	\end{tabular}
\end{table}

Based on the front-end MASAC waypoints searching module, we integrate the minimal snap trajectory optimization method with safety zone constraints. The final hybrid motion planner can generate dynamic-feasible, collision-free, and energy-optimal trajectories for multiple mobile robots navigating cooperative motion planning under no-communication and partial observation conditions. Besides, the velocity and the acceleration profiles of generated trajectories are smooth and continuous. This feature is conducive to the design of the low-level tracking controllers.  The hyperparameters configuration of the back-end trajectory optimizing module is summarized in Table \ref{Table_Backend}.

We take the 3-robot fully interactive scenario as an example. We preset the initial speed, initial acceleration, terminal speed, and terminal acceleration to the zero states. Meanwhile, we select the classical minimal snap method as the comparison benchmark. We aggregate the final trajectory generation results into the following Fig. \ref{Fig_FinalTraj_Comparsion}.
\begin{figure}[!t]
	\centering
	\includegraphics[width=3.5 in]{FinalTraj_Comp.png}
	\caption{The final trajectory of cooperative motion planning process of multi-robots. The minimal snap-based trajectory optimizer at the back-end can produce a smoother robot trajectory and a more reasonable speed arrangement scheme.}
	\label{Fig_FinalTraj_Comparsion}
\end{figure}

\begin{figure}[!t]
	\centering
	\includegraphics[width=3.4 in]{AfterTrajOptim_paj.png}
	\caption{The visualization of multiple kino-dynamic profiles of one of the robots in the task scenario. The results suggest that after optimizing by our back-end module, the acceleration and jerk curves of the trajectory are smoother, more continuous,  and have better executability.}
	\label{Fig_AfterTrajOptim}
\end{figure}

It can be found that the final trajectory generated by the hybrid motion planner with safety zone constraints is better than that of the classical minimal snap method. Attributing to the implicit time allocation mechanism, the back-end trajectory optimization module helps the hybrid motion planner to generate a more reasonable speed arrangement scheme. This facilitates the correction of anomalous segments of the robot trajectory under the no-communication and local observation conditions.

Fig. \ref{Fig_AfterTrajOptim} presents more precisely the changes of the position profile,  the acceleration profile, and the jerk profile before and after introducing our back-end trajectory optimizer. The dashed lines represent the fold trajectory of the output of the front-end waypoints searching module, and the solid line is the final executable trajectory. It can be found that the position profile is fine-tuned by our back-end trajectory optimizer. Moreover, this hybrid planner can output smoother, continuous, and differentiable acceleration action curves than the original discontinuous and mutable acceleration curves. This performance ensures the stability and smoothness of robots during autonomous motion. Also, within the input safety range, the more stable acceleration change process helps to protect the actuator of the robot and thus save energy consumption.

\begin{table}[!t]
	\centering
	\caption{Performance comparison of different back-end trajectory generation methods. }
	\label{Table_Backend}
	\begin{tabular}{ccccc}
		\toprule[1pt]
		\multirow{2}{*}{\textbf{Methods}} & \multicolumn{4}{c}{\textbf{Differenct Performance Indicators}}             \\ \cline{2-5}
		& $\textbf{cost}_{\text{comfort}}$&$\textbf{cost}_{\text{smooth}}$&$\textbf{cost}_{\text{energy}}$&$\textbf{L}_{\text{distance}}$ \\ \hline
		Original Wpts                     & 1.053e-2          & -28.955           & 5.489            & 2.081            \\
		Cubic Spline                      & 1.886e-5          & -51.558           & 14.298           & 2.945            \\
		Minimal Snap                      & 1.501e-5          & \textbf{-109.914}          & 0.278            & 2.945            \\ \hline		\textbf{Ours}                              & \textbf{1.201e-6}          & -107.171          & \textbf{0.018}            & \textbf{2.824}            \\
		\bottomrule[1pt]
	\end{tabular}
\end{table}

Furthermore, we introduce four performance metrics to quantify the final performance of our method and conduct multiple experiments. The average results are summarized in Table \ref{Table_Backend}. The indicator for evaluating the  straightness of the final trajectory is as follows:
\begin{equation}
	\textbf{cost}_{\text{comfort}} = \sum_{i=1}^{n-1}(x_{i-1}+x_{i+1}-2x_i)^2+(y_{i-1}+y_{i+1}-2y_i)^2
\end{equation}
where $[x_{i-1},x_{i},x_{i+1}]$ and $[y_{i-1},y_i,y_{i+1}]$ are the horizontal and vertical coordinates of any three adjacent points $[P_{i-1},P_{i},P_{i+1}]$ in the final trajectory. $n$ represents the number of trajectory points after discretization. Also, the indicator for evaluating the smoothness of the final trajectory can be represented as follows:
\begin{equation}
	\begin{aligned}
	&\textbf{cost}_{\text{smoothness}}=  \\
	&-\sum_{i=1}^{n-1} \frac{\left(x_{i}-x_{i-1}\right) \left(x_{i+1}-x_{i}\right)+\left(y_{i}-y_{i-1}\right) \left(y_{i+1}-y_{i}\right)}{\|P_iP_{i-1}\|_{2}{\|P_{i+1}P_{i}\|_{2}}}
	\end{aligned}
\end{equation}
$\textbf{cost}_\text{energy}=\sum_{i=1}^{n-1}(\bigtriangleup \bm{a})^2$ represents the force change during the robot movement, which can be used to measure energy consumption level. $\textbf{L}_\text{distance}$ represents the distance length.  Meanwhile, we select the ``front-end $+$ cubic spline" method and the ``front-end $+$  minimal snap" method as the benchmarks for comparison. The final results show that the output trajectories of our hybrid motion planner have better performance among the several groups given in this paper.

\section{Conclusion}
In this paper, we propose a model-free multi-robot hybrid motion planner based on the MASAC-based waypoints searching method and the minimal snap with safety zone constraints trajectory optimizer.  This planer can output smooth, continuous, and dynamic feasible cooperative trajectories under non-communication and local observable conditions. In the front-end of the planner, we utilize MASAC with autotune exploration temperature coefficients to train multiple robots offline to learn to search available joint discrete waypoints.  In the back-end of the planner,  we construct a minimal snap optimization objective and introduce dynamic and safety constraints to revise and improve known discrete waypoints. By solving a  quadratic programming problem, we can obtain the final collision-free multi-robot executable trajectories. We set multi-group multi-robot motion planning experiment scenarios and select several mainstream baselines and an ideal algorithm under perfect perception assumption as the comparison benchmarks. The final simulation results verify the superior performance of our method.  Among multiple performance metrics, our method is closer to the ideal algorithm among several baselines.  Moreover, the final results also show that the back-end optimizer can successfully improve the quality of the final cooperative planning trajectories.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{Introduction}
Dynamical mass generation (DMG) \cite{Brauner:2005hw,Maris:1999nt,Fischer:2003rp,Aguilar:2005sb,Bowman:2005vx,Aguilar:2010cn,Cloet:2013jya,Mitter:2014wpa,Binosi:2016wcx,Ayala:2006sv,Libanov:2005vu,Benes:2008ir} is, by definition, a non-perturbative phenomenon. Historically, it has mostly been relevant to QCD physics \cite{Cornwall:1981zr}, particularly to the lightest family of quarks due to limited capacity of QCD interactions to generate large dynamical masses. In the new physics, the phenomenon is equally relevant to masses considerably lighter than the electroweak scale. As existence of scalar sector in nature has become a possibility after the experimental finding of the Higgs boson \cite{pdg,Carena:2002es,Aad:2012tfa,Chatrchyan:2012xdj,MalbertionbehalfoftheCMSCollaboration:2018eqs}, ascertaining the extent of dynamical masses of new scalar fields due to different interactions \cite{Bezrukov:2021mio,Bezrukov:2013fka,Mukhanov:2005sc,Lee:2017qve,Bento:2000ah,Bertolami:2016ywc,Munoz:2017ezd,Martin:1997ns,Lee:2017qve} can not be underestimated.
\par
Yukawa interaction vertex is among the phenomenologically interesting vertices \cite{Schwartz:2013pla,Kaku:1993ym,Peskin:1995ev}. A model containing the vertex and the standard model (SM) Higgs field serves as an important avenue to explore the scalar sector. Two Higgs doublet models (2HDM) \cite{Gunion:2002zf,Cabrera:2020lmg,Altmannshofer:2020shb} are among widely studied models for various reasons. The model has its own strength in scalar sector, i.e. it presents two important scenarios once the SM Higgs \cite{Aad:2012tfa,Chatrchyan:2012xdj} is introduced in the theory. First, the model presents a mock up scenario of fermions interacting with the SM Higgs. Second, the model presents an opportunity to investigate how same particles from two different families interact with each other. A highly interesting aspect in both scenarios is that the physics of dark matter in scalar sector arising from possible Higgs-ultralight scalar ($m_{s}=O(10^{-22})$ eV) interactions \cite{Rindler-Daller:2013zxa,Harko:2014vya,Chavanis:2011zi,Huang:2013oqa,Marsh:2015xka,Hui:2016ltb,Primack:2009jr,Hu:2000ke} can also be studied in the model. It immediately justifies the investigation related to dynamical mass generation in the model.
\par
This paper is an extension of Wick Cutkosky (WC) model \cite{Darewych:1998mb,Sauli:2002qa,Efimov:2003hs,Nugaev:2016uqd,Darewych:2009wk} to incorporate two complex doublet fields, termed as the SM Higgs and the second Higgs fields for convenience, which interact with each other only by a real singlet scalar field. The paper is a continuation of the study of a two Higgs doublet model \cite{Mufti:2021lit}. One of the motives for the study is to understand extent of generation of dynamical masses for a scalar singlet field in different regions of the parameter space of the model. The scalar singlet acts as a mediator field between the two different, and otherwise non-interacting, complex doublet fields under the Yukawa interactions. The renormalized \footnote{The terms renormalized mass and physical mass are interchangeably used throughout the paper for the SM Higgs. The term scalar field is reserved only for the singlet scalar.} masses of the two complex doublet scalars are kept at their physical masses in order to study the phenomenon of DMG from the perspective of phenomenology. The physical mass of the SM Higgs is kept at its experimentally known value \cite{pdg} while the renormalized mass of the second Higgs is chosen for different scenarios mentioned above. A certain advantage of fixing masses is reduction of the parameter space to be explored. Hence, It effectively renders the model only a 3 dimensional parameter space in the Lagrangian. It greatly facilitates in training an algorithm \cite{Alpaydin14}, which is not a novel approach in quantum field theory \cite{Bachtis:2021xoh,Halverson:2020trp,Akutagawa:2020yeo}, over the sample of calculated dynamical masses in different regions of the parameter space, the details are given in the next section.
\par
The approach of Dyson-Schwinger equations (DSEs) \cite{Schwinger:1951ex,Schwinger:1951hq,Swanson:2010pw,Roberts:1994dr,Rivers:1987hi} is used for the study. The DSEs for the three field propagators are considered while the interaction vertices are fixed at their tree level form up to certain renormalization dependent terms, the details are given in the next section.
\par
There exist other renormalizable vertices to further extend the model \cite{Schwartz:2013pla,Kaku:1993ym,Peskin:1995ev}. However, such extensions may not be suitable for investigations using the approach of DSEs due to possible effects of further truncations and ansatz \cite{Roberts:1994dr} in the model. An extended model is to be reported somewhere else \cite{Mufti:2022abc} for which another non-perturbative method \cite{Ruthe:2008rut} is employed.
\section{Technical Details}
The Euclidean version of the Lagrangian \footnote{The technical details, a significant part of which can also be found somewhere else, are generously included to keep the section self-contained.} with the counter terms is given by
\begin{equation} \label{Lagrangian:eq}
\begin{split}
L = \frac{1}{2}(1+A) \partial_{\mu} \phi(x) \partial^{\mu} \phi(x) + \frac{1}{2} (m_{s}^{2}+B) \phi^{2}(x) + (1+\alpha) \partial_{\mu} h^{\dagger}(x) \partial^{\mu} h(x) \\ + (m_{h}^{2}+\beta) h^{\dagger}(x) h(x) + (1+a) \partial_{\mu} H^{\dagger}(x) \partial^{\mu} H(x) + (m_{H}^{2}+b) H^{\dagger}(x) H(x) \\ + (\lambda_{1}+C_{1}) \phi(x) h^{\dagger}(x) h(x) + (\lambda_{2}+C_{2}) \phi(x) H^{\dagger}(x) H(x)
\end{split}
\end{equation}
where $A$, $B$, $\alpha$, $\beta$, $a$, $b$, $C_{1}$, and $C_{2}$ are coefficients due to the counter terms in the Lagrangian \footnote{The bare scalar singlet mass ($m_{s}$) is kept in the equation \ref{Lagrangian:eq} for the sake of clarity.}. The real singlet scalar field is represented by $\phi(x)$, $h(x)$ is designated for the SM Higgs boson while $H(x)$ represents the second Higgs boson. The resulting DSEs for the field propagators are given below:
\begin{equation} \label{hdse:eq}
\begin{split}
D_{h}^{-1}(p)=(1+\alpha) p^{2} + m^{2}_{h} (1+\alpha) + 2 (1+A) (1+\alpha) (1+a) \sigma_{h} + \\ (\lambda_{1}+C_{1}) \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \Gamma_{1}(-p,q) D_{h}(q-p)
\end{split}
\end{equation}
\begin{equation} \label{Hdse:eq}
\begin{split}
D_{H}^{-1}(p)=(1+a) p^{2} + m^{2}_{H} (1+a) + 2 (1+A) (1+\alpha) (1+a) \sigma_{H} + \\ (\lambda_{2}+C_{2}) \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \Gamma_{2}(-p,q) D_{H}(q-p)
\end{split}
\end{equation}
\begin{equation} \label{sdse:eq}
\begin{split}
D_{s}^{-1}(p)=(1+A) p^{2} + m^{2}_{s} (1+A) + 2 (1+A) (1+\alpha) (1+a) \sigma_{s} + \\ (\lambda_{1}+C_{1}) \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{h}(q) \Gamma_{1}(q,-p) D_{h}(q-p)+ \\ (\lambda_{2}+C_{2}) \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{H}(q) \Gamma_{2}(q,-p) D_{H}(q-p)
\end{split}
\end{equation}
where the following definitions are used:
\begin{subequations} \label{mterms:eq}
\begin{align}
\beta = \alpha m^{2}_{h} + 2(1+A) (1+\alpha) (1+a) \sigma_{h}  \\
b = a m^{2}_{H} + 2(1+A) (1+\alpha) (1+a) \sigma_{H}  \\
B = A m^{2}_{s} + 2(1+A) (1+\alpha) (1+a) \sigma_{s}
\end{align}
\end{subequations}
$\sigma_{h}$, $\sigma_{H}$, $\sigma_{s}$ are the terms to be determined during a computation. Due to their nature, above definitions do not impose any constraints on the equations. The definition of the two vertices during computations is given below:
\begin{subequations} \label{vers:eq}
\begin{align}
\Gamma_{1}(u,v)=(1+A) (1+\alpha) (1+a) \tilde{\Gamma}_{1}(u,v) \\
\Gamma_{2}(u,v)=(1+A) (1+\alpha) (1+a) \tilde{\Gamma}_{2}(u,v)
\end{align}
\end{subequations}
Hence, the DSEs for the three field propagators become
\begin{equation} \label{hfdse:eq}
\begin{split}
D^{-1}_{h}(p)=(1+\alpha) [\ p^{2} + \frac{m^{2}_{h,r}}{(1+\alpha)} + (\lambda_{1}+C_{1})(1+A)(1+a) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{1}(-p,q) D_{h}(q-p) ]\
\end{split}
\end{equation}
\begin{equation} \label{Hfdse:eq}
\begin{split}
D^{-1}_{H}(p)=(1+a) [\ p^{2} + \frac{m^{2}_{H,r}}{1+a} + (\lambda_{2}+C_{2})(1+A)(1+\alpha) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{2}(-p,q) D_{H}(q-p) ]\
\end{split}
\end{equation}
\begin{equation} \label{sfdse:eq}
\begin{split}
D^{-1}_{s}(p)=(1+A) [\ p^{2} + 2 (1+a) (1+\alpha) \sigma_{s} + (\lambda_{1}+C_{1})(1+a)(1+\alpha) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{h}(q) \tilde{\Gamma}_{1}(q,-p) D_{h}(q-p) + (\lambda_{2}+C_{2})(1+a)(1+\alpha) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{H}(q) \tilde{\Gamma}_{2}(q,-p) D_{H}(q-p) ]\
\end{split}
\end{equation}
where in equation \ref{hfdse:eq} the renormalized mass for the SM Higgs ($ m_{h,r} $) is fixed at 125.09 GeV during the entire study, while the renormalized mass of the second Higgs boson is fixed during each computation \footnote{The definition of the squared physical mass of the SM Higgs is $m^{2}_{h,r}=m^{2}_{h}+\beta$, and the definition of the squared renormalized mass of the second Higgs is $m^{2}_{H,r}=m^{2}_{H}+b$.}. Equations \ref{hfdse:eq}-\ref{sfdse:eq} are the three DSEs considered for the study. The bare mass of the scalar field is set to zero in \ref{sfdse:eq} in order to investigate the dynamical masses of the singlet scalar field. Lastly, the quantities $\tilde{\Gamma}_{1}(u,v)$ and $\tilde{\Gamma}_{2}(u,v)$ are fixed at $\lambda_{1}$ and $\lambda_{2}$, respectively \cite{Roberts:1994dr}. However, in the current investigation the vertices can still change depending upon the contributions from the coefficients in the counter terms, see equations \ref{vers:eq}.
\par
For each of the propagators, the following renormalization conditions are used.
\begin{equation} \label{hcond:eq}
D_{h}^{ij}(p)  |_{p^{2}=m^{2}_{h,r}} = \frac{\delta ^{ij}}{p^{2}+m^{2}_{h,r}} |_{p^{2}=m^{2}_{h,r}}
\end{equation}
\begin{equation} \label{Hcond:eq}
D_{H}^{ij}(p)  |_{p^{2}=m^{2}_{H,r}} = \frac{\delta ^{ij}}{p^{2}+m^{2}_{H,r}} |_{p^{2}=m^{2}_{H,r}}
\end{equation}
\begin{equation} \label{scond:eq}
D_{s}(p)  |_{p=1} = \frac{1}{p^{2}} |_{p=1}
\end{equation}
The following two conditions are also imposed to numerically compute the correlation functions and the other quantities which are introduced for the counter terms.
\begin{equation} \label{hleast:eq}
\begin{split}
\int_{-\Lambda}^{\Lambda} (\ -D^{-1}_{h}(p)+ (1+\alpha) [\ p^{2} + \frac{m^{2}_{h,r}}{(1+\alpha)} + (\lambda_{1}+C_{1})(1+A)(1+a) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{1}(-p,q) D_{h}(q-p) ]\ )\ ^{2} dp =0
\end{split}
\end{equation}
\begin{equation} \label{Hleast:eq}
\begin{split}
\int_{-\Lambda}^{\Lambda} (\ -D^{-1}_{H}(p)+(1+a) [\ p^{2} + \frac{m^{2}_{H,r}}{1+a} + (\lambda_{2}+ C_{2})(1+A)(1+\alpha) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{2}(-p,q) D_{H}(q-p) ]\ )\ ^{2} dp =0
\end{split}
\end{equation}
Equations \ref{hleast:eq} - \ref{Hleast:eq} are indeed the implementation of the least square method with errors $E_{1}$ and $E_{2}$ defined below.
\begin{equation} \label{herr:eq}
\begin{split}
E_{1}= \int_{-\Lambda}^{\Lambda} (\ -D^{-1}_{h}(p)+ (1+\alpha) [\ p^{2} + \frac{m^{2}_{h,r}}{(1+\alpha)} + (\lambda_{1}+C_{1})(1+A)(1+a) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{1}(-p,q) D_{h}(q-p) ]\ )\ ^{2} dp
\end{split}
\end{equation}
\begin{equation} \label{Herr:eq}
\begin{split}
E_{2}=\int_{-\Lambda}^{\Lambda} (\ -D^{-1}_{H}(p)+(1+a) [\ p^{2} + \frac{m^{2}_{H,r}}{1+a} + (\lambda_{2}+C_{2})(1+A)(1+\alpha) \\ \int_{-\Lambda}^{\Lambda} \frac{d^{4}q}{(2\pi)^{4}} D_{s}(q) \tilde{\Gamma}_{2}(-p,q) D_{H}(q-p) ]\ )\ ^{2} dp
\end{split}
\end{equation}
With imposition of these constraints, the problem at hand becomes that of optimization in which solutions are sought which satisfy equations \ref{hleast:eq}-\ref{Hleast:eq}.
\par
An additional condition given below is also imposed in order to ensure positivity of renormalized squared dynamical scalar mass and evade unwanted numerical fluctuations which may arise due to the difference between the fixed renormalized masses and the dynamically generated masses, and the ever-present limitation in momentum resolution.
\begin{equation} \label{mass2cond:eq}
\begin{split}
m^{2}_{s,r} = (1 + A) (\ m^{2}_{s} + 2(1+\alpha) (1+a) \sigma_{s} )\ \geq 0
\end{split}
\end{equation}
\par
In order to further suppress numerical fluctuations, the SM Higgs is expanded in the form given below:
\begin{equation} \label{hexp1:eq}
D^{ij}_{h}(p)= \delta^{ij} \frac{1}{c(p^{2}+d+f(p))}
\end{equation}
with $f(p)$ given by
\begin{equation} \label{hexp2:eq}
f(p) = \frac{\displaystyle \sum_{l=0}^{N} a_{l} p^{2l}}{\displaystyle \sum_{l=0}^{N} b_{l} p^{2l}}
\end{equation}
In equations \ref{hexp1:eq}-\ref{hexp2:eq}, $c$, $d$, $a_{l}$, and $b_{l}$ are the coefficients to be determined during a computation. A similar expansion with different coefficients is used for the second Higgs propagator. Beside stability, these expansions are also time efficient while performing renormalization and updating the SM and the second Higgs propagators.
\par
The computation starts with $\sigma_{H}=\sigma_{s}=C_{1}=C_{2}=0$, i.e. with no contribution by the counter terms to the renormalized masses and the two Yukawa couplings. Both Higgs propagators are also rendered their respective tree level structures. For the SM Higgs in equations \ref{hexp1:eq}-\ref{hexp2:eq}, $c=1$ and $d=m^{2}_{h}$ \footnote{$m^{2}_{h}=m^{2}_{h,r}$ is used throughout the study.} while all the other coefficients are zero. A similar setup of coefficients is used for the second Higgs. The terms $1+\alpha$ and $1+a$ are calculated from the renormalization conditions in equations \ref{hcond:eq} and \ref{Hcond:eq}. The scalar propagator assumes the values using the equation \ref{sdse:eq} and the quantity $1+A$ is calculated by the renormalization condition \ref{scond:eq} \footnote{The scalar propagator is calculated without the term $1+A$ and then then renormalization condition sets the value of the term $1+A$.}.
\par
An iteration involves updating of the correlation functions and parameters. During an iteration, first, the $\sigma_{s}$, $C_{1}$, and $C_{2}$ are updated as in the mentioned ordered. The update of each of these quantities is performed using Newton Raphson's method with the criteria imposed by the least square method in equations \ref{hleast:eq}-\ref{Hleast:eq}. The updated value is accepted only when both of the errors $E_{1}$ and $E_{2}$ reduce, see equations \ref{herr:eq}-\ref{Herr:eq}.
\par
It is followed by the SM Higgs propagator for which the coefficients in equations \ref{hexp1:eq}-\ref{hexp2:eq} are updated with the above mentioned criteria of acceptance. Upon each change, the SM Higgs propagator is calculated from equations \ref{hexp1:eq}-\ref{hexp2:eq} and renormalized using equation \ref{hcond:eq}.
\par
Lastly, the second Higgs propagator is updated using the same procedure as described above for the SM Higgs, but using equation \ref{Hcond:eq} for renormalization. Upon every change, the scalar propagator is calculated from equation \ref{sfdse:eq} and renormalized using equation \ref{scond:eq}, as mentioned earlier.
\par
A computation concludes only when there are either no further improvements in the quantities such that $E_{1}$ and $E_{2}$ are further reduced, or both of these errors are equal or below the preset value of tolerance. The value of tolerance is set at $10^{-20}$.
\par
Gauss quadrature algorithm is used for numerical integration in the DSEs. The algorithms are developed in C++ environment.
\par
The solutions presented are unique in the sense that order of updates performed on the propagators and other quantities practically does not effect the quantities being computed.
\par
It is assumed that the model is not trivial. The assumption is supported by the fact that, despite that $\phi^{4}$ theory \cite{Hasenfratz:1988kr,Gliozzi:1997ve,Weber:2000dp} is found trivial \cite{Jora:2015yga,Aizenman:1981zz,Weisz:2010xx,Siefert:2014ela,Hogervorst:2011zw}, Higgs interaction with gauge bosons does not render the model trivial \cite{Maas:2013aia,Maas:2014pba}.
\par
A peculiar feature of the model is the possibility of negative Hamiltonian due to certain paths \cite{Rivers:1987hi}. However, it has been noticed in a different investigation of scalar interactions \cite{Mufti:2022abc,Ruthe:2008rut} that in the presence of SM Higgs mass the corresponding histories do not contribute during simulations \footnote{In the potential part of the model, the dominant contribution comes from the term containing the SM (squared) Higgs mass ($125.09$ GeV) while only $\phi(x)$ may force the potential term below zero. Hence, $\phi(x)$ must take significantly large values at all points of the spacetime to effectively compete the term containing $m_{h}$. Numerically, it is not favored in Monte Carlo simulations due to precision related issues and that the fields are rendered values from certain (usually Gaussian) distributions. Another argument is that a Monte Carlo simulation proceeds towards smaller action values which, depending upon the model, are favored by the fields with relatively lower values.}. Presence of large $m_{H}$ in the model further diminishes the possibility that the potential term will cause negative Hamiltonian. Hence, it is assumed here that the results are not effected by the above mentioned feature of the potential term.
\par
Once the scalar dynamical masses in the parameter space are calculated, knowledge from machine learning \cite{Alpaydin14} is employed to represent the results in an attempt to understand how scalar dynamical masses manifest in the parameter space and estimate critical coupling value if possible. There are three free parameters in the model. Hence, scalar dynamical mass is expanded in terms of the second Higgs mass $m_{H}$, and the two Yukawa couplings ($\lambda_{1}$ and $\lambda_{2}$) as given below:
\begin{equation} \label{sDMGexpansion:eq}
m_{s,f} (m_{H},\lambda_{1},\lambda_{2}) = \sum_{i=0}^{i=6} \sum_{j=0}^{j=6} \sum_{k=0}^{k=6} a_{ijk} m_{H}^{i} \lambda_{1}^{j} \lambda_{2}^{k}
\end{equation}
where $m_{s,f}$ is the function representing scalar dynamical mass in the parameter space. The error value $E_{f}$ is described by the following expression:
\begin{equation} \label{sdmgerror:eq}
E_{f} =\frac{1}{N}\sqrt{\sum_{i=0}^{i=N}(m_{s,r} - m_{s,f})_{i}^{2}}
\end{equation}
where $N$ is the number of explored points in the parameter space. The model is studied on 64 points for each of the cutoff values set at $10$ TeV and $100$ TeV. The procedure of training at a fixed cutoff value starts with all weights $a_{ijk}$ set to zero. The results from the DSE calculations are successively introduced to the algorithm from higher to lower couplings, hence the results from DSE calculations being examined increase on every introduction. The weights are varied using Newton Raphson's method \footnote{Other methods are also used for this purpose, see \cite{Alpaydin14} for instance.} such that the error value in equation \ref{sdmgerror:eq} decreases. If the value does not decrease, the change in the weight is not accepted. All weights are examined during each iteration. Once, the algorithm has reached the point when it either can not further improve the weights or the error value has reached its tolerance, further results are introduced to the algorithm and weights are reexamined, hence the algorithm is retrained. The tolerance is set at the value $10^{-16}$. In short, scalar dynamical mass is fitted by the expansion in equation \ref{sDMGexpansion:eq} while successively introducing the results from DSE calculations and readjusting the weights $a_{ijk}$. Throughout the training, $0 \leq m_{s,f} (m_{H},\lambda_{1},\lambda_{2}) $ is imposed. In order to remove any personal bias, no value of scalar dynamical mass was taken as outliers.
\par
As is commonly known, training algorithms on finite data (and duration) may not render extreme accuracy. In addition, the scalar dynamical masses and other calculated quantities (from DSEs) may also suffer from numerical fluctuations. However, a relatively smoother description of the data immensely helps in understanding the overall features of the model and making nontrivial estimates. Hence, the algorithm training was employed to assist in our conclusions regarding the critical coupling value in the model.
\section{Field Propagators}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{h1prs1.eps}
\caption{\label{fig:h1prs1} The SM Higgs propagators with $m_{H}=0.001$ GeV and $m_{H}=1.0$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{h1prs2.eps}
\caption{\label{fig:h1prs2} The SM Higgs propagators with $m_{H}=100$ GeV and $m_{H}=1000$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{h2prs1.eps}
\caption{\label{fig:h2prs1} The second Higgs propagators with $m_{H}=0.001$ GeV and $m_{H}=1.0$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{h2prs2.eps}
\caption{\label{fig:h2prs2} The second Higgs propagators with $m_{H}=100$ GeV and $m_{H}=1000$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{sprs1.eps}
\caption{\label{fig:sprs1} The scalar propagators with $m_{H}=0.001$ GeV and $m_{H}=1.0$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{sprs2.eps}
\caption{\label{fig:sprs2} The scalar propagators with $m_{H}=100$ GeV and $m_{H}=1000$ GeV are plotted for cutoff values at $10$ TeV and $100$ TeV. The parameters in the legend are given as $(m_{H},\lambda_{1},\lambda_{2})$ with all of the parameters mentioned in GeV. For the same cutoff, every two consecutive propagators are $1.0$ TeV apart on the momentum axis for the sake of clarity.}
\end{figure}
Unlike the vertices in the model, the field propagators are rendered considerable freedom to have momentum dependence while conforming to the renormalization conditions. Scalar singlet field propagator has the highest freedom since it takes contributions from interactions with both Higgs field. However, such freedom may also allow contributions from numerical fluctuations of both Higgs propagators as well as other involved quantities in the model.
\par
The SM Higgs propagators are given in figures \ref{fig:h1prs1}-\ref{fig:h1prs2}. An immediate observation is that cutoff effects are not significant throughout the explored parameter space, particularly for the cases with $\lambda_{1} < 1$. The deviations are mostly due to the multiplicative constant which is computed using the renormalization condition \ref{hcond:eq}.
\par
In the parameter space with $\lambda_{1}$ as low as $10^{-6}$, there is no significant dependence of the propagators. However, for higher couplings slight enhancement is observed. Considerable changes occur only at $\lambda_{1}=1.0$. The propagators are found to have similar qualitative behavior for all the renormalized masses of the second Higgs which indicates that the second Higgs field does not influence the SM propagators which could have been possible through the scalar propagators. Overall, there is no significant changes in the SM Higgs propagators in the parameter space of the model. Hence, it is deduced that if there is any non-trivial structure in the phase space, it does not influence the SM Higgs propagators \cite{Maas:2013aia} and possibly the other field propagators since such effects can also translate to other quantities by the coupled DSEs.
\par
In the model, the second Higgs field differs from the SM Higgs field due to different renormalized masses and couplings. Hence, particularly for masses in the vicinity of the SM Higgs mass, it is expected to have similar behavior within numerical fluctuations. The propagators are shown in figures \ref{fig:h2prs1}-\ref{fig:h2prs2}. There are no considerable cutoff effects, as is the case with the SM Higgs propagators. However, a certain dependence on $m_{H}$ is observed in the propagators which is relatively stronger for $m_{H} < m_{h}$ and weakens as $m_{H}$ approaches $m_{h}$. The two Higgs propagators have similar behavior for $m_{H} \simeq m_{h}$ which validates the implementation of the algorithms. For $m_{H} = 1$ TeV  dependence on couplings is lost since the bare mass of the second Higgs dominates the contributions in the propagators, hence, rendering it a tree level structure up to the renormalization dependent term. The propagators are suppressed for such a large renormalized mass which is expected for the case of dominant tree level contribution in the propagator.
\par
The scalar singlet propagators are shown in figures \ref{fig:sprs1}-\ref{fig:sprs2}. As argued in above, cutoff effects are evident on the propagators \footnote{There may also be contributions from the numerical interpolation performed during renormalization of the scalar propagators. In order to suppress these effects, resolution in momentum is kept at its highest feasible value.}. The overall behavior is that for higher couplings, particularly $\lambda_{2}$ since $m_{h}$ is kept fixed, the propagators are enhanced depending upon the second Higgs mass $m_{H}$. As $m_{H}$ increases, this effect tends to disappear in favor of a tree level structure, as is the case for the other field propagators. The scalar propagators suffer strongly from the cutoff effects despite that both Higgs propagators are not effected to such an extent. It implies that the cutoff effects must show up in at least one of the other calculated quantities other than the two Higgs propagators.
\section{Dynamical Scalar Masses}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{smasses10TeV.eps}
\caption{\label{fig:smasses10TeV} Dynamical scalar masses are plotted against second Higgs mass $m_{H}$ for various couplings at $10$ TeV cutoff. The couplings (in GeV) are shown in the caption as $\lambda_{1},\lambda_{2}$. Every consecutive pair of couplings is slightly displaced along x-axis for clarity. The error bars represent difference between the values of scalar dynamical masses obtained by computation (plotted values) and by training of the algorithms.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{smasses100TeV.eps}
\caption{\label{fig:smasses100TeV} Dynamical scalar masses are plotted against second Higgs mass $m_{H}$ for various couplings at $100$ TeV cutoff. The couplings (in GeV) are shown in the caption as $\lambda_{1},\lambda_{2}$. Every consecutive pair of couplings is slightly displaced along x-axis for clarity. The error bars represent difference between the values of scalar dynamical masses obtained by computation (plotted values) and by training of the algorithms.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{fittedsDMG10.eps}
\caption{\label{fig:smassfit10TeV} Dynamical scalar masses obtained by training algorithms are shown against the two couplings $\lambda_{1}$ and $\lambda_{2}$, for various second Higgs masses $m_{H}$ in GeV and at cutoff $\Lambda=10$ TeV indicated on the caption.}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{fittedsDMG100.eps}
\caption{\label{fig:smassfit100TeV} Dynamical scalar masses obtained by training algorithms are shown against the two couplings $\lambda_{1}$ and $\lambda_{2}$ for various second Higgs masses at 100 TeV, for various second Higgs masses $m_{H}$ in GeV and at cutoff $\Lambda=10$ TeV indicated on the caption.}
\end{figure}
Since the scalar propagator also contains its dynamically generated mass, the arguments related to the cutoff effects and numerical artifacts translate into the dynamical masses. However, there are certain vivid features in the masses shown in the figures \ref{fig:smasses10TeV}-\ref{fig:smasses100TeV}.
\par
Considering the dynamical masses for 100 TeV cutoff in figure \ref{fig:smasses100TeV}, the extent of the mass production is relatively higher for the second Higgs mass in MeVs and high second coupling for a number of cases. For very small couplings, large difference between the computed value and $m_{s,f}$, possibly due to numerical fluctuations, suggests that these two points could have been taken as outliers. Hence, it is concluded with confidence that the dynamically generated scalar mass in the model is restricted well within 200 MeVs. It is clear that the Yukawa interactions in the scalar sector have a shortcoming when it comes to producing dynamical masses in GeVs.
\par
It is also evident from figures \ref{fig:smasses10TeV}-\ref{fig:smasses100TeV} that cutoff effects are far less severe for the cutoff in hundreds of TeVs. For higher cutoff, the model produces scalar masses with less dependence on the second Higgs mass. Furthermore, the production of scalar mass is significantly lower if one of the couplings is as low as $10^{-6}$ GeV, unless the other coupling reaches $1.0$ GeV, see figure \ref{fig:smasses100TeV}. Hence, the model does posses a demarcation over the coupling values in the vicinity of $1.0$  and $10^{-6}$ GeV. For the case of both couplings at $10^{-6}$ GeV, scalar mass practically looses dependence on the two couplings. It is taken as a sign of existence of critical coupling value between $10^{-3}$ GeV and $10^{-6}$ GeV.
\par
An interesting feature in the low coupling region is that for higher cutoff values the model produces smaller masses which also tend to decrease with the two couplings. One may expect that the behavior persists until the critical coupling value. If this is indeed the case, the model may be useful to study ultralight scalar interactions \cite{Rindler-Daller:2013zxa,Harko:2014vya,Chavanis:2011zi,Huang:2013oqa,Marsh:2015xka,Hui:2016ltb,Primack:2009jr,Hu:2000ke}, though it may undoubtedly be a daunting task from the perspective of numerical precision.
\par
Scalar dynamical masses obtained by the trained algorithm are plotted in figures \ref{fig:smassfit10TeV}-\ref{fig:smassfit100TeV}. The weights of the expansion in equation \ref{sDMGexpansion:eq} are given in appendix A. An immediate observation is existence of a unique value of the critical coupling in the model below $10^{-3}$ GeV. Determining exact value is hampered by the ever-present limitation of the data in machine learning. However, it is evident that the model strongly favors a critical coupling in the region $10^{-6} < \lambda_{i} < 10^{-3}$ GeV.
\par
Despite the limitations in training of algorithms and numerical fluctuations in the exposed data, the weights $a_{ijk}$ are found to have a particular behavior. Firstly there are strong contributions for the terms with $k = 0$ and $i,j \neq 0$, while there is practically no contributions by $a_{ijk}$ for $i,j=0$ irrespective of the value of $k$. Furthermore, increasing the cutoff suppresses $a_{ijk}$ for most of the $i$ and $k$ values at $j=0$. Since the model has significant cutoff effects which becomes milder as the cutoff is raised, it is expected that scalar dynamical mass may not receive significant contributions from terms involving $j=0$ at high cutoff values. Since the mass of the SM Higgs boson is kept fixed throughout the investigation, these observations suggest that the (magnitude of) contributions to dynamical scalar mass favor interactions with the SM Higgs boson mass.
\section{Conclusion}
In the presence of an SM Higgs, the two Higgs doublet model has the capacity to dynamically produce scalar mass significantly larger than of the lightest quarks and leptons. As the cutoff is raised above $100$ TeV, stability in the mass with magnitude below $200$ MeV ensues against the second Higgs mass. It is expected that at a cutoff much higher than $100$ TeV, a single value of dynamical mass is favored by the model irrespective of the second Higgs mass.
\par
However, the dynamical mass has sensitivity to the couplings which diminishes as the couplings are reduced below $10^{-3}$ GeV, with the minimum mass being produced. It strongly implies existence of a critical coupling between $10^{-3}$ GeV and $10^{-6}$ GeV. It presents an opportunity to investigate the model to understand new physics involving the particles considerably lighter than $1$ GeV, such as ultra-light scalars. At the same time, it invites for further study of richer models which contain higher renormalizable vertices in an attempt to explore existence of critical couplings and extent of mass generation due to diversity in interactions.
\par
The role of cutoff effects can not be neglected in the model. Scalar propagator and scalar dynamical mass suffer the most, while for masses larger than $100$ GeV the field propagators are relatively less effected as was observed in the two Higgs propagators. Since dynamical masses are less than $1$ GeV, numerical fluctuations and the cutoff effects hamper in finding an accurate description of how the masses behave in the parameter space. However, considering sufficient number of points in the parameter space a mathematical description for dynamical masses was still found which concurs with the deduction regarding the critical coupling in the model.
\par
A model with the capacity to dynamically produce masses considerably larger than the lightest quarks and leptons can certainly not be ruled out in scalar interactions. The study invites extensions, such as by including richer interactions or other fields, to have a better understanding of how (particularly light) scalars play role at the fundamental level in our universe.
\section{Acknowledgments}
This work was supported by Lahore University of Management Sciences, Pakistan for developing the algorithms and performing computations.
\section{Appendix A} The values of the weights $a_{ijk}$ in equation \ref{sDMGexpansion:eq} for cutoff values at $\Lambda=10$ TeV and $\Lambda=100$ TeV are given below:
\begin{center}
\begin{longtable}{ | m{0.5cm} | m{0.5cm}| m{0.5cm} | m{3.5cm} | m{3.5cm} | }
\hline
$i$ & $j$ & $k$ & $a_{ijk}$ (10 TeV) & $a_{ijk}$ (100   TeV) \\
\hline
0 & 0 & 0 & 0.062470119000000 & -0.005102152199999 \\
\hline
1 & 0 & 0 & 0.432639286999998 & 0.067154687499999 \\
\hline
0 & 1 & 0 & 0.431424060300000 & 0.276498992400000\\
\hline
0 & 0 & 1 & 0.000127544300000 & 0.000023271900000\\
\hline
2 & 0 & 0 & -0.196267245799999 & 0.826964182800007\\
\hline
1 & 1 & 0 & 3.483278381999988 & -2.646594633199990\\
\hline
1 & 0 & 1 & 0.001256076799999 & 0.000786048000000\\
\hline
0 & 2 & 0 & -2.426385300799987 & -0.005208820599999\\
\hline
0 & 1 & 1 & -0.001780857000000 & -0.000014844700000\\
\hline
0 & 0 & 2 & -0.0000001083 & -0.000000023100000\\
\hline
3 & 0 & 0 & -0.595301560000001 & -0.158007649599999\\
\hline
2 & 1 & 0 & 0.311572103699999 & 5.428065021000011\\
\hline
2 & 0 & 1 & -0.001414899099999 & -0.002876003599999\\
\hline
1 & 2 & 0 & -2.967471651800008 & -3.556474306000007\\
\hline
1 & 1 & 1 & -0.011901846299999 & 0.013548386599999\\
\hline
1 & 0 & 2 & -0.0000018441 & -0.000000603000000\\
\hline
0 & 3 & 0 & 3.360127912999990 & -0.928262618399999\\
\hline
0 & 2 & 1 & 0.004300518299999 & -0.002201626699999\\
\hline
0 & 1 & 2 & 0.000000999000000 & -0.000000020900000\\
\hline
0 & 0 & 3 & 0.0 & 0.0\\
\hline
4 & 0 & 0 & 1.045801962500008 & -0.402565733399999\\
\hline
3 & 1 & 0 & -1.420411216799999 & -0.572378273000001\\
\hline
3 & 0 & 1 & 0.003927822499999 & -0.000554230599999\\
\hline
2 & 2 & 0 & -1.122774422700004 & -0.529377973400000\\
\hline
2 & 1 & 1 & 0.000976498800000 & -0.005964260699999\\
\hline
2 & 0 & 2 & 0.0000004915 & 0.000001401500000\\
\hline
1 & 3 & 0 & 0.322474038600001 & 3.267544194000007\\
\hline
1 & 2 & 1 & 0.005696769799999 & -0.001591530000000\\
\hline
1 & 1 & 2 & 0.0000093254 & -0.000008751100000\\
\hline
1 & 0 & 3 & 0.0 & 0.0\\
\hline
0 & 4 & 0 & -3.398481779199986 & 1.658322491499994\\
\hline
0 & 3 & 1 & -0.002802221199999 & 0.001725826300000\\
\hline
0 & 2 & 2 & -0.0000000495 & 0.000000259100000\\
\hline
0 & 1 & 3 & 0.0 & 0.0\\
\hline
0 & 0 & 4 & 0.0 & 0.0\\
\hline
5 & 0 & 0 & -0.831871660400005 & 0.132814269100000\\
\hline
4 & 1 & 0 & 0.949120657999995 & -1.788138425800004\\
\hline
4 & 0 & 1 & -0.003439259299999 & 0.001322541799999\\
\hline
3 & 2 & 0 & -1.903333068200005 & -4.520240236500010\\
\hline
3 & 1 & 1 & -0.005681518600000 & 0.006113348800000\\
\hline
3 & 0 & 2 & 0.0000018543 & 0.000000118100000\\
\hline
2 & 3 & 0 & 5.737533950000011 & 2.994909360000007\\
\hline
2 & 2 & 1 & 0.004100454999999 & -0.003434122199999\\
\hline
2 & 1 & 2 & 0.0000045994 & -0.000004740900000\\
\hline
2 & 0 & 3 & 0.0 & 0.0\\
\hline
1 & 4 & 0 & -4.607335835299982 & 3.560205072699987\\
\hline
1 & 3 & 1 & -0.0017640611 & 0.000632370000000\\
\hline
1 & 2 & 2 & -0.000002936 & 0.000001718100000\\
\hline
1 & 1 & 3 & 0.0 & 0.0\\
\hline
1 & 0 & 4 & 0.0 & 0.0\\
\hline
0 & 5 & 0 & 2.476695333800004 & -1.438056151699995\\
\hline
0 & 4 & 1 & -0.001170340799999 & 0.000608526900000\\
\hline
0 & 3 & 2 & -0.0000001277 & 0.000000238100000\\
\hline
0 & 2 & 3 & 0.0 & 0.0\\
\hline
0 & 1 & 4 & 0.0 & 0.0\\
\hline
0 & 0 & 5 & 0.0 & 0.0\\
\hline
6 & 0 & 0 & 0.250543694400000 & -0.303845693199999\\
\hline
5 & 1 & 0 & -0.872865771300006 & -0.570790623500001\\
\hline
5 & 0 & 1 & -0.0014202491 & -0.000099293799999\\
\hline
4 & 2 & 0 & -2.629992650700005 & -4.544265494000010\\
\hline
4 & 1 & 1 & -0.002726043499999 & 0.004424348300000\\
\hline
4 & 0 & 2 & 0.0000007263 & 0.000000419100000\\
\hline
3 & 3 & 0 & 4.539023046400008 & -1.631924930300001\\
\hline
3 & 2 & 1 & 0.004210737 & 0.000127814799999\\
\hline
3 & 1 & 2 & 0.0000038225 & -0.000004546200000\\
\hline
3 & 0 & 3 & 0.0 & 0.0\\
\hline
2 & 4 & 0 & -2.361741260400014 & 7.959941882000007\\
\hline
2 & 3 & 1 & 0.0004258372 & -0.003437821899999\\
\hline
2 & 2 & 2 & -0.0000105292 & 0.000010411400000\\
\hline
2 & 1 & 3 & 0.0 & 0.0\\
\hline
2 & 0 & 4 & 0.0 & 0.0\\
\hline
1 & 5 & 0 & 2.65929292970001 & -3.122536580600026\\
\hline
1 & 4 & 1 & 0.007337953699999 & -0.008389909299999\\
\hline
1 & 3 & 2 & -0.0000050914 & 0.000004083200000\\
\hline
1 & 2 & 3 & 0.0 & 0.0\\
\hline
1 & 1 & 4 & 0.0 & 0.0\\
\hline
1 & 0 & 5 & 0.0 & 0.0\\
\hline
0 & 6 & 0 & -0.408070357499999 & 0.605073464699998\\
\hline
0 & 5 & 1 & 0.0012935581 & -0.001141694000000\\
\hline
0 & 4 & 2 & -0.0000006819 & 0.000000429700000\\
\hline
0 & 3 & 3 & 0.0 & 0.0\\
\hline
0 & 2 & 4 & 0.0 & 0.0\\
\hline
0 & 1 & 5 & 0.0 & 0.0\\
\hline
0 & 0 & 6 & 0.0 & 0.0\\
\hline
\end{longtable}
\end{center}
\bibliographystyle{plain}