Research on the Calculation Model and Control Method of Initial Supporting Force for Temporary Support in the Underground Excavation Roadway of Coal Mine

Abstract

This paper proposes a temporary support system for improving the efficiency and safety of underground roadway excavation in coal mines. Firstly, this study establishes a calculation model for the initial supporting force of the excavation of roadway temporary support and a gray system-based automatic prediction model for the initial supporting force level, based on the mechanism of temporary support controlling the roof. These models enable the prediction of the required initial supporting force at different locations along the roadway’s temporary support area, thereby providing a basis for controlling the initial supporting force of the temporary support system. To achieve efficient and adaptive control of the initial supporting force of temporary supports at different locations, this study designs a support force controller based on Simulated Annealing Particle Swarm Optimization Proportional-Integral-Derivative (SAPSO-PID). This study establishes a mathematical model for the hydraulic cylinder pressure system controlled by the temporary support overflow valve and conducts a stability analysis and model verification. The study constructs a simulation control system for the initial supporting force based on SAPSO-PID using the combined simulation platform of AMESim and Matlab. The simulation results demonstrate that the proposed support force control system efficiently achieves adaptive control of the initial supporting force of temporary supports. An experimental system in the underground roadway of a coal mine is constructed to validate the results of the simulation analysis.

1. Introduction

The extraction of coal in underground mines presents a significant risk to the safety and well-being of miners [1]. Miners persistently encounter one of the highest level of occupational hazards among various industrialized nations globally. Achieving sustainable development in industrialized nations can be facilitated through the implementation of automation technologies in mining operations [2,3]. According to the International Energy Agency, global coal consumption has increased by 65% since 2000, with coal accounting for over 40% of all energy generation [4]. While coal is expected to continue its dominant role in the energy sector, the mining of deeper coal seams (below −800 m) has become increasingly challenging due to the depletion of shallow buried reserves [5]. The excavation of coal mine roadways is a critical and hazardous process in underground coal mines. The difficult excavation conditions and practical challenges faced during this operation severely hinder the development of intelligent and rapid roadway excavation techniques [6]. Therefore, achieving intelligent and efficient tunneling of roadways is imperative for ensuring high production rates and efficiency in coal mines, and is also an urgent research priority.

Restricted by the complex operating environment of roadways, the current excavation process in coal mines follows a step-by-step approach, as depicted in Figure 1. This process consists of six major steps: coal chute excavation, coal cutting and removal, transloading and transportation of coal, drilling machine operations, placement of anchor nets, and installation of permanent anchorage. However, the current excavation process is characterized by the use of numerous pieces of equipment that are independent of each other, a lack of coordination, fragmented workflow, and manual sequential operations. Particularly, the support process that follows excavation is time-consuming and labor-intensive, and the excavation and support operations cannot be conducted in a continuous manner as they rely heavily on manual intervention. Consequently, this situation hampers excavation efficiency and the rate of successful initiation [7,8,9].

This paper addresses the issues associated with the existing sequential excavation process by proposing the design of temporary supports for excavation lanes and establishing a designated area for temporary support, referred to as the advanced support area. This approach enables a new parallel operation process that encompasses excavation, temporary support, drilling, and anchoring. Figure 2 illustrates the complete set of equipment required for excavating, supporting, and anchoring, as well as the new parallel process. This set includes a roadheading machine, temporary support, transfer vehicle, drilling and anchoring machine, and other auxiliary components. The operational process of the parallel coal chute excavation, support, and anchoring is as follows: the roadheading machine, positioned ahead of the coal tunneling workface, performs slotting and section cutting, while the temporary support transfer vehicle moves the temporary support to the front of the roadheading. Subsequently, the temporary support transfer vehicle returns, allowing the drilling and anchoring machine to carry out permanent anchorage. By repeating these steps, the parallel and coordinated operation of excavation, temporary support, and anchoring is achieved, thereby enhancing excavation efficiency [10].

The temporary support structure of the tunneling roadway is depicted in Figure 3, comprising the roof beam, transverse expansion beam, guard plate, base, balance jack, and support column. The roof beam is connected to the transverse telescopic beam using a pushing cylinder, while the guard plate is hinged with the balancing jack. The support column, on the other hand, is hinged to the transverse telescopic beam via a ball vice. As a supporting element, the support roof beam directly contacts the roof plate of the roadway and actively bears the load of the surrounding rock. Likewise, the slide shoe base establishes direct contact with the roadway floor, allowing it to withstand the pressure exerted by the surrounding rock load transmitted through the support column and subsequently transferring that pressure to the roadway floor.

During temporary support of tunneling roadways, the roof support can be classified into three stages: initial bracing and increasing resistance, passive increasing resistance, and constant resistance overflowing. The initial bracing and increasing resistance stage involves active support of the roadway roof, with the support force during this stage referred to as the initial bracing force. The initial bracing force is a crucial parameter in temporary support, playing a significant role in roof control [11]. In tunneling roadway temporary support, the magnitude of the initial bracing force exerted by temporary brackets affects the effectiveness of roof plate support. Insufficient initial bracing force reduces support safety and increases the risk of early roof separation. Conversely, a high initial bracing force enables quick resistance buildup when pressure is exerted on the roof, minimizing top sinking and reducing the likelihood of roof collapse at the working face. However, excessive initial bracing force can lead to the coal tunnel roof crushing, causing secondary damage and escalating safety risks [10,12]. Therefore, it is particularly important for research on tunnel excavation to determine the appropriate initial support force for temporary supports and establish an efficient initial support force control system. Currently, research on temporary bracing for excavation in underground coal mines is still in its early stages. The research primarily focuses on optimizing the structural design of temporary bracing. Liu et al. [13] developed a coupled model of temporary support-roof plate system, suggested a temporary brace for roadways, and investigated its supportive performance. Ding et al. [14] examined the design of a self-moving temporary brace for roadways. Ma et al. [15] utilized the Optimal Space-Filling (OSF) method to optimize the structure of temporary support devices for coal mines, which plays a crucial role in the design of rectangular braces. Xue et al. [16,17] introduced a control method for the support force of temporary braces based on neural network PID. They obtained the input curve of surrounding rock pressure variation through simulation fitting using FLAC3D 6.0 software. Additionally, they developed a neural network PID control algorithm for adaptively tracking the control of surrounding rock pressure in the support force.

Several researchers have investigated the initial bracing force and regulatory methods for temporary bracing in coal mining faces, yielding useful research findings that can serve as references. Wang et al. [11] introduced the design principle of non-equal strength coupling support and temporary support, known as “low initial support, high work resistance,” which enhanced the adaptability of temporary support to high-stress and impact-prone roadways. Cheng et al. [18] developed an intelligent sensing system to monitor the state of braces and roofs. This system utilizes extensive monitoring data gathered from electro-hydraulically controlled hydraulic braces, enabling the analysis of the interaction between brace working resistance and surrounding rock and roof pressures. Ding et al. [19] devised a novel type of overrun support for general mining roadways, featuring high adaptability, support strength, and ample operating space. This design caters to the challenging conditions found in deep general mining working faces, characterized by severely unbalanced distribution of excavation, support, and anchorage duration. The authors further examined the stress–strain characteristics of the overrun support. Zhao et al. [20] presented a technique for the sectional adjustment of the hydraulic support setting load, aimed at pre-cracking and weakening coal walls in the context of comprehensive mining operations in hard and thick coal seams. Chen et al. [21] conducted a study on maintaining the spatial position of support braces during coal mining operations. However, their research did not encompass the calculation of initial bracing force or the investigation of regulatory systems.

A literature review revealed that no scholar has conducted a systematic study on the calculation and control of the initial support force of temporary support in underground coal mine roadways. In order to promote the safety and efficiency of temporary support, based on the above analysis, this study takes the control of the initial support force of roadway temporary support as the research object. It establishes a calculation model for the initial support force of underground coal mine roadway temporary support and a gray system automatic prediction model for the initial support force level, which achieves the demand prediction of the initial support force at different positions in the temporary support zone of the roadway. A controller for the initial support force, based on SAPSO-PID, was designed, along with a control system for the temporary support’s initial support force. The mathematical model of the control system for the hydraulic cylinder overflow valve of the temporary support was constructed and, subsequently, a model analysis and stability verification were conducted. A simulation control system for the initial support force based on SAPSO-PID was developed using the integrated simulation platform of AMESim and Matlab/Simulink, followed by a simulation analysis [22,23]. Furthermore, an experimental setup was developed in an underground coal mine to conduct the experiments. The simulation and experimental results clearly indicate that the proposed control system for the initial support force can effectively achieve efficient and adaptive control over the temporary support’s initial force in underground coal mine roadways. This research serves as a valuable reference for the calculation and control of the initial support force in temporary support for an underground mine roadway.

The subsequent sections of this paper are outlined below, providing a summary of the content examined and the methodologies employed in each chapter. In Section 2, a calculation model is established for the initial support force of temporary support in underground coal mine roadways. The range of support force levels is divided, and a gray system prediction model is developed for the automatic prediction of the initial support force level. In Section 3, a SAPSO-PID-based controller is designed for the initial support force, and a control system based on SAPSO-PID is established. In Section 4, the mathematical model of the overflow valve control system for temporary support hydraulic cylinder is constructed, and the model analysis and stability verification are conducted. A joint simulation platform using AMESim and Matlab/Simulink is built for the simulation analysis of the SAPSO-PID-based initial support force control system. The experimental analysis of the initial support force control system is conducted in underground coal mines in Section 5. Finally, corresponding conclusions are provided in Section 6.

2. Calculation Model for Initial Support Strength and Prediction Model for Supporting Strength

2.1. Calculation Model of Initial Support Force

Based on the different conditions during excavation, considering the stability of the working face and the coal wall in the roadway section, the initial bracing force of the brace must overcome the combined load of the direct top of the roadway and the self-weight of the brace roof beam (known as the equivalent direct top load). Additionally, once the brace reaches its rated working resistance, it can effectively handle the dynamic load caused by the basic top pressure of the control roadway roof, resulting in reduced roof sinking in the coal wall. It is essential that the initial bracing force and the rated working resistance maintain a reasonable proportional relationship (typically between 0.6 to 0.8). To minimize top plate sinkage and control the deformation of the coal wall under pressure, adjusting the relevant coefficient can increase the initial bracing force, enabling the brace to reach the rated working resistance more rapidly [10,12]. The setting force calculation model for considering the unstable state of the working face under direct roof and basic roof pressure is presented in Figure 4.

Assuming that the roof of the roadway fractures ahead of the coal wall, the equivalent direct top load should consider the portion of the coal wall that fractures in advance. Thus, the initial bracing force P₀ of the brace controlling the equivalent direct top must fulfill the following condition:

$$\left\{\begin{array}{c}{P}_0\ge R_1+R_2\\ R_1=l_s\left[\left(l_x+l_k\right)+\frac{1}{2}{h}_1\cot \alpha \right]h_1\gamma \\ R_2=l_{s}l{h}_2\gamma \end{array}\right.,$$

(1)

where R₁ is the lower direct top self-weight, kN; R₂ is the upper equivalent direct top self-weight, kN; h₁ is the thickness of the lower direct roof, m; h₂ is the upper equivalent direct roof thickness, m; h₃ is the basic roof key layer thickness, m; ∑h is the equivalent direct roof total thickness, m; l_s is the width of the support, m; l is the upper direct roof breaking length, m; l_k is the control roof distance of the support, m; l_x is the direct roof leading coal wall breaking distance; α is the direct roof breaking angle, °; and γ is the average bulk density of bedrock, kN/m³.

By rectifying the aforementioned Equation (1), we can determine the expression for the initial support force P₀ required to control the equivalent direct top in the unstable state as follows:

$$P_0\ge l_s\gamma \left\{\left[\left(l_x+l_k\right)+\frac{1}{2}{h}_1\cot \alpha \right]h_1+l{h}_2\right\}.$$

(2)

2.2. Classification of Initial Support Force Levels

In actual coal mine excavation sites, supports in various locations along the roadway frequently experience variations in cut-off disturbances and stress disruptions. As a result, the required magnitude of the initial support force varies between different supports. The pressure value in the lower cavity of the hydraulic cylinder of the brace support column serves as an indicator for determining the magnitude of the initial support force needed for the tunnel envelope and the upper control zone of the tunnel is classified into three distinct grading ranges using the grade zone division method. These ranges include the “stable support zone,” the “transition support zone,” and the “high resistance interference zone,” as demonstrated in Figure 5. The division of numerical regions in different ranges is derived from the supporting data of engineering sites under various conditions, combined with the experience of the research team. Such division can meet the requirements of most site conditions. The “stable support zone” is characterized by minimal cut-off disturbances whereas the “transition support zone” experiences frequent cut-off interferences and the “high resistance interference zone” is characterized by significant cut-off interferences and frequent disturbances requiring a high initial bracing force to prevent premature separation of the roof rock [10,11].

When characterizing the initial support force of the roadway using the pressure value in the lower chamber of the support column hydraulic cylinder, it is crucial to consider the duration of this pressure value, which is an important parameter that should not be disregarded. To mitigate the influence of fluid characteristics and other factors, this study determines the magnitude of the required initial support force for the roadway by integrating the average value of the pressure in the lower chamber of the hydraulic cylinder over time. The integration time selected for this purpose is 3 s.

Assuming that the bracket is positioned at the end of the roadway excavation during the support process, where the influence of cutting interference is minimal and the initial bracing force remains low (P₀ = 800 kN), the integration of time can yield a specific value for the pressure in the lower chamber of the hydraulic cylinder, denoted as P_t = 53 MPa·s. Similarly, the corresponding interval nodes for each level of initial bracing force for the bracket can be determined (1100, 1300, and 1600 kN) and this parameter can be used to establish

Table 1. Hydraulic cylinder chamber pressure as a function of the bracing force.

P0/(kN)	800	1100	1300	1600
Pt/(MPa·s)	53	73	86	110

as the reference standard for the bracket’s initial bracing force.

The adaptive fuzzy identification controller is designed based on the data provided in Table 1. The hydraulic cylinder’s lower cavity pressure, represented as P_t, is utilized as the input with a range of [53, 110] MPa·s. The initial bracing force rating, denoted as Q, is given as the output and has a value range of [1, 3]. Under these circumstances, a single-input single-output fuzzy identification controller (F₁) is established. A linear relationship exists between the initial bracing force and the pressure in the hydraulic cylinder’s lower cavity during the support process. Additionally, there is a proportional coefficient adjustment between the initial bracing force and the rated working resistance, which yields a numerical range indicating the level of the initial bracing force. Consequently, the Trapmf type affiliation function curve is selected, and its corresponding function is depicted in Figure 6.

The output variable for the initial bracing force level (Q) of the fuzzy controller is determined using a switch selector for the corresponding force range. The desired value of the initial bracing force, falling within the predefined range of 800–1600 kN, is illustrated in Figure 7.

2.3. Prediction of Brace Initial Support Strength Based on Grey System

Based on the completed range of initial support force levels described in the previous section, the positions of individual groups of temporary supports will shift as the excavation progresses in parallel. This necessitates an automatic real-time adjustment of the initial support force levels for each group of supports as they move dynamically.

The pressure monitored over time in the hydraulic cylinder cavity of the supporting column, which is related to the interaction between the support and the surrounding rock of the roadway, can be represented as a time series {X⁽⁰⁾(k) = [X⁽⁰⁾(1),X⁽⁰⁾(2),…,X⁽⁰⁾(n)], (k = 1,2,…,n)}. Due to the influence and constraints of multiple factors, it is challenging to consider the behavioral sequence of the system in the general mine roof envelope as having a straightforward exponential growth characteristic. Therefore, the gray system prediction model, specifically the DNGM (1,1) model, is applied using the direct estimation method of whitening differential equation parameters to predict and estimate the support strength within the controlled roof distance of the excavated tunnel [24].

Let X⁽⁰⁾ be a sequence defined as X⁽⁰⁾ = (x⁽⁰⁾(1), x⁽⁰⁾(2),···, x⁽⁰⁾(n), where x(0)(k) ≧ 0, k = 1,2,···,n. The sequence X⁽⁰⁾ is then cumulatively summed once to produce the sequence X⁽¹⁾, and the adjacent mean value generation sequence Z⁽¹⁾ of X⁽¹⁾ is given by:

$$\left\{\begin{array}{c}{x}^{\left(1\right)}\left(k\right)={{\displaystyle \sum }}_{i=1}^k{x}^{\left(0\right)}\left(k\right),k=1,2,\cdots ,n\\ z^{\left(1\right)}\left(k\right)=\frac{1}{2}\times \left(x^{\left(1\right)}\left(k\right)+x^{\left(1\right)}\left(k-1\right)\right),k=2,3,\cdots ,n\end{array}\right..$$

(3)

Then the whitening differential equation of the DNGM (1,1) model can be expressed as:

$$\frac{d{x}^{\left(1\right)}}{dt}+a{x}^{\left(1\right)}=bt+c,$$

(4)

where $\widehat{a}={[a,b,c]}^T$ are the parameters to be determined.

The time response of the DNGM (1,1) model can be obtained by the method of solving the Chi-square differential equation as:

$$x^{\left(1\right)}\left(t\right)=\left(x^{\left(1\right)}\left(1\right)-\frac{b}{a}+\frac{b}{a^2}-\frac{c}{a}\right)e^{-a\left(t-1\right)}+\frac{b}{a}t-\frac{b}{a^2}+\frac{c}{a}.$$

(5)

The final reduced equation of Equation (5) is:

$$\begin{array}{cc}\hfill {\hat{x}}^{\left(0\right)}\left(t\right)& ={\hat{x}}^{\left(1\right)}\left(t\right)-{\hat{x}}^{\left(1\right)}\left(t-1\right)\hfill \\ & =\left(1-e^a\right)\left(x^{\left(0\right)}\left(1\right)-\frac{b}{a}+\frac{b}{a^2}-\frac{c}{a}\right)e^{-a\left(t-1\right)}+\frac{b}{a}\hfill \end{array}.$$

(6)

From Equation (6), we have:

$$x^{\left(1\right)}\left(t+1\right)=e^{-a}{x}^{\left(1\right)}\left(t\right)+\frac{b}{a}\left(1-e^{-a}\right)t+\left(1-e^{-a}\right)\left(\frac{c}{a}-\frac{b}{a^2}\right)+\frac{b}{a}.$$

(7)

Order:

$$\alpha =e^{-a},\beta =\frac{b}{a}\left(1-e^{-a}\right),\gamma =\left(1-e^{-a}\right)\left(\frac{c}{a}-\frac{b}{a^2}\right)+\frac{b}{a}.$$

(8)

Equation (7) can be expressed as follows:

$$x^{\left(1\right)}\left(t+1\right)=\alpha x^{\left(1\right)}\left(t\right)+\beta t+\gamma .$$

(9)

Let $\widehat{\alpha },\widehat{\beta },\widehat{\gamma }$ be the parameter estimate of the above Equation (9), and replace x⁽¹⁾(t + 1) with the analog value ${\widehat{x}}^{\left(1\right)}\left(t+1\right)=\widehat{\alpha }{x}^{\left(1\right)}\left(t\right)+\widehat{\beta }t+\widehat{\gamma }$; the sum of squared errors can be obtained:

$$S={{\displaystyle \sum }}_{t=1}^{n-1}{[x^{\left(1\right)}\left(t+1\right)-\widehat{\alpha }{x}^{\left(1\right)}\left(t\right)-\widehat{\beta }t-\widehat{\gamma }]}^2.$$

(10)

According to the principle of least squares, the smallest S should satisfy:

$$\left\{\begin{array}{c}\frac{\partial S}{\partial \widehat{\alpha }}=-2{\displaystyle {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}}\left[x^{\left(1\right)}\left(t+1\right)-\widehat{\alpha }{x}^{\left(1\right)}\left(t\right)-\widehat{\beta }t-\widehat{\gamma }\right]x^{\left(1\right)}\left(t\right)=0\\ \frac{\partial S}{\partial \widehat{\beta }}=-2{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}\left[x^{\left(1\right)}\left(t+1\right)-\widehat{\alpha }{x}^{\left(1\right)}\left(t\right)-\widehat{\beta }t-\widehat{\gamma }\right]t=0\\ \frac{\partial S}{\partial \widehat{\gamma }}=-2{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}\left[x^{\left(1\right)}\left(t+1\right)-\widehat{\alpha }{x}^{\left(1\right)}\left(t\right)-\widehat{\beta }t-\widehat{\gamma }\right]=0\end{array}.\right.$$

(11)

Simplification:

$$\left\{\begin{array}{c}\widehat{\alpha }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}{(t)}^2+\widehat{\beta }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)+\widehat{\gamma }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)={\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}\left(t+1\right)x^{\left(1\right)}\left(t\right)\\ \widehat{\alpha }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)+\widehat{\beta }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{t}^2+\widehat{\gamma }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t={\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t+1\right)\\ \widehat{\alpha }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)+\widehat{\beta }{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t+\widehat{\gamma }\left(n-1\right)={\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)\end{array}.\right.$$

(12)

Solve the above system of non-homogeneous equations according to Cramer’s law such that:

$$B=\left|\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}{(t)}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{t}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t& n-1\end{array}\right|$$

(13)

$$B_1=\left|\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)x^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t+1\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{t}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t& n-1\end{array}\right|$$

(14)

$$B_2=\left|\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}{(t)}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)x^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t+1\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)& n-1\end{array}\right|$$

(15)

$$B_3=\left|\begin{array}{ccc}{\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}{(t)}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)x^{\left(1\right)}\left(t\right)\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{t}^2& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t{x}^{\left(1\right)}\left(t+1\right)\\ {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t\right)& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}t& {\displaystyle {\displaystyle \sum }_{t=1}^{n-1}}{x}^{\left(1\right)}\left(t+1\right)\end{array}\right|$$

(16)

the solution of the non-simultaneous system of equations can be obtained as:

$$\hat{\alpha }=\frac{B_1}{B},\hat{\beta }=\frac{B_2}{B},\hat{\gamma }=\frac{B_3}{B}.$$

(17)

Substituting the solution of the equation into Equation (8) yields the estimates of the parameters a, b, and c as:

$$\widehat{a}=-ln\widehat{a}, \widehat{b}=\frac{\widehat{a}\beta }{1-\widehat{\alpha }}, \widehat{c}=\frac{\widehat{a}\widehat{\gamma }-\widehat{b}}{1-\widehat{\alpha }}+\frac{\widehat{b}}{\widehat{a}}.$$

(18)

The estimated parameters are used to substitute into the reduced Equation (5), resulting in simulated and predicted values for the original series. Based on the algorithm and principles of the direct estimation of the whitening differential equation parameters mentioned above, a gray system prediction model can be established. The historical monitoring values of the pressure in the hydraulic cylinder cavity of the support column are used as input. The prediction results of the gray system are presented in

Table 2. Running results of gray system prediction model.

Number	Actual Data	Forecast Data	Residual	Relative Error
1	60.5	60.371	0.129	0.213%
2	60.1	60.2837	0.184	0.306%
3	60	60.042	0.042	0.070%
4	59.8	59.7454	−0.055	0.092%
5	59.9	59.3813	−0.519	0.866%
6	58.5	58.9344	0.434	0.741%
7	58.4	58.3859	−0.014	0.024%

. The average relative simulation error of the predicted data in the table compared to the actual measured data is 0.3302%, which is less than 0.01%. This indicates a high level of confidence in the reliability of the predictions. This system enables the automatic prediction of the initial support force for the temporary support.

3. Design of SAPSO-PID Initial Support Force Controller

After determining the initial bracing force value through the gray system approach, it becomes essential for the initial bracing force controller to accurately and adaptively track the target set value. Because of the complex and time-varying characteristics of the temporary bracing hydraulic control system, the traditional PID controller faces difficulties in achieving accurate tracking of the nonlinear systems. Particle Swarm Optimization (PSO), a popular optimization algorithm in the field of computation, offers an effective solution to this problem. The combination of PSO with the PID control algorithm allows for the dynamic optimization of PID control parameters, resulting in improved tracking error accuracy and control performance [25,26]. However, during the iteration and accumulation of search attempts, the PSO algorithm can easily become trapped in local optima. Hence, this paper presents the Simulated Annealing (SA) algorithm, known for its strong global search capabilities, as a means to optimize the PSO algorithm [27,28,29]. This section concentrates on the design and implementation of the SAPSO-PID initial bracing force controller.

3.1. PSO-PID Controller

PID control is widely employed as a fundamental control method due to its simplicity, safety, reliability, flexibility, efficiency, and ease of adjustment.

The PID control algorithm expression is:

$$G\left(s\right)=K_P+K_I\cdot \frac{1}{s}+K_D\cdot s,$$

(19)

where K_P is the proportional gain; K_I is the integral gain; and K_D is the differential gain.

In the Particle Swarm algorithm model, a search collective X_i = [x_i1,x_i2,···,x_iD]^T(i = 1,2,···,N) is formed by N particles in a D-dimensional target search space, and each particle’s own state is described by a set of position vectors p_i = [p_i1,p_i2,···,p_iD]^T(i = 1,2,···,N) and velocity vectors V_i = [v_i1,v_i2,···,v_iD]^T(i = 1,2,···,N) are described, in which the velocity of the particle directly affects the search distance of the particle at each step in the search space and can be adjusted according to the advantages and disadvantages of other particles and its own fitness. The optimal individual extremum searched by each particle is recorded as p_best = [p_bi1,p_bi2,···,p_biD]^T(i = 1,2,···,N), and the global optimal extremum searched by the whole particle population is recorded as g_best = [g₁,g₂,···,g_D]^T. Then, the particle properties are updated by Equation (20) as follows:

$$\left\{\begin{array}{l}{v}_{ij}\left(t+1\right)=\varpi \cdot v_{ij}\left(t\right)+c_1{r}_1\left[p_{bij}\left(t\right)-x_{ij}\left(t\right)\right]+c_2{r}_2\left[g_{bj}\left(t\right)-x_{ij}\left(t\right)\right]\hfill \\ x_{ij}\left(t+1\right)=x_{ij}\left(t\right)+v_{ij}\left(t+1\right)\hfill \\ \varpi ={\varpi }_{\max }-\frac{\left({\varpi }_{\max }-{\varpi }_{\min }\right)\cdot t}{T_{\max }}\hfill \end{array}\right.,$$

(20)

where i denotes the i-th particle; j denotes the j-th dimension of the particle; t denotes the current iteration number; v_ij(t) denotes the j-th dimensional flight velocity component of particle i when evolving to generation t; x_ij(t) denotes the j-th dimensional position component of particle i when evolving to generation t; p_bij(t) denotes the j-th dimensional individual optimal position p_best component of particle i when evolving to generation t; p_bij(t) denotes the j-th dimensional component of the optimal position g_best of the whole particle population at evolution to generation t; c₁,c₂ are acceleration factors; r₁,r₂ are random numbers of [0, 1]; T_max denotes the maximum evolutionary generation; $\varpi$_max denotes the maximum inertia weight; and $\varpi$_min denotes the minimum inertia weight.

The optimal design of PID controller parameters using the Particle Swarm algorithm follows a specific optimization process consisting of the following steps:

(1)

Generate a search particle swarm (either initialized or updated).

(2)

Assign particles in the swarm to the middle control parameters K_P, K_I, and K_D of the PID controller sequentially.

(3)

Link the Particle Swarm algorithm to the PID control model of the initial bracing force control system and calculate the performance indices associated with the control parameters.

(4)

Subsequently, the performance metric is used as the fitness value for each particle in the control algorithm, followed by a determination of whether it is possible to exit the algorithm by reaching the optimal value and breaking out of the loop.

3.2. Establishment of SAPSO Algorithm

The core of the simulated annealing algorithm lies in the thermodynamic process of solid cooling. Annealing refers to the process where the energy of an object reaches its minimum value as the temperature decreases. The Metropolis criterion is the central criterion in the simulated annealing algorithm, characterized by accepting not only superior solutions during the temperature decrease process but also inferior solutions based on the probability distribution of the temperature variable. This increases the likelihood of the algorithm escaping local optima. Therefore, combining the Particle Swarm algorithm with the Simulated Annealing algorithm can effectively reduce the probability of the Particle Swarm algorithm becoming trapped in local optima.

The Simulated Annealing algorithm requires setting an initial temperature T based on the initial state of the population during the initial iteration stage. The temperature linearly decreases by a certain degree at each iteration, forming an optimization search process that encompasses searching for new solutions and gradual cooling. The Metropolis criterion is employed to ascertain whether a disturbance-generated new solution replaces the global optimal solution. The expression for the Metropolis criterion is as follows:

$$P=\left\{\begin{array}{cc}1& f\left(p_{new}\right)<f\left(p_{old}\right)\\ \exp \left(-\frac{f\left(p_{new}\right)-f\left(p_{old}\right)}{T}\right)& f\left(p_{new}\right)\ge f\left(p_{old}\right)\end{array}\right..$$

(21)

Here, f(p_old) represents the fitness value of the particle in its initial position state, f(p_new) represents the fitness value of the particle at its new position, and P denotes the acceptance probability for the system transitioning from the initial state to the new position state.

The initial temperature is determined based on the characteristics of the initial particle, and it subsequently decreases by a cooling factor α after each iteration. The temperature is decreased according to the following function:

$$T\left(k+1\right)=\alpha \ast T\left(k\right),k=0,1,2\dots$$

(22)

When the temperature T is determined, the control algorithm establishes the conditions for the transition of particle state updates as follows:

(i)

If f(p_old) < f(p_new), indicating a decrease in energy after the state change, the particle accepts p_new as its current state.

(ii)

If f(p_old) < f(p_new), it denotes an increase in energy after the state change, indicating that the particle diverges from the optimal position. Subsequently, the control algorithm employs the rand (0, 1) function to generate a random number for comparison with the P value. If the P value is greater, the control algorithm accepts p_new as the current state. In contrast, if the P value is smaller, the current state remains as p_old.

The principles and characteristics of the Simulated Annealing algorithm and the Particle Swarm algorithm were used to fuse the two algorithms. The flow of the resulting algorithm, called the Simulated Annealing Particle Swarm algorithm, is presented in Figure 8.

By following the algorithm principles and steps outlined above, a Simulated Annealed Particle Swarm Optimized PID initial force controller was developed in Matlab. The advantage of Matlab 2020 software lies in its concise and comprehensible programming syntax, as well as its extensive library of encapsulated functions and algorithms, enabling convenient and efficient scientific computations and engineering analyses. The controller was incorporated into Simulink using S-functions of Matlab, providing the foundation for performance simulations and experimental applications.

3.3. The Composition of the Initial Bracing Force Control System

The initial bracing force control system based on SAPSO-PID consists of an input module, an initial bracing force calculation, prediction, and classification module, a valve-controlled hydraulic cylinder pressure control system module, an SAPSO-PID initial bracing force controller module, and a pressure feedback module. The input module includes all the initial data. The composition structure of the initial force control system is depicted in Figure 9.

4. Simulation Model Building and Result Analysis of Initial Support Force Control System

To assess the effectiveness and advancement of the initial bracing force controller proposed in this paper, a simulation control system for the initial bracing force was set up on the server. This allowed for preliminary simulation experiments to be conducted prior to performing downhole experiments.

4.1. Modeling and Stability Verification of Temporary Support Hydraulic Control System

In the context of simulation analysis, the initial step involves mathematically modeling the hydraulic control system of the temporary support. This mathematical modeling allows for a comprehensive understanding and accurate representation of the system’s behavior during simulations.

During the process of mathematically modeling the hydraulic pressure control system of the temporary bracket, it is permissible to make appropriate assumptions considering the internal and external parameters, as well as the motion characteristics of the temporary bracket. Simultaneously, certain functional components can be excluded, while focusing on the comprehensive modeling of significant components. This approach enables effective simplification of the system’s complexity.

In the hydraulic cylinder pressure control system structure, the simulation primarily focuses on analyzing the dynamic relationship among the proportional relief valve, hydraulic cylinder, and load. As a result, functional components such as the hydraulic check valve and proportional directional valve can be omitted during the mathematical modeling process. Additionally, the hydraulic cylinder that supports the column is responsible for both raising the top plate support and lowering the unloading drop frame. However, hydraulic cylinder cavity pressure control is only required during the rising action to achieve initial bracing force control. Consequently, the modeling process disregards the lowering action of the bracket. In addition, the hydraulic cylinder pressure control system considers both internal and external parameters, as well as system characteristics. Firstly, the connection pipe is assumed to be short, neglecting the effects of pressure loss in the hydraulic pipeline, fluid mass effect, and pipeline dynamic characteristics. Secondly, temporary assumptions are made regarding the internal and external leakage of the hydraulic cylinder and each valve, viewing them as laminar flow. The fluid temperature is set to a constant value.

The following will combine the above conditions to analyze the mathematical model of the pilot proportional relief valve, hydraulic cylinder, and feedback sensor in the hydraulic cylinder pressure control system.

4.1.1. Establishment of the System Model

(1)

Pilot Operated Proportional Relief Valve Modeling

The pilot-operated proportional relief valve comprises a proportional amplifier, a proportional solenoid, an electromagnetic force, a pilot valve spool displacement conversion link, and a pilot valve spool displacement to output pressure conversion link.

The role of the proportional amplifier is to process, compute, and amplify the input voltage signal into a current signal, which is then used to control the proportional solenoid. In essence, the proportional amplifier can be described as a link that establishes proportionality.

$$K_{pa}=\frac{I\left(s\right)}{U\left(s\right)},$$

(23)

where K_pa is the gain coefficient of the proportional amplifier, A/V; I(s) is the current of the proportional amplifier, A; and U(s) is the input voltage of the proportional amplifier, V.

The proportional solenoid plays a crucial role in generating thrust on the armature, causing the movement of the pilot valve cone and consequent alteration in the magnitude of the pilot valve opening pressure. Due to its high frequency response, the proportional solenoid can be likened to a proportional link.

$$K_{ps}=\frac{F\left(s\right)}{I\left(s\right)},$$

(24)

where K_ps is the current-force gain coefficient of the proportional electromagnet, N/A; F(s) is the corresponding output thrust of the electromagnet, N.

The mechanical balance equation of the pilot valve spool:

$$F-P_4{A}_4=m_{sp}\frac{d^{2}x}{d{t}^2}+B_s\frac{dx}{dt}+K_{es}x+F_{ss1}+F_{cf},$$

(25)

where P₄ is the outlet pressure of the pilot valve, N/m³; A₄ is the effective area of the pilot valve spool, m²; P₄A₄ is the algebraic sum of the hydraulic pressure on the effective area, N; m_sp is the mass of the pilot valve spool, kg; x is the displacement of the spool, m; B_s is the equivalent damping factor of the valve, N/(m·s⁻¹); K_es is the equivalent spring stiffness, N/m, which is the sum of the mechanical spring the sum of stiffness and steady-state hydrodynamic stiffness; F_ss1 is the steady-state hydrodynamic force, N; and F_cf is the Coulomb friction force on the spool, N. Where P₄A₄, F_ss1, and F_cf are usually negligible, the above equation is collated and the Laplace transformation is completed, the displacement of the pilot valve spool under the action of electromagnetic force can be obtained as:

$$x\left(s\right)=\frac{1}{\left(m_{sp}{s}^2+B_{s}s+K_{es}\right)}F\left(s\right).$$

(26)

The output pressure of the main valve relies on both the displacement of the pilot valve spool and the flow rate within the lower volume chamber of the main valve. By analyzing the force equation and flow continuity equation of the main valve, and applying suitable simplifications, it is possible to derive the expression that establishes the relationship between the displacement of the pilot valve spool and the output pressure, as indicated by Equation (27). The specific process of derivation will not be depicted in this context.

$$P_1=\frac{K_{r0}\left(1+\frac{s}{{\omega }_{rv}}\right)Q_4+K_{rs}x}{\frac{1}{{\omega }_{r0}^2}{s}^2+\frac{2{\xi }_{r0}}{{\omega }_{r0}}s+1},$$

(27)

where P₁ is the pressure at the lower end of the main valve, N/m²; x is the displacement of the spool, m; Q₄ is the flow rate of the lower volume chamber of the flow master valve, m³/s; ω_rv is the turning frequency of the main valve spool movement, Hz; and K_rs, K_r0, ξ_r0 and ω_r0 are coefficients related to the intrinsic frequency of the main valve of the flow valve, the flow coefficient of each structure of the valve, and the frequency associated with the valve motion.

By organizing the transfer functions of the aforementioned links in series using Equations (23), (24), (26), and (27), the block diagram of the transfer function for the pilot-operated proportional relief valve is illustrated in Figure 10.

The expression for the output pressure of the pilot-operated proportional relief valve is as follows:

$$P_1\left(s\right)=\left[\frac{K_{pa}{K}_{ps}{K}_{rs}}{m_{sp}{s}^2+B_{s}s+K_{es}}U\left(s\right)+K_{r0}\left(1+\frac{s}{{\omega }_{rs}}\right)Q_4\right]\frac{1}{\frac{1}{{\omega }_{r0}^2}{s}^2+\frac{2{\xi }_{r0}}{{\omega }_{r0}}s+1}.$$

(28)

(2)

Hydraulic Cylinder Modeling

Based on the motion characteristics of the hydraulic cylinder supporting the column, the mechanical equilibrium equation of the hydraulic cylinder can be derived as follows:

$$F_H=A_H{P}_H=m_p\frac{d^{2}z}{d{t}^2}+B_p\frac{dz}{dt}+K_{mp}z+F_{ef},$$

(29)

where F_H is the driving force generated by the hydraulic cylinder, N; A_H is the effective pressure cross-sectional area of the hydraulic cylinder, m²; PH is the cavity pressure of the hydraulic cylinder, N/m², because the relief valve is closer to the hydraulic cylinder, P_H can be approximated as equal to P₁; mp is the mass of the hydraulic cylinder plunger, kg; Bp is the equivalent damping factor of the hydraulic cylinder, N/(m·s⁻¹); K_mp is the equivalent spring stiffness, N/m; z is the displacement of the hydraulic cylinder plunger, m; and F_ef is the equivalent external load pressure, N, mainly including the gravity of the top beam of the support, etc.

Meanwhile, the flow equation of the hydraulic cylinder is:

$$Q_{Hc}=A_H\frac{dz}{dt}+C_{hc}{P}_H+\frac{V_t}{4{\beta }_{ec}}\frac{d{P}_H}{dt},$$

(30)

where Q_Hc is the hydraulic cylinder flow rate, m3/s; C_hc is the internal leakage coefficient of the hydraulic cylinder, L/min/bar; V_t is the working volume cavity of the hydraulic cylinder, m³; and β_ec is the effective volume modulus of elasticity of the hydraulic cylinder, MPa.

Combining Equations (29) with (30) and performing the Laplace transformation, we obtain:

$$Q_{Hc}=\left(\frac{A_H^{2}s}{m_p{s}^2+K_{mp}}+\frac{V_t}{4{\beta }_{ec}}+C_{hc}\right)P_H-\frac{A_{H}s}{m_p{s}^2+K_{mp}}{F}_{ef}.$$

(31)

(3)

Pressure Feedback Sensor Modeling

Under normal conditions, the pressure sensor can be simplified as a proportional element, and its transfer function can be expressed as follows:

$$K_f=\frac{U\left(s\right)}{A_H{P}_H\left(s\right)}=\frac{U\left(s\right)}{F\left(s\right)},$$

(32)

where K_f is the amplification factor of the pressure sensor, V/N.

By associating Equations (28), (31), and (32) and applying appropriate simplifications, we can derive the system transfer function relationship for controlling the hydraulic cylinder pressure using a relief valve. Figure 11 represents the relationship of the transfer function.

4.1.2. Control System Model Stability Analysis

To analyze the performance and validate the stability of the controlled object model within the system, the data for each system component was initially examined. Subsequently, the model parameters were selected considering the actual operating conditions. The stability analysis was conducted using the Nyquist criterion, while the system’s performance was evaluated through the Byrd diagram, unit step response, and unit impulse response.

Table 3. The relief valve control hydraulic cylinder pressure system parameters to take the value of the table.

Parameters	Unit	Value
Kpa	A/V	25
Kps	N/A	6.2
Kes	N/m	2910
msp	kg	0.002
A4	mm2	19.625
Kss	N/m	12,000
ωrv	Hz	3.4
AH	mm2	4.3
mp	kg	250
Bp	N/(m·s−1)	300
Kmp	N/m	3.2 × 106

presents the parameter values of the main components of a hydraulic cylinder system controlled by a specific model of safety valve. In practical applications, these parameter values may vary depending on the models of both the valve and the hydraulic cylinder.

From the relevant parameters listed in Table 3, the calculations yield the following values: ω_r0 = 18.54 Hz, ω_rv = 187.5 Hz, ω_rA = 0.244 Hz, and ξ_r0 = 1.832. These parameters are then substituted into the transfer function shown in Figure 11. Utilizing the Matlab programming environment, the system’s Nyquist diagram, Bode diagram, pole-zero diagram, and unit step response and unit impulse response curves are plotted, as depicted in Figure 12.

It can be seen from Figure 12 that the number of turns around the point (−1, j0) in the Nyquist diagram of the open-loop system is 0, and the number of poles and zeros in the right half of the pole-zero diagram is 0. According to the Nyquist stability criterion, it can be seen that the pressure control closed-loop system of the overflow valve-controlled hydraulic cylinder is stable; in the Bode diagram of the system, it can be observed that the phase margin is greater than 40° and the amplitude margin is greater than 20 dB, indicating the stability of the system. It can be seen from the jump response graph that the system stabilizes after 41.6 s, and it can be seen from the impulse response curve that the system will return to a stable state after 12.6 s.

Hence, the hydraulic cylinder pressure control system model, presented in this study and governed by a temporary direct overflow valve, exhibits stability. Additionally, it is evident that the inherent response speed and sensitivity of the hydraulic cylinder control system are sluggish and inadequate. To enhance both the response speed and sensitivity, an SAPSO-PID initial bracing force controller was introduced into the control system during simulation.

4.2. Simulation Model Building

To ensure a more realistic simulation, the control system for the simulation was constructed using the integrated simulation platform of AMESim and Matlab/Simulink [22,23]. The AMEsim 2020 software enables multi-domain system modeling and simulation, supporting the simulation of diverse systems encompassing mechanical, hydraulic, and control domains. It aids engineers in conducting system-level modeling, simulation analysis, and performance evaluation, contributing to the optimization of system design. The pressure system of the valve-controlled hydraulic cylinder was then modeled in AMESim, utilizing the mathematical model established in the preceding section. Figure 13 illustrates the depicted model.

The gray system, PID controller, and SAPSO-PID controller are implemented in Matlab/Simulink using programming techniques. Instead of the roadway roof load, the elastic load and damping load are applied, with a stiffness coefficient of 10,000 N/mm and a damping coefficient of 10 N·s/m. Based on the previous mathematical model, the gray system calculates and predicts the initial bracing force, which is set as the input signal of the simulation control system to achieve different levels of initial bracing force [30,31]. The joint simulation model of AMESim and Matlab/Simulink is illustrated in Figure 14.

4.3. Analysis of Simulation Results

The predicted initial bracing force levels from the gray prediction system for various time periods were utilized as input signals for the controller. Figure 15 and Figure 16 exhibit the dynamic response and error curves of the hydraulic cylinder’s initial bracing force control using both the PID controller and SAPSO-PID controller. In comparison to the PID controller, the SAPSO-PID initial force controller showcases improved response time by 85.6%, 81.4%, and 81.6% across the three initial bracing force levels, as well as reduced overshoot by 9.31%, 47.3%, and 46.7%, respectively. Moreover, the SAPSO-PID controller exhibits a superior tracking response performance with smooth amplitudes and minimal excitation, thus mitigating the chattering phenomenon caused by excessive input switching.

To further validate the dynamic response of the SAPSO-PID initial force controller in response to periodic signals, the system was subjected to periodic sine, sawtooth, and square wave input signals. The resulting dynamic response curves are presented in Figure 17.

In Figure 17a, when a sinusoidal signal was applied, the SAPSO-PID controller exhibited a reduction in system response time by 90.9% and a decrease in overshoot by 25.8% when compared to the conventional PID controller. In Figure 17b, when a sawtooth wave signal is used as input, the SAPSO-PID control algorithm achieves a remarkable reduction in system response time by 86.3% and an 87.3% decrease in overshoot, without accounting for the coarse error resulting from system delay. Similarly, in Figure 17c, the SAPSO-PID control algorithm exhibits a 60.98% reduction in the system’s response time and an 82.1% decrease in overshoot when a square wave signal is applied. These simulation results highlight the superior dynamic performance of the SAPSO-PID controller proposed in this study, particularly when operating under various working period signals. Please refer to

Table 4. Summary of simulation results.

Algorithm	Input Signal		Overshoot (kN)	Response Time (s)
PID	signal level	1	0.934	1.04
		2	0.167	0.59
		3	0.199	0.49
	Sine wave		0.597	1.98
	Sawtooth		0.314	4.02
	Square		0.848	0.41
SAPSO-PID	signal level	1	0.847	0.15
		2	0.088	0.11
		3	0.106	0.09
	Sine wave		0.443	0.18
	Sawtooth		0.040	0.55
	Square		0.152	0.16

for specific simulation data.

5. Experiment

5.1. Experimental Protocol

To validate the efficacy of the proposed temporary support initial bracing force controller, an experimental system was constructed in the excavation tunnel of Qianjuntai coal mine, which is under the administration of China Jingmei Group. Figure 18 displays the configuration of the experimental system. The temporary brace, as proposed by the researchers, was substituted with a joint temporary brace designed with a similar structure. The intelligent emulsion pumping station was employed as the power source for brace support. The intelligent emulsion pumping station is a vital equipment component within the hydraulic support system, primarily intended for emulsion supply to the hydraulic support. Emulsion, created by blending water and emulsifiers, functions in the hydraulic support system by providing lubrication, cooling, sealing, and corrosion prevention. The emulsion pump generates adequate pressure to transport the emulsion within the hydraulic support system. The pipeline system and valves regulate the flow and distribution of the emulsion, thus ensuring optimal support force and posture. In contrast to conventional manually adjusted emulsion pumping stations, the control system implemented in this study enables automated monitoring and adjustment of the pumping station’s operational status, responding to changes in roof support pressure.

Within the control box of the intelligent emulsion pumping station, the programmable controller utilizes the Controller Area Network (CAN) bus to establish communication protocol with the host computer and collection board. This allows for automated control of the valve port. The manipulation of the oil pipe’s output pressure, which is connected to the temporary support, facilitates regulation of the pressure in the lower chamber of the hydraulic cylinder of the column. As a result, the initial support force of the temporary support can be precisely adjusted.

The B&R Automation Studio 4.7 software was utilized for configuring, programming, monitoring, and diagnosing the temporary support initial bracing force control algorithm and parameters [32]. An analog signal, representing the range of initial bracing levels calculated based on geological conditions for each group of temporary supports, is input to the host computer to simulate the actual bracing requirements.

5.2. Experimental Results

Throughout the experiment, the pressure sensor measured the pressure in the lower chamber of the hydraulic cylinder of each group of temporary support columns, which can be considered an indication of the initial bracing force magnitude. Three groups were established with varying levels of initial bracing force. These levels were input to the host computer, which facilitated the adjustment and monitoring of the initial bracing force. The experimental findings are displayed in Figure 19.

Based on the experimental results, the SAPSO-PID initial bracing force controller demonstrated an improved average rise time compared to the conventional PID control across different initial bracing force levels. The SAPSO-PID controller achieved a reduction of 5.85%, 45.73%, and 14.83% in average rise time for the three different force levels, respectively. Moreover, the SAPSO-PID controller showcased significant improvements in system error fluctuation compared to the conventional PID control, exhibiting a reduction of 91.50%, 40.75%, and 71.97% in pressure error, respectively. These findings highlight the efficacy of the pressure feedback-based initial bracing force control system proposed in this paper, emphasizing its ability to achieve efficient adaptive control of the initial bracing force.

6. Conclusions

(1)

Developed a temporary support design for underground coal mine excavation tunnels and investigated the roof control mechanism of the support system. Created calculation and prediction models for the initial bracing force levels of the temporary support in the excavation tunnels, enabling accurate forecasting of bracing force requirements at different locations within the tunnels.

(2)

Designed an initial bracing force controller based on the SAPSO-PID approach and established a control system model specifically for the initial bracing force. Constructed a mathematical model to represent the pressure system of the pilot-operated relief valve-controlled hydraulic cylinder and subjected it to analysis and stability verification.

(3)

An integrated simulation platform combining AMESim and Matlab was employed to develop a simulation model for the initial bracing force control system utilizing the SAPSO-PID methodology. The simulation results demonstrate the outstanding performance of the SAPSO-PID initial bracing force control system in terms of response time and error control, highlighting its effectiveness in achieving efficient and adaptive control of the initial bracing force.

(4)

Carried out experiments within underground coal mine excavation tunnels where structurally similar temporary support systems based on the SAPSO-PID initial bracing force control were employed. The experimental findings demonstrate that the implemented bracing force control system effectively and adaptively manages the initial bracing force of the temporary support structures in the excavation tunnels. The temporary bracing control system utilized in the experiment demonstrates its universal applicability in coal mine excavation tunnels. It enables the adjustment of temporary bracing dimensions based on the tunnel cross-section dimensions in different projects while preserving the effectiveness of the control system remains intact.

The research presented in this paper contributes to the advancement of intelligent construction methods in the field of underground mining and the development of unmanned coal mining equipment technologies within the coal mining industry.