Reinforcement Learning for the Face Support Pressure of Tunnel Boring Machines

1. Introduction

Face stability is critical in shallow tunnels to avoid collapse. In mechanised tunnels, the face support is provided by the tunnel boring machine (TBM), e.g. slurry (SPB) or earth pressure balance (EPB) shields. An estimate of the support pressure is required for safe and efficient construction.

The problem of face stability can be solved with analytical, numerical and experimental approaches. The analytical methods are mainly based on the limit state analysis [1,2,3,4,5,6,7,8,9]. Face stability can also be investigated by numerical analysis [10,11,12,13,14,15,16,17]. Alternatively, the problem can be studied experimentally at 1g [13,14,18,19,20,21] or with centrifuge model tests [22,23,24].

Recently, machine learning has emerged as a promising technique for predictive assessment in geotechnical engineering, in general [25,26,27,28,29,30,31], and in tunnelling, in particular [32,33,34,35,36,37,38]. Some promising research domains for machine learning in tunnelling are the geological prognosis ahead of the face, the interpretation of monitoring results, automation and maintenance [32]. At present, however, research appears to be focussed on the following topics: prediction of TBM operational parameters [34,39,40,41,42,43,44,45,46], penetration rate [47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63], pore-water pressure [64], ground settlement [65,66,67], disc cutter replacement [68,69,70], jamming risk [71,72] and geological classification [73,74,75,76]). Few authors estimated the face support pressure of TBMs with machine learning [35,52].

The previous studies developed models that are trained after tunnel construction and thus pertain to the domain of supervised learning, the machine learning paradigm where the predictions are based on labelled datasets [77]. This study proposes a method to determine the face support pressure of TBMs with reinforcement learning. Within the realm of artificial intelligence, reinforcement learning is a generic framework for representing and solving control tasks in which the actions performed by the algorithm at one time are influenced by past decisions. Reinforcement learning algorithms are trained by incentivising them to accomplish a goal. Although reinforcement learning research is still in its infancy, there have been recently some exciting developments, most notably those showcased by the Google’s DeepMind research group [78,79,80,81,82], which demonstrated the ability to master Atari video games and popular board games such as chess and go. The transition of reinforcement learning from academia to industry is gaining momentum. Noteworthy examples are dynamic job shop scheduling [83], memory control optimisation [84,85], personalised web services [86,87], self-driving cars [88], algorithmic trading [89], natural language summaries [90] and healthcare applications such as dynamic treatment regimes and automated medical diagnosis [91]. Reinforcement learning industrial products include Google Translate, Apple’s Siri and Bing Voice Search [92].

Unfortunately, only few such studies are found in geotechnical engineering, especially in tunnelling. In particular, Erharter and Marcher presented a framework for the application of reinforcement learning to NATM tunnelling in rock [93] and to the TBM disc cutter replacement [68]. Zhang et al. [94] employed reinforcement learning to predict tunneling-induced ground response in real time.

In this study, the capability of the reinforcement learning algorithm to choose the best sequence of face support pressures is investigated by adapting the algorithm used by the DeepMind research group [80] and testing it on random geologies. The novelty of our method resides in the reinforcement learning approach, where the machine has no previous knowledge of the environment which it explores and is educated through the rewards defined by the user, as well as in the simulation of the environment with the finite difference method (FDM). This study shows that our model is capable of optimising the face support pressure, provided that a sufficient number of episodes are played.

In the following, the proposed method is presented (Section 2). The method is tested in environments of growing complexity, from analytical calculations to numerical analysis (Section 3). Its performance, limitations and possible improvements are discussed in Section 4. Finally, Section 5 concludes the paper.

2. Methods

In this section, the reinforcement learning algorithm is described. It is implemented with the interpreted high-level general-purpose programming language Python [95]. The algorithm is tested against analytical calculations (Section 2.1) and numerical analysis (Section 3.4).

One of three basic machine learning paradigms alongside supervised and unsupervised learning, reinforcement learning involves learning optimal actions that maximize a numerical reward [96]. In other words, reinforcement learning deals with the way (the policy) an intelligent agent (an algorithm) takes the actions that maximise a user-defined reward in a particular setting (the state of the environment).

The policy defines how the algorithm behaves in a certain situation. More precisely, it connects the states of the environment to the actions to be taken. As such, the policy is the core of a reinforcement learning agent, given that it determines its behaviour [96]. The well-established “epsilon greedy strategy”, one of the oldest policies [96], is selected in this study. An action is considered “greedy” if it is expected to return the maximum reward for a given state. However, since the environment is unknown a priori, an initial exploration of the environment is necessary to determine these actions. This exploration begins with the first TBM excavation where the face support pressure is randomly chosen at every round. The randomness decreases after every episode as the environment knowledge is exploited.

At each excavation round i, the agent chooses a random action $n_{\mathrm{a}}$ with probability $\epsilon$ and the action associated with the highest expected reward $Q_{\max }(S_{i+1},a)$ with probability $1-\epsilon$. $\epsilon$ is initialised at 1, which corresponds to a completely random selection, and decremented by $1/N$ after each episode, where N is the total number of excavations (the episodes). For N episodes, $\epsilon$ decreases by $1/N$ per episode until it reaches 0. In mathematical terms, let $r\in \mathbb{R}\cap (0,1)$ be a random number and $\epsilon =1-j/N$ at episode j, the agent takes the action $A_i$ at the state i according to equation (1), where $n_{\mathrm{a}}\in \mathbb{N}\cap (0,4)$ is a random integer.

$$A_{\mathrm{t}}\left\{\begin{array}{cc}{Q}_{\max }(S_{i+1},a)\hfill & \mathrm{if}r\ge \epsilon \hfill \\ n_{\mathrm{a}}\hfill & \mathrm{if}r<\epsilon \hfill \end{array}\right.$$

(1)

Hence, the random choice of the face support pressure is initially the dominant pattern that is slowly, but steadily, abandoned over time as the agent gains some experience of the environment in which it operates. In other words, the face support pressure $p_{\mathrm{f}}$ becomes more of a “conscious” choice based on the rewards collected and represented by the value function $Q(s,a)$, which returns the reward expected for the action a in the state s.

The rewards are user-defined and determined by the support pressure, the excavation rounds and the surface settlement. They reflect the definition of efficient construction: a safe process (consisting in the minimisation of surface settlement) executed with the least possible effort (determined by the lowest possible support pressure). The actual reward values are irrelevant, as long as the algorithm can maximise the expected reward, given the input data. The relative weight of the rewards, however, has an impact on the results [97].

A $+1$ reward is collected at each excavation round. The lower the face support pressure applied by the TBM, the lower the building effort. Hence, a reward corresponding to $-\frac{p_{\mathrm{f}}}{200\mathrm{kPa}}$ is assigned to every action. The surface settlement at every round causes a $-1$ reward for each mm of additional settlement. A $-100$ reward is assigned if the surface settlement is larger than 10 cm (in the analytical environment) or if the calculations diverges (in the finite difference environment). These outcomes terminate the episode (game over). A $+100$ reward is assigned at excavation completion. Each completed or terminated tunnel excavation defines an episode. The rewards are listed in Table 1.

The “playground” of the agent is named environment [98]. The actions are chosen and the rewards are received by the agent in the environment. In the present case, two classes of environments are created. First, an environment consisting of simple analytical calculations of the required face support pressure and expected settlements is considered (Section 2.1). The second environment is simulated numerically via the FDM as described in Section 3.4.

2.1. Analytical training environment

In the following, the initial training environment of the algorithm is described. For a tunnel with diameter $D=10$ m, a random 2000 m long geological profile is generated (Figure 1). The pseudorandom number generator is initialised with a constant seed so that the results can be replicated. The 2000 m length corresponds to the break-even point for the choice of mechanised over NATM tunnelling [99]. At first, this geology is kept constant and is thus the same for all the 100 training episodes. In the second instance, different random geological profiles are created at each episode (Section 3.2). The soil cover C is randomly initialised in the interval $(0.5,3)D$. The soil cover to diameter ratio is capped at $C/D=3$. At every 2 m (the excavation step or round length), a random slope is selected from the interval $(-1,+1)$, corresponding to the interval of slope angles $(-{45}^{\circ },+{45}^{\circ })$. The soil unit weight $\gamma$, friction angle $\phi$, cohesion c and Young’s modulus E are randomly initialised in the intervals shown in Table 2. The soil property values slightly change from their initial values at every excavation step (Table 1). The soil properties are re-initialised with a $1\%$ probability at every 2 m to simulate soil stratification.

By moving across states of these environments, the agent collects the rewards and stores them in the so-called value function $Q(s,a)$. The state representation is described in the next section.

2.2. State representation

The state is a representation of the information necessary and relevant to take an action. It is not the physical state of the environment, but rather a representation of the information for the algorithmic decision-making [100].

Since there are no rules to determine the state variables, domain knowledge, i.e. “the knowledge about the environment in which the data is processed” [101], must be used. According to the literature [4,7,12,16,22], face stability primarily depends on the soil unit weight, cohesion, friction angle and the depth of the overburden. Furthermore, soil settlement depends on the soil Young’s modulus E and the stress release [102]. Hence, the soil properties $\gamma$, c, $\phi$, E directly ahead of the tunnel face and the overburden C are normalised by dividing them by their maxima ${\gamma }_{\max }$, $c_{\max }$, ${\phi }_{\max }$, $E_{\max }$, $C_{\max }$ (Table 2) and selected as the state variables. The stress release is determined by the choice of the support pressure.

Generally, TBMs cannot estimate the soil properties ahead of the tunnel face [103]. However, boreholes are retrieved prior to soil excavation and analysed in the laboratory to obtain the material properties for the engineering design. Unfortunately, due to the soil heterogeneity, the property values cannot perfectly match reality, but are rather mean values that can be nonetheless used as a first approximation.

2.3. Face support pressure and settlement

TBMs provide face support pressures up to approximately 200 to 300 kPa in soft soils [104]. At every excavation step, the proposed model searches the optimal support pressure within the interval $(50,250)$ kPa. According to the guidelines [4], based on a limit equilibrium approach for drained conditions, the required support pressure is calculated with

$$p_{\mathrm{f}},\mathrm{req}=\frac{1}{A}\frac{(G+P_{\mathrm{V}})(sin\vartheta -cos\vartheta tan{\phi }^{\prime })-2T-c^{\prime }\frac{D^2}{sin\vartheta }}{sin\vartheta tan{\phi }^{\prime }+cos{\phi }^{\prime }}$$

(2)

where $A=\pi \frac{D^2}{4}$ is the cross-sectional area of the tunnel, G is the self weight of the sliding wedge, $P_{\mathrm{V}}$ is the vertical load from the soil prism, and T is the shear force on the vertical slip surface (Figure 2). The critical value of the sliding angle $\vartheta$ that maximises $p_{\mathrm{f},\mathrm{req}}$ is searched iteratively with the Python package scipy [105]. Note that the guidelines differentiate between $c_1^{\prime }$ and ${\phi }_1^{\prime }$ above the tunnel tunnel and $c_2^{\prime }$ and ${\phi }_2^{\prime }$ at the level of the tunnel. In this simulation, however, $c_1^{\prime }=c_2^{\prime }=c^{\prime }$ and ${\phi }_1^{\prime }={\phi }_2^{\prime }={\phi }^{\prime }$.

The proposed model estimates the support pressure $p_{\mathrm{f}}$ at every step to maximise the expected outcome according to the following criteria and compares it with the required pressure $p_{\mathrm{f},\mathrm{req}}$. In this simplified environment, the surface settlement occurs if $p_{\mathrm{f}}<p_{\mathrm{f},\mathrm{req}}$ and the stress release $\lambda$ is calculated as the ratio of the provided to the required tunnel face support pressure, so that $\lambda =1-\frac{p_{\mathrm{f}}}{p_{\mathrm{f},\mathrm{req}}}$. The corresponding soil settlement is then calculated as follows [102].

$$u=\frac{\lambda \gamma }{E}{\left(\frac{D}{2}\right)}^2$$

(3)

In reality, the settlement occurs not only due to the stress release at the tunnel face, but also due to other factors, such as the overcutting and ring gap [102]. Furthermore, an experience factor $K<1$ depending on ground stress and conditions, and tunnel geometry is generally considered in equation 3.

Since the model has no prior information about the required support pressure $p_{\mathrm{f},\mathrm{req}}$, the support pressure is chosen randomly during the first episode. Then, the model learns the optimal support pressure based on the rewards collected, as explained in the next section.

2.4. Deep Q-Network

In this section, the algorithm to maximise the rewards collected during the TBM excavation is elucidated. Since the material properties have continuous values, the number of possible states as described in Section 2.2 is infinite. Hence, the model cannot have a complete knowledge of the expected reward associated with the actions at each state and the knowledge of the value function $Q(s,a)$ is thus incomplete. This is where the Deep Q-Network (DQN) comes into play.

The DQN uses deep neural networks to approximate the value function $Q(s,a)$ [106]. This algorithm was developed by DeepMind, the Google research group, and was able to play six Atari games at record levels [78,80,81,82,107]. DQN is a particular approach to Q-learning, a method of learning optimal action values, whose main idea is to predict the value of a state-action pair, compare this prediction to the observed rewards and update the parameters of the algorithm to improve future predictions. Q-learning algorithms are formally described by equation (4).

$$Q(S_i,A_i)=Q(S_i,A_i)+{\lambda }_{\mathrm{r}}[R_{i+1}+\Gamma ·Q_{\max }(S_{i+1},a)-Q(S_i,A_i)]$$

(4)

In the present example, the DQN is trained to choose the best tunnel face support pressure based on the state variables. The DQN receives the state vector $S_i$ as the input and delivers the vector of the expected rewards associated with every action as the output. After the action $A_i$ is taken and the agent moves to the new state $S_{i+1}$, the reward $R_i$ is observed. The algorithm is run again with $S_{i+1}$ as the input and returns the action with the highest value $Q_{\max }(S_{i+1},a)$. Finally, the learning algorithm is updated to reflect the actual reward, using the mean squared error as loss function and minimising the difference between the predicted and target prediction of $Q(S_i,A_i)+{\lambda }_{\mathrm{r}}[R_{i+1}+\Gamma ·Q_{\max }(S_{i+1},A)-Q(S_i,A_i)]$.

In equation (4), ${\lambda }_{\mathrm{r}}$ and $\Gamma$ are the hyperparameters that influence the algorithm learning process. ${\lambda }_{\mathrm{r}}$ is the learning rate. The algorithm makes only small updates at each step with a low value of ${\lambda }_{\mathrm{r}}$ and large ones with a high value. The hyperparameter $\Gamma$ is the discount factor and controls how much our agent discounts future rewards when making a decision. Given that the discount factor is less than 1 (Table 3), future rewards are discounted more than immediate rewards.

The deep neural network is implemented using the PyTorch library [108] of the Facebook AI Research lab and its higher-level interface nn, as explained in [109]. Figure 3 shows the architecture of the neural network. The input layer is the state vector whose elements are ${\gamma }_i/{\gamma }_{\max }$, $c_i/c_{\max }$, ${\phi }_i/{\phi }_{\max }$, $E_i/E_{\max }$ and $C_i/C_{\max }$. The support pressure $p_{\mathrm{f}}$ is comprised in the discrete interval of natural integers $(50,250)$ kPa and is assigned in 50 kPa steps. Hence, the output layer is a vector with five elements, one for each possible support pressure (50, 100, 150, 200 and 250 kPa). The first and second hidden layers have each 50 neurons. This architecture has been previously applied to slope stability and tunnelling [36,110] and is the result of a trade-off between learning capabilities and computational effort.

However, DQNs are generally prone to training instability [78]. Hence, the following sections describe two common techniques used to stabilise the DQN, namely experience replay [111] and target memory [78].

2.4.1. Experience replay

Since the soil properties are variable in nature, the algorithm cannot simply memorize a sequence of support pressures $p_{\mathrm{f}}$. It must find the maximum support pressure that cause little to no surface settlement regardless of the geological setting.

However, reinforcement learning is essentially a trial-and-error process, which could lead to unstable solutions. This instability is particularly evident if learning is slow, which occurs when the rewards are sparse. This problem is deemed “catastrophic forgetting” [112] and is a common issue of gradient descent-based training methods in online training, i.e. when backpropagation is performed after each state. The essence of catastrophic forgetting is the push-pull between very similar state-action pairs that results in the inability to train the algorithm.

To this end, experience replay can be used to accelerate the learning process [111] by adding batch updating to the online learning scheme, as follows:

• In state i, the algorithm takes action a, and observes the new state $s_{i+1}$ and reward $r_{i+1}$,
• it stores this as a tuple $(s,a,s_{i+1},r_{i+1})$ in a list and
• continues to store each experience until the list is filled to a specific length, called “memory size” ($s_{\mathrm{m}}$).
• Once the experience replay memory is filled, a subset with a predefined batch size ($s_{\mathrm{b}}$) is randomly selected.
• The algorithm iterates through this subset and calculates the value updates for each subset; it stores these in the target array Y and the state s of each memory in X.
• Finally, it uses X and Y as a mini-batch for training and overwrites the old values in the experience replay memory when the array is full.

Thus, in addition to learning the value of the actions, a random sample of past experiences is also used for training. A further refinement is represented by the target memory.

2.4.2. Target memory

In DQNs, instabilities can arise by updating the parameters after each action. This instability is caused, inter alia, by the correlation between subsequent observations and can be overcome by updating the value function $Q(s,a)$ only after a prescribed number of episodes [78].

In the present problem, the rewards are “sparse”, i.e. they are significantly higher at the end of each episode than after every action. To solve this problem, the Q-network is duplicated. The parameters of this copy, i.e. the target $\widehat{\mathrm{Q}}$-network, are not up-to-date and lag behind the regular Q-network [78]. The $\widehat{\mathrm{Q}}$-network is implemented as follows:

• The Q-network is initialised with parameters ${\theta }_{\mathrm{Q}}$ and
• the $\widehat{\mathrm{Q}}$-network is a copy of the Q-network with distinct parameters ${\theta }_{\mathrm{T}}$. At first, ${\theta }_{\mathrm{T}}={\theta }_{\mathrm{Q}}$.
• The epsilon-greedy strategy is used with the Q-values of the Q-network to select the action a.
• The reward and new state $r_{i+1}$ and $s_{i+1}$ are observed.
• The $\widehat{Q}$-values of the $\widehat{\mathrm{Q}}$-network are set to $r_{i+1}$ at the end of the episode or to $r_{i+1}+{\lambda }_{\mathrm{r}}·{\widehat{Q}}_{\max }\left(S_{i+1}\right)$ otherwise.
• The $\widehat{Q}$-value is backpropagated through the Q-network (not the $\widehat{\mathrm{Q}}$-network).
• After a certain number of iterations, called synchronisation frequency (f), ${\theta }_{\mathrm{T}}$ is again set equal to ${\theta }_{\mathrm{Q}}$.

Up to this point, all the features of the algorithm 1 have been presented. In the next section, the results of the analytical environment are presented. The values of the hyperparameters ${\lambda }_{\mathrm{r}}$, $\Gamma$, f, $s_{\mathrm{m}}$ and $s_{\mathrm{b}}$ are chosen based on the sensitivity analysis of Section 3.1.

Algorithm 1:Workflow of the Reinforcement Learning algorithm

3. Results of the analytical environment

In this section, the results of the analytical environment are presented. The cumulative reward collected at every episode is shown in Figure 4a. The orange line represents the moving average of 10 episodes and the orange band is the range of $\pm 1$ standard deviations. Since the algorithm is trained on the geology of Figure 1 only, the reward increases steadily, barely oscillating around the moving average. In Figure 4b, the cumulative rewards obtained from the $1^{\mathrm{st}}$, ${25}^{\mathrm{th}}$, ${50}^{\mathrm{th}}$ and ${100}^{\mathrm{th}}$ episodes are highlighted. As expected, the reward increases more rapidly for the last episodes. Also, approximately at the “chainages” $x=200$ and 700 m, slight drops in the cumulative reward occur due to the abrupt changes in geology (Figure 1).

The support pressures and the resulting settlement at the episodes 1, 50 and 100 are shown at the left side and right sides of Figure 5, respectively. As the support pressure is randomly chosen at the first episode, Figure 5a exhibits a chaotic behaviour. As the agent learns the optimal support pressures, this behaviour is mitigated and the support pressure becomes more stable (Figure 5c). Figure 5c also proves that the algorithm is able to choose the support pressure in such a way that it almost encases the required pressure. Furthermore, even when $p_{\mathrm{f}}<p_{\mathrm{f},\mathrm{req}}$, the resulting settlement is negligible (Figure 5f). The model sensitivity to changes in the values of the hyparameters is studied in the next section.

3.1. Sensitivity analysis

Contrary to the parameters of the deep neural networks depicted in Figure 3 that are learnt by the algorithm, hyperparameters are set by the user. In Section 2, a number of hyperparameters have been introduced, such as the discount factor $\Gamma$, learning rate ${\lambda }_{\mathrm{r}}$, synchronisation frequency f, the memory $s_{\mathrm{m}}$ and batch $s_{\mathrm{b}}$ sizes.

The model sensitivity to these hyparameters is summarised in Table 3. Obviously, the results are not very sensitive to the discount factor and $\Gamma =0.15$ returns the highest cumulative reward. A low value of $\Gamma$ implies that the action taken at one state has little impact on the actions taken in the following states. This seems plausible for tunnelling where the settlement caused by an excessively low support pressure cannot be possibly offset by a later support pressure increase. The optimal learning rate is a matter of numerical stability. In this study, ${\lambda }_{\mathrm{r}}=0.001$ appears to be the optimal value and, if the learning rate is set at ${\lambda }_{\mathrm{r}}=0.01$, the cumulative reward cannot be properly maximised. The cumulative reward is not very sensitive to the synchronisation frequency where the optimal value is 5 episodes. The memory size has a more pronounced effect and the optimal value is 10. Finally, the batch size of 2 episodes provides the best results.

In summary, the proposed model is not very sensitive to the discount factor and the synchronisation frequency, it is starkly affected by the learning rate and the memory size in this particular problem. After the hyperparameters are optimised, the algorithm is tested against random geologies in the next section.

3.2. Random geologies

Prompted by the initial success in optimising the cumulative reward with constant geology, the algorithm is further tested against random geologies. Hence, the depth of the overburden and the material properties are no longer fixed as depicted in Figure 1, but change at every episode. Figure 6 exemplarily shows the profiles of the cohesion values at the episodes 2, 3, 4 and 5.

Clearly, since the environment changes at every episode, the cumulative reward shown in Figure 7a does not increase steadily as in Figure 4a, but swings markedly. The model is able to forecast the support pressure increases to accommodate for larger $p_{\mathrm{f}},\mathrm{req}$ in different settings (Figure 8). Furthermore, even when the support pressure provided fall below the required one, the resulting settlement is lower than about 4 mm.

The previous results are obtained with no additional model training. This means that the initial value of $\epsilon$ is set to zero. By starting from different values of $\epsilon$, it is shown in the following that no additional exploration is needed. Note that, strictly speaking, even if ${\epsilon }_0=0$ the Q-value is slightly updated as rewards are collected along the way. The agent, however, would perform the same action, given a certain state, until its expected rewards falls below that of other actions. The results of this analysis are shown in Table 4 where the mean cumulative reward and standard deviation are listed for different ${\epsilon }_0$. The highest mean cumulative reward and the lowest standard deviation are achieved for ${\epsilon }_0=0$. Therefore ${\epsilon }_0=0$ is selected as the optimum value, meaning that no additional exploration is required to optimise the algorithm performance in the random analytical environment. The numerical environment is described in the next section.

3.3. Effect of the Number of Episodes

When more episodes are played, the environment is explored more thoroughly. Hence, it is expected that the maximum cumulative reward increases with the number of episodes. Figure 9 shows that the maximum cumulative reward grows asymptotically up to approximately 683 after 400 episodes. It also shows that, the maximum cumulative reward obtained after 50 episodes already represents approximately 90% of the maximum reward achieved after 400 episodes.

3.4. Finite Difference Environment

Due to the simplified analytical formulations, the previous environment is not very realistic. It is, however, useful to demonstrate the capability of the model and to optimise its hyperparameters. Tunnels are often simulated with numerical analysis [113,114,115,116]. Hence, a more realistic finite difference environment is outlined in the following based on the finite-difference program FLAC^3D [117].

The mesh grid is 100 m wide, between 27 and 44 m high, and 100 m long. The tunnel diameter is 10 m and the distance between the tunnel axis and the bottom is 15 m. The soil surface slope changes at every 10 m of projected distance. The grid consists of 16,932 gridpoints and 15,300 elements with target dimensions of $2\times 2\times 1.5$ m (Figure 10). The linear elastic material law with Mohr-Coulomb failure criterion is considered. The soil properties are assigned within the intervals of Table 2 to randomly inclined layers. The property values are listed in Table 5. The displacements are fixed in the horizontal direction at the vertical boundaries and in the vertical direction at the bottom. The excavation is performed by removing the elements corresponding to the excavated soil. The support provided by the tunnel lining is simulated by fixing the displacement at the excavation boundaries. The support pressure is provided by applying a linearly increasing external pressure onto the tunnel face (Figure 11). This pressure is equal to 50, 100, 150, 200 or 250 kPa at the tunnel axis and increases linearly with depth according to the unit weight of the support medium ${\gamma }_{\mathrm{sm}}=12$ kN/m³. The support medium consists of excavated soil and additives for EPB machines or slurry suspension for SPB and mixshield machines [118]. The surface settlement is measured at the gridpoints situated on the surface every 10 m.

In the analytical environment, the surface settlement is calculated at every excavation round. In reality, the cumulative surface settlement depends on the previous excavation steps. Hence, the settlement increase is considered in the finite difference environment to account for the effect of the tunnel face support pressure at each step. The settlement reward for state t is thus expressed as the maximum settlement difference between the excavation steps $max(\Delta u_{t,j})$ of all the j-th measuring point along the soil surface.

Finally, the game over condition is not triggered by excessively large settlements, such as in the analytical environment, but by the divergence of the numerical solution.

Within the framework of transfer learning, the application of a machine learning model to a similar but unrelated problem [119], it is worth studying if and to what extent the algorithm used for the analytical environment can be deployed on the finite difference environment. The cumulative reward of $64.5$ is obtained in this environment if the support pressure is predicted with the algorithm trained in the analytical environment. This result is compared by retraining the model with ${\epsilon }_0=1$ in the finite difference environment. As it takes approximately 30 minutes to complete one episode, this training is computationally costly. Therefore, the model is trained only for 50 episodes or about 90% of the expected peak cumulative reward according to Figure 9.

The cumulative reward obtained at each episode is shown in Figure 12. For the first 20 episodes, the cumulative reward oscillates across approximately 65, corresponding to the reward obtained by using the model trained in the analytical environment. The cumulative reward oscillates across 75 between episodes 20 and 37 and stabilises at about 90 from episode 40.

Figure 13 shows the support pressures chosen along the chainage at each episode. This visualisation shows how the initial randomness begins to vanish at episode 20 where the first patterns start to emerge. In particular, the agent learns that the optimal support pressure between chainage 65 and 100 m is 250 kPa. It is also interesting to notice, that, between episodes 20 and 37, the agent preferably chooses a support pressure of 150 kPa between chainage 0 and 65 m. Starting from episode 37, it is evident that the optimal support pressure for the first 65 m is 100 kPa. This is consistent with Figure 12, where the cumulative reward starts increasing after episode 20 and increases again after episode 37.

Figure 14 shows the settlement after the last episode is completed. The surface settlement is up to 1 cm in the first 65 m of excavation and up to 2 cm in the last 35 m. Hence, it is evident that the actions chosen by the agent keep the settlement within reasonable limits.

From the previous analysis, it seems that the model developed in the analytical environment can be transferred to the finite difference environment provided that the model is retrained. These results are discussed in the next section.

4. Discussion

The results show that our model can optimise the support pressure by simultaneously control surface settlement within a reasonable threshold in both analytical and finite difference environments. Implementing model training in an analytical environment is relatively simple and a large number of episodes can be completed fairly fast. Moreover, this class of environments allows for hyperparameter tuning. On the other hand, reinforcement learning training in the finite difference environment (or, more generically, in numerical environments) is rather costly, see also [36]. Therefore, transfer learning is employed for hyperparameter tuning. As shown in the previous section, the model architecture and hyperparameters can be generally transferred to the finite difference environment, on the condition that retraining is performed starting from ${\epsilon }_0=1.0$.

Some limitations of this study can be highlighted and a strategy to amend them is outlined in the following. Firstly, albeit its use in engineering design, the finite difference environment cannot completely match reality. This is especially true in light of the simplifications considered in this study, such as the linear-elastic material law with Mohr-Coulomb failure criterion, the simulation of the tunnel lining as a zero displacement boundary condition, the absence of the ring gap and mortar and the deterministic soil properties values. These limitations can be overcome as follows:

• The adoption of more advanced constitutive models, the simulation of the lining with shell elements and the simulation of the ring gap and mortar [120,121]. It is perhaps worth noting that different types of segments (in terms of concrete class and reinforcement) and ring gap mortar pressures are chosen in practice. Hence, two additional agents could be implemented to predict the segment types and mortar pressures.
• The consideration of the spatial variability of soil properties with random fields, by varying the soil properties according to certain statistical distributions and correlation lengths [122]. Since random fields further complicate the environment, more advanced reinforcement learning algorithms might be adopted, such as the C51 [123]. Also, the definition of the state variables can be improved, e.g. by considering the soil properties in more than one point at each epoch.

The results match the expectation that the agent can be trained to predict the tunnel face support pressure. However, it is striking that the agent does not appear to need any additional training when deployed on random geologies (Section 3.2). This feature could be also theoretically tested with the finite difference environment. However, hyperparameter tuning with this environment is still computationally costly.

The results show that the DQN algorithm can be successfully used to control the tunnel support pressure and adapt to changes in the soil properties, such as variations in unit weight, cohesion, friction angle and Young’s modulus. One added value of the DQN in this context is that it can be used to develop more efficient and effective control strategies for maintaining tunnel face stability compared to traditional methods. The DQN has the capability to generalise to new situations, which can be useful in the case of changes in soil properties or overburden height. Furthermore, the DQN algorithm allows for an efficient use of the available data as it is not heavily dependent on its quality, which is a common problem with traditional methods.

5. Conclusions

In this study, the Deep Q-Network reinforcement learning algorithm is applied to control the tunnel face support pressure during excavation. The algorithm is tested against analytical as well as numerical environments. The analytical environment is used for hyperparameter tuning. The optimised model is used in the numerical environment.

It is found that:

• The algorithm is capable of predicting the tunnel face support pressure that ensures stability and minimise settlements among a prescribed range of pressures. The algorithm can adapt to geological (soil properties) or geometrical (overburden) changes.
• An analytical environment is used to optimise the algorithm. The optimal hyperparameters are found as $\Gamma =0.15$ (discount factor), ${\lambda }_{\mathrm{r}}={10}^{-3}$ (learning rate), $f=5$ (synchronisation frequency), $s_m=10$ (memory size) and $s_b=2$ (batch size). These hyperparameter values are effective also in the numerical environment.
• Although the algorithm is trained in a static environment with constant geology, it is also effective with random geological settings. In particular, it is found that using the algorithm trained with a constant geology can be used for random geologies without retraining.
• The maximum cumulative reward plateaus after 400 training episodes and about 90% of the peak performance is reached after 50 episodes.
• The algorithm proves effective both in the analytical and in the more realistic numerical environment. Training is more computationally costly in the numerical environment. However, the hyperparameter values optimised in the analytical environment can be efficiently adopted.

Future research can consider more refined environments (in terms of constitutive models, simulation of the lining, ring gap and mortar and random fields), provide more advanced state definitions (by considering the soil property values of various points) and use more refined reinforcement learning algorithms.

In spite of some limitations of this method, this study shows that the tunnel face support pressure can be estimated by an intelligent agent for design and possibly even during building operations.