Multi-agent reinforcement learning using echo-state network and its application to pedestrian dynamics

In environments with several agents, RL agents are required to cooperate or compete with other agents. Such environments appear in various fields, thus, MARL has been studied intensively[17]. In particular, methods of deep reinforcement learning (DRL) have been applied to MARL, as well as single-agent RL[19, 18, 20].

The simplest method to implement MARL is to allow the agents to learn their policies independently. However, independent learners cannot access the experiences of other agents, and consider other agents as a part of environment. This results in the non-stationarity of the environment from the perspective of each agent. To address this problem, many algorithms have been proposed. For example, the multi-agent deep deterministic policy gradient (MADDPG) algorithm adopts the actor-critic method, and allows the critic to observe all agents’ observations, actions, and target policies[27]. This algorithm was improved and applied to the pedestrian dynamics by Zheng and Liu[6]. Gupta et al. introduced a parameter sharing method wherein agents shared the parameters of their neural networks[22]. Their method exhibited high learning efficiency in environments with homogeneous agents.

This study adopted the parameter sharing method because it can be easily extended to a form suitable to ESN. Specifically, we considered one or two groups of agents that shared their experiences and policies but received rewards separately.

As mentioned earlier, ESN is a type of reservoir computing that uses RNN[23, 24, 25]:

where $\mbox{\boldmath$X$}^{\mathrm{in}}(t)$ and $\mbox{\boldmath$X$}^{\mathrm{out}}(t)$ are the input and output vectors at time $t$, respectively, $f$ is the activation function, and $\alpha$ is the constant called leaking rate. In the case of ESN, elements of matrices $\mbox{\boldmath$W$}^{\mathrm{in}}_{s}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$, and $\mbox{\boldmath$W$}^{\mathrm{res}}$ are randomly fixed, and only the output weight matrix $\mbox{\boldmath$W$}^{\mathrm{out}}$ is trained. Through this restriction of the trainable parameters, the computational cost is considerably lowered compared with the deep learning. In addition, the exploding/vanishing gradient problem does not occur unlike the training of RNN using backpropagation. For example, in the case of supervised learning, training is typically conducted by the Ridge regression that minimizes the mean square error between the output vector $\mbox{\boldmath$X$}^{\mathrm{out}}$ and the training data $\mbox{\boldmath$X$}^{\mathrm{train}}$. This calculation is considerably faster than the backpropagation used in the deep learning.

However, in the case of RL the training data should be collected by agents themselves. Hence, the training method is not as easy as in the case of the supervised learning. Many previous studies have proposed different algorithms to realize RL using ESN. For example, Szita et al. used the SARSA method[12], Oubbati et al. adopted ESN as a critic of the actor-critic method[13], Chang and Futagami used a covariance matrix adaptation evolution strategy (CMA-ES)[14], and Zhang et al. improved the least squares policy iteration (LSPI) method[15]. Previous studies have also utilized ESN to MARL. For example, Chang et al. trained ESN employing a method similar to that of a deep Q network(DQN), and applied it to the distribution of the radio spectra of telecommunication equipment[16]. In this study, we adopted the LSPI method, as in Ref. [15]; however, we used a simpler algorithm. Note that the absence of the exploding/vanishing gradient problem in the ESN comes from the point that the weight matrices involved in this problem are fixed. Hence, although the above-mentioned training algorithms including LSPI are more complex than those of supervised learning, such problem still does not appear.

Typically, the recurrent weight $\mbox{\boldmath$W$}^{\mathrm{res}}$ is a sparse matrix with spectral radius $\rho<1$. To generate such a matrix, first, a random matrix $\mbox{\boldmath$W$}_{0}^{\mathrm{res}}$ must be created whose components $W_{0,\mu\nu}^{\mathrm{res}}$ are expressed as:

where $\mathcal{P^{\mathrm{res}}}$ is a predetermined probability distribution and $p_{s}^{\mathrm{res}}$ is a hyperparameter called sparsity. Regularizing this matrix such that its spectral radius is expressed as $\rho$, we obtain $\mbox{\boldmath$W$}^{\mathrm{res}}$:

where $\rho\left(\mbox{\boldmath$W$}_{0}^{\mathrm{res}}\right)$ is the spectral radius of $\mbox{\boldmath$W$}_{0}^{\mathrm{res}}$. According to Ref. [24], the input weights $\mbox{\boldmath$W$}^{\mathrm{in}}_{s}$ and $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$ are usually generated as dense matrices. However, we also rendered them sparse, as has been explained later.

The LSPI method is an RL algorithm that was proposed by Lagoudakis and Parr[26]. It approximates the state-action value function $Q(s,a)$ as the linear combination of given functions of $s$ and $a$, $\mbox{\boldmath$\phi$}(s,a)=\left(\phi_{1}(s,a),...,\phi_{N}(s,a)\right)^{T}$:

where $\mbox{\boldmath$w$}=(w_{1},...,w_{N})^{T}$ is the coefficient vector that is to be trained. Using this value, action $a_{t}$ of the agent at each time step $t$ is decided by the $\epsilon$-greedy method:

where $s_{t}$ is the state at this step. Usually, $\epsilon$ is initially set as a relatively large value $\epsilon_{0}$, and gradually decreases.

To search the appropriate policy, the LSPI method updates $w$ iteratively. Let the $n$-th update be executed after the duration $\mathcal{E}_{n}$. Then the form of $w$ after the $n_{l}$-th duration, $\mathcal{E}_{n_{l}}$, is expressed as:

where

$r_{t}$ is the reward at this step, and $\gamma\in(0,1)$ is the discount factor. We let $\mbox{\boldmath$\phi$}(s_{t+1},a_{t+1})=0$ if the episode ends at time $t$. The constant $\lambda\in(0,1)$ is a forgetting factor that is introduced because the contributions of old experiences collected under the old policy do not reflect the response of the environment under the present policy. This factor has also been adopted by previous studies such as Ref. [15].

In this study, we considered the partially observable Markov decision process (POMDP), a process wherein each agent observed only a part of the state. As the input of the $i$-th agent, we set its observation $\mbox{\boldmath$o$}_{i,t}$, the candidate of the action $a$, and the bias term:

where $\mbox{\boldmath$\eta$}_{a}$ is the one-hot expression of $a$ and activation function $f$ is ReLU. We expressed the number of neurons of the reservoir as $N^{\mathrm{res}}$, i.e. $N^{\mathrm{res}}\equiv\mathrm{dim}\mbox{\boldmath$X$}_{i}^{\mathrm{res}}$. The initial condition of reservoir was set as $\mbox{\boldmath$X$}_{i}^{\mathrm{res}}(-1)=0$. We assumed that all agents shared the common input and reservoir matrices $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$, and $\mbox{\boldmath$W$}^{\mathrm{res}}$. For later convenience, $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$, and $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$, that is, the parts of the input matrix are denoted separately. The action $a_{t}$ was decided by the $\epsilon$-greedy method, Eq. (9), under this $Q_{i}(\mbox{\boldmath$o$}_{i;t},a)$. Here, $Q_{i}$ should be calculated for all possible candidate of action $a$. Using the accepted action, $\mbox{\boldmath$X$}_{i}^{\mathrm{res}}$ is updated as follows:

In the actual calculation of Eq. (13), it is convenient to calculate $\tilde{\mbox{\boldmath$X$}}_{i}^{\mathrm{res}}(t;a)$ for all candidates of action $a$ at once using the following relation:

Here, $|\mathcal{A}|$ is the number of candidates of the action, and $\mathrm{rep}_{|\mathcal{A}|}\left(\mbox{\boldmath$X$}_{t}^{\mathrm{in}}\right)$ is a matrix constructed by repeating the column vector, $\mbox{\boldmath$X$}_{t}^{\mathrm{in}}$, $|\mathcal{A}|$-times. To derive the last line of Eq. (19), we should remind that $\mbox{\boldmath$\eta$}_{a}$ is the one-hot expression:

Comparing Eqs. (8), (17) and (18), we can apply the LSPI method by substituting $\mbox{\boldmath$w$}^{T}=\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$ and $\mbox{\boldmath$\phi$}=\left(\mbox{\boldmath$X$}_{i}^{\mathrm{res}}(t)^{T},1\right)^{T}$. Thus, taking the transpose of Eq.(10), $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$ can be calculated by the following relation

where

and

Here, we divided the agents into one or several groups, and let all members of the same group share their experiences when we calculated $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$. In Eqs. (22), (23) and (24), $G_{i}$ indicates the group to which the $i$-th agent belongs. In Eq. (23), we let $\hat{\mbox{\boldmath$X$}}_{j}^{\mathrm{res}}(t+1)=0$ if the episode ended at time $t$, as in Section 2.3. The second term of the right-hand side of Eq. (23) was introduced to avoid the divergence of the inverse matrix of Eq. (22). The constant $\beta$ is a small hyperparameter, and $I$ is the identity matrix. As the duration $\mathcal{E}_{n}$, we adopted one episode in this study. In other words, the output weight matrix $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$ was updated at the end of each episode. The agents of the same group shared the common output weight matrix by Eq. (22). Considering that all agents had the same input and reservoir matrices, as explained above, agents of the same group shared all parameters of neural networks. Thus, this is a form of parameter sharing suitable for ESN. As mentioned in Section 2.2, Ref. [15] also applied LSPI method to train ESN. They calculated the output weight matrix using the recursive least squares(RLS) method to avoid the calculation of the inverse matrix of Eq. (22), which incurred a computational cost of $O\left((N^{\mathrm{res}})^{3}\right)$. However, we did not adopt RLS method, because an additional approximation called mean-value approximation is required for this method. In addition, when each duration $\mathcal{E}_{n}$ is long and the parameter sharing method is introduced, the effect of the computational cost mentioned above on the total computation time is limited.

In the actual calculation, we made the matrices

and

, and performed the operation

at the end of each episode $\mathcal{E}_{n}$. Here, we used the fact that the reward after the end of the episode is zero: $r_{j;t_{\mathrm{max}}}=0$. However, the term $\hat{\mbox{\boldmath$X$}}_{j}^{\mathrm{res}}(t_{\mathrm{max}})\neq 0$ should not be ignored. After this manipulation, we calculated $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$ using Eq. (22), multiplying forgetting factors to matrices $\tilde{A}_{G}$ and $\tilde{B}_{G}$:

and updated the value of $\epsilon$ used for the $\epsilon$-greedy method as:

if $\epsilon$ was larger than the predetermined threshold value $\epsilon_{\mathrm{min}}$. Training is dependent on the information in past state, action, and reward only through Eqs. (49) and (50). Thus, matrices $Y_{1},Y_{2},$ and $R$ can be deleted after the update of $\tilde{A}_{G}$ and $\tilde{B}_{G}$ using these equations. Note that, in principle, we can update $\tilde{A}_{G}$ and $\tilde{B}_{G}$ by adding the increments every time step even if we do not prepare matrices $Y_{1},Y_{2},$ and $R$. However, the number of calculations of large matrices multiplication should be reduced to lower the computational cost, at least when the calculation is implemented by Python. The algorithm explained above is summarized in Algorithm 1.

We constructed a grid-world environment wherein each cell could have three states: “vacant,” “wall,” and “occupied by an agent.” At each step, each agent chose one of the four adjacent cells: up, down, right, and left, and attempted to move into it. This move succeeded if and only if the candidate cell was vacant and no other agents attempted to move there. Each agent observed $11\times 11$ cells centering on the agent itself as 2-channel bitmap data. Here, agents including observer itself are indicated as pairs of numbers $(1,0)$, walls are as $(0,1)$, and vacant cells are as $(0,0)$. Hence, the observation of $i$-th agent, $\mbox{\boldmath$o$}_{i;t}$, was $11\times 11\times 2=242$-dimensional vector. Although several previous studies have proposed a method that combined ESN with untrained convolutional neural networks when it is used for image recognition[14, 28, 29], we input the observation directly into ESN. Instead, we improved the sparsity of input weight matrix $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$ to reflect the spatial structure. Specifically, we first called the $3\times 3$, $7\times 7$, and $11\times 11$ cells centering on agent as $\mathcal{A}_{1}$, $\mathcal{A}_{2}$ and $\mathcal{A}_{3}$ respectively, as in Fig. 1. Then, we changed the sparsity of the components of $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$ depending on which cell they combined with. Namely, components of $W_{o,\mu\nu}^{\mathrm{in}}$ are generated as:

and the sparsity $p_{s}^{\mathrm{in}}(\nu)$ is expressed as follows:

where the hyperparameters $p_{s1}^{\mathrm{in}}$, $p_{s2}^{\mathrm{in}}$, and $p_{s3}^{\mathrm{in}}$ obey $p_{s1}^{\mathrm{in}}\leq p_{s2}^{\mathrm{in}}\leq p_{s3}^{\mathrm{in}}$. Thus, information on cells near the agent combined with the reservoir more densely than that on cells far from the agent. Further, we let $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$ be a sparse matrix with sparsity $p_{sb}^{\mathrm{in}}$, but made $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$ dense. In addition, $\mbox{\boldmath$W$}^{\mathrm{res}}$ was generated as a sparce matrix with its spectral radius $\rho<1$, by the process explained in Section 2.2. Here, $\mathcal{P}^{\mathrm{in}}_{o}$, $\mathcal{P}^{\mathrm{in}}_{a}$, $\mathcal{P}^{\mathrm{in}}_{b}$, and $\mathcal{P}^{\mathrm{res}}_{0}$, the distribution of nonzero components of $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$, and $\mbox{\boldmath$W$}^{\mathrm{res}}_{0}$ were expressed as gaussian functions $\mathcal{N}\left(0,\left(\sigma^{\mathrm{in}}_{o}\right)^{2}\right)$, $\mathcal{N}\left(0,\left(\sigma^{\mathrm{in}}_{a}\right)^{2}\right)$, $\mathcal{N}\left(0,\left(\sigma^{\mathrm{in}}_{b}\right)^{2}\right)$, and $\mathcal{N}\left(0,\left(\sigma^{\mathrm{res}}_{0}\right)^{2}\right)$, respectively. The list of hyperparameters are presented in Table 1.

As the specific tasks, we considered the following two situations, changing the number of agents, $n_{\mathrm{agent}}$. Note that we did not execute any kinds of pretraining in the simulations.

The learning curves of task. I are shown in Fig. 5. In these graphs, the green curve is the mean value of all agents’ rewards, and the red and blue curves indicate the rewards of the best and worst-performing agents. In each case, we executed 8 independent trials and averaged the values of the graphs, and painted the standard errors taken from these trials in pale colors. As shown in Fig. 5, the performance of each agent appeared to worsen with increasing $n_{\mathrm{agent}}$. This is because each agent was prevented from proceeding by other agents. In the case that $n_{\mathrm{agent}}\geq 24$, the direct route was excessively crowded and certain agents had no choice but to go to the detour, as shown in snapshots and density colormaps in Figs. 6 and 7. (b)-(d). Thus, the difference in the performance between the best and worst agents increased compared to the case where $n_{\mathrm{agent}}=12$. In addition, seeing the density colormaps of Fig. 7, some agents tend to get stuck in the bottom-right corner of the detour when $n_{\mathrm{agent}}$ is large. Considering that there is no benefit for agents in staying at this corner, it is thought to be the misjudgment resulting from the low information processing capacity of ESN.

From the definition of reward in this study, the total reward that one agent gets during one episode is equal to the displacement along the direction it tries to go. Hence, the average velocity of agents, $\bar{v}$, can be calculated as

Note that the upper bound of this value is 1, which can be realized when all agents always proceed in the direction they want to. In addition, considering that simulations of this study impose the periodic boundary condition and the number of agents $n_{\mathrm{agent}}$ does not change, the average density $\bar{\rho}$ is expressed as

Here, note that the phrase “the number of cells that agents can walk into” means the area that is not devided from agents by walls, and does not mean the number of cells that a certain agent can move into at a certain step. In the case of task. I, for example, this number can be counted as 192, by Fig. 2. Using these equations, we can plot the fundamental diagram[30], the graph that draw the relation between velocity $\bar{v}$ and density $\bar{\rho}$, as Fig. 8. These data are averaged over the 151-250 episodes of 8 independent trials. In this graph, we estimated the error bars by calculating the standard error of the mean over 8 independent trials, but these error bars are smaller than the points. Seeing Fig. 8, $\bar{v}$ has the value near its upper bound when $\bar{\rho}$ is small, and gradually decreases with increasing $\bar{\rho}$. This decrease itself is observed in many situations of pedestrian dynamics. However, considering that agents get stuck in the corner under present calculation, the improvement of the algorithm and neural network is thought to be required before quantitative comparisons with previous studies.

The learning curves of task. II are shown in Fig. 9. Here, as discussed in Section 4.1, we averaged the values of the rewards over 8 independent trials, and painted the standard errors among these trials in pale colors. According to Fig. 9, the training itself was successful except when $n_{\mathrm{agent}}=64$. The corresponding snapshots and density colormaps shown in Figs. 10 and 11 show that in the case of $n_{\mathrm{agent}}\leq 48$, agents of the same group created lanes and avoided collisions with agents in the opposite direction. Such lane formation was also observed in various previous studies that used mathematical models[2], experiments[21, 31, 32], and other types of RL agents [4, 8]. Note that the lanes had different shapes depending on the seed of random values, i.e. the right-proceeding agents walked on the upper side of the corridor in some trials, and on the lower side in others, for example. Hence, the colormap of Fig. 11 indicates the time average of one representative trial. To show the dependency on the random seed, we also drew the colormaps of the case that $n_{\mathrm{agent}}=32$ with different seeds as Fig. 12.

When $n_{\mathrm{agent}}=64$, two groups of agents failed to learn how to avoid each other and collided. In this case, each agent cannot move forward after the collision, and the reward they obtained was limited to that for the first few steps. Note that in this task, each agent cannot access the experiences of agents of the other group. Thus, thinking from the perspective of MARL, there is a possibility that the non-stationarity of the environment (in the viewpoint of agents) worsened the learning efficiency, as with the case of the independent learners. However, from the perspective of soft matter physics, this stagnation of the movement under high density resembles the jamming transition observed in wide range of granular or active matters[3, 33, 34, 35]. Nevertheless, to investigate whether it is really related with the jamming transition, simulations with larger scale are required.

As in task. I, we can plot the fundamental diagram shown in Fig. 13 using Eqs. (56) and (57). This graph shows the sudden decrease of the average velocity between $n_{\mathrm{agent}}=48$ and 64 (i.e. $\bar{\rho}=0.3$ and 0.4), reflecting the stagnation of the movement discussed above. Similar behavior is also reported in previous researches that considered the jamming transition in mathematical models of pedestrians[3].

As explained in Section 3.1, we adopted the parameter sharing method in our calculations. In this section, we calculate the two tasks discussed above considering the case that the parameter sharing was not adopted, and compare the results with those of previous sections to investigate the effect of this method. Specifically, equations

and

were used to calculate $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$ in this section. Thus, each agent learned independently. To observe the effect of replacing Eq.(22) by Eq.(58), input and reservoir weight matrices, $\mbox{\boldmath$W$}^{\mathrm{in}}_{o}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$, $\mbox{\boldmath$W$}^{\mathrm{in}}_{b}$, and $\mbox{\boldmath$W$}^{\mathrm{res}}$, were still shared by all agents, and the hyperparameters were the same as listed in Table 1.

The learning curves of tasks. I and II are shown in Figs. 14 and 15, respectively. In these graphs, scores when using the parameter sharing method averaged over the 151–250 episodes of 8 independent trials are drawn as the dotted lines. A comparison of the learning curves with these dotted lines revealed that the performance when the parameter sharing was not used was almost the same as or slightly lower than that under the parameter sharing. Hence, the parameter sharing method did not result in the drastic improvement of the learning. The advantage of this method is rather the reduction in the number of inverse matrix calculations, at least in the case of our tasks. We actually calculated the case of task. II under $n_{\mathrm{agent}}=64$, but did not plot the graph because agents failed to make lanes like the case of Section 4.2.

As we explained in Section 3.2.2, in the case of task. II, we divided the pedestrians into two groups by the direction they want to proceed in. Here, the training of each group was executed independently. However, we can also implement the MARL agents even in the case that the parameters of neural networks are shared between two groups. Specifically, we first define a two-dimensional one-hot vector $\mbox{\boldmath$\eta$}_{g}=(1,0)^{T}$ or $(0,1)^{T}$, which represents the information on each agent’s group. Adding this vector to agents’ observation and modifying eq. (13) as follows:

the learning process can be executed even if parameters are shared between two groups. To share the parameters, equations

and

are used to calculate $\mbox{\boldmath$W$}_{i}^{\mathrm{out}}$, instead of eqs. (22), (23), and (24). Here, we let $\mbox{\boldmath$W$}^{\mathrm{in}}_{g}$, the input weight matrix corresponding to $\mbox{\boldmath$\eta$}_{g}$, be a dense matrix, and used the same standard deviation as $\mbox{\boldmath$W$}^{\mathrm{in}}_{a}$. Namely, the distribution of each component of $\mbox{\boldmath$W$}^{\mathrm{in}}_{g}$ is $\mathcal{N}\left(0,\left(\sigma^{\mathrm{in}}_{a}\right)^{2}\right)$. In this section, we investigate whether this method work. The hyperparameters themselves were the same as listed in Table 1.

The result is shown in Fig. 16. In this figure, the dotted lines are the scores of Section 4.2, the case that two groups do not share their parameters, averaged over the 151–250 episodes of 8 independent trials. Seeing Fig. 16 (a) and (b), in the case that $n_{\mathrm{agent}}=32$ and 48, the parameter sharing between two groups slightly improved or hardly changed the mean and best scores compared with those under the training without parameter sharing between different groups, whereas the worst score was lowered. It is thought that the additional information, $\mbox{\boldmath$\eta$}_{g}$, led to the confusion of the neural network and wrong action choice, especially under subtle situations that is difficult to evaluate $Q_{i}(\mbox{\boldmath$o$}_{i;t},a)$. The result under $n_{\mathrm{agent}}=64$ shown in Fig. 16 (c) is interesting. As we saw in Section 4.2, agents that did not share their parameters between different groups failed to learn to make lanes in all trials in this case. However, agents of this section succeeded in making lanes in some trials. We can think two possible factors that cause this difference. One is the number of agents choosing wrong actions. Blockage of agents in this task is caused when they learn wrong strategy that “go straight and earn rewards for first few steps”, before they learn to make lanes. As we discussed above, agents of this section are thought to be sometimes confused by the additional information, $\mbox{\boldmath$\eta$}_{g}$. Even if agents learn the wrong strategy, blockage is eased compared with that of Section 4.2 because such “confused” agents do not simply go straight. The other possible factor is the sharing of the experience. Agents of this section share the experience between different groups by the parameter sharing, whereas those of Section 4.2 share it only among the same group. This experience sharing is thought to promote the efficient learning of correct strategy that “agents of different groups should go to different lanes”. It is difficult to identify which of these two factors plays the significant role, only by this task.

In this section, we calculated the same tasks using representative DRL algorithms, DQN[9], advantage actor-critic(A2C)[10], and proximal policy optimization (PPO)[11] methods, and compare the performance with our algorithm. In these calculations, we prepared the neural network with the same size as the reservoir of our ESN, i.e., the network has 1 hidden layer (of perceptron) with 1024 neurons. Each agent observed 2-channel bitmap data with $11\times 11$ cells centering on itself and agents of the same group shared the experience by parameter sharing, as in our algorithm. Parameter sharing between different groups like Section 4.4 is not introduced. Note that in the A2C and PPO methods, actor- and critic-networks shared the parameters except for the output weight matrices.

As the optimizer, we adopted Adam[36]. The codes were implemented using Pytorch (torch 2.2.0), and we used the default values of this library for the initial distribution of the parameters of neural network and the hyperparameters of the optimizer except for the learning rate. The parameters related to environment or reward, such as $t_{\mathrm{max}}$ and $\gamma$, were the same as our method. In addition, we used the same value of $\epsilon$ as our method for DQN. Other hyperparameters used for this section are listed in Table 2. Here, the learning rate for PPO was set smaller than those of other methods except the cases of task II. with $n_{\mathrm{agent}}=48$ and 64 to stabilize the learning. In addition, we let the minibatch size of DQN for task II. with $n_{\mathrm{agent}}=48$ and 64 larger than other tasks. The reason why we let these tasks be the exception is explained later.

The result is shown in Figs. 17 and 18. Here we plotted the mean value of all agents’ rewards for each algorithm, and all graphs are averaged over 8 independent trials. Seeing this figures, performance of our method on easy tasks, i.e. task. I with $n_{\mathrm{agent}}=12$ and task. II with $n_{\mathrm{agent}}=32$ resembles that of DQN, and slightly worse than the policy gradient methods. In addition, the learning speed of ours is slower than these methods. These points are thought to result from the $\epsilon$-greedy method, which decreases $\epsilon$ gradually and keeps this constant larger than the threshold value, $\epsilon_{\mathrm{min}}$. However, the performance on the task. I with $n_{\mathrm{agent}}=40$ is explicitly worse than all DRL algorithms. As we discussed in Sec. 4.1, the low information processing capacity of ESN is thought to be the cause of this phenomenon.

In task. II with $n_{\mathrm{agent}}=48$, the performance of DRL methods were worse than our method. It is because DRL methods failed to learn to make lanes and stagnated in some trials, whereas our method succeeded in learning in all trials. In this task, we investigated the performance of each DRL method for four values of the learning rate, $\eta=1\times 10^{-3},2.5\times 10^{-4}$, $5\times 10^{-5}$, and $1\times 10^{-5}$, because it was possible that choice of appropriate hyperparameter led to improve the learning. The result is shown in Fig. 19. As we can see from this figure, stagnation happened in almost all learning rates of every method. Note that the case of PPO with $\eta=1\times 10^{-3}$ succeeded in learning to make lanes in every trial exceptionally, however, $\eta$ of this case was so large that the learning became unstable. In Fig. 18 (b), we chose the learning rate $\eta=1\times 10^{-3}$ for DQN (the red curve of Fig. 19 (a)), and $\eta=2.5\times 10^{-4}$ for A2C and PPO (the green curves of Fig. 19 (b) and (c)), because these values seemed to perform better than other values in each method. We actually executed the similar calculation for DQN in the case that minibatch size is 64, the same size used for other tasks. However, all trials for each of four learning rate $\eta$ failed to learn to make lanes in this case. Hence we increased the minibatch size to improve the sampling efficiency. The similar investigation was also executed for the case of task. II with $n_{\mathrm{agent}}=64$, i.e. we calculated the performance of each DRL method for four learning rates for this case. However, in this case, only PPO with $\eta=1\times 10^{-3}$ succeeded or began to secceed in making lanes in some trials, while the other methods and other learning rates failed in all trials.

It is difficult to interpret the result of task. II because there are many unknown points on the algorithm combining ESN and LSPI method. One possibility is that agents choosing wrong actions cause the unclogging, as we discussed in Section 4.4. Indeed, ESN has lower information processing capacity and choose the wrong action more frequently than DRL methods, as we saw in task. I with $n_{\mathrm{agent}}=40$. This point is thought to led to the unclogging and final success in learning under $n_{\mathrm{agent}}=48$. In addition, PPO with $\eta=1\times 10^{-3}$, which succeeded in making lanes in all trials of $n_{\mathrm{agent}}=48$ and some trials of $n_{\mathrm{agent}}=64$, has unstable learning curves because of large learning rate, as we explained above. It is thought that such unstable behavior of agents also results in the unclogging. The other possible factor of the success of our method under $n_{\mathrm{agent}}=48$ is that sampling efficiency of the LSPI method is higher than the gradient descent method, because it utilizes all experiences gathered by agents completely by the linear regression. However, our method failed in making lanes under $n_{\mathrm{agent}}=64$, and the performance of this case was reversed by that of PPO with $\eta=1\times 10^{-3}$. Considering this point, effects of the above mentioned two factors on the final performance are thought to be complex. Hence, these effects should be studied also in other RL tasks we did not consider.

We also measured the computational time of whole simulation in task. I at $n_{\mathrm{agent}}=1$ and 12 for each method, and summarize the result in Table 3 (a). The calculation was executed by a laptop equipped with 12th Gen Intel(R) Core(TM) i9-12900H, and application software other than resident one was closed during the measurement of time. The real computational time for each simulation shown in Table 3 (a) includes the time spent on the calculation not related to the neural network, such as updating of the environment. Hence, to evaluate the computational cost of the calculation on the neural network including both training and forward propagation, we calculated the computational time of the case when agents choose their action randomly, and showed the difference of each simulation from this case in Table 3 (b). Here, in the random-action case, the neural network itself was not implemented. According to this table, computational cost of our method is less than half of that of each DRL method in both cases of $n_{\mathrm{agent}}=1$ and 12. However, as for the difference of time between these two cases, which corresponds to increase in computational time due to increase in agents, DQN with its minibatch size 64 is better than other methods and our method is second to this. This is because the time spent for the calculation of backpropagation for DQN under fixed minibatch size does not depend on $n_{\mathrm{agent}}$, whereas minibatch size for A2C, replay-memory size for PPO, and the size of matrices used in our method, $Y_{1},Y_{2}$, and $R$ given by eqs. (47) and (48), are all proportional to $n_{\mathrm{agent}}$. Considering the possibility that fixing minibatch size of DQN results in the worsening of sampling efficiency, we cannot say that DQN is more efficient algorithm than ours under a large number of agents, only from Table 3 (b). Indeed, in the case of difficult tasks such as task. II with $n_{\mathrm{agent}}=48$, we should increase the minibatch size for DQN at the cost of the computational time, as we explained above. Hence, the computational cost of our method under large $n_{\mathrm{agent}}$ is lower than those of DRL methods, whose sampling efficiency do not reduce as $n_{\mathrm{agent}}$ increases.