The maximum Delaunay triangulation area method is a planting location determination approach based on Delaunay triangulation. In a Delaunay triangulation network, edge lengths represent the distances between neighboring trees, while nodes correspond to individual trees. Planting locations are determined by identifying the incenter of the largest triangle formed by the nodes, thereby reflecting both the presence of canopy gaps and the overall stand distribution pattern. To account for the influence of surrounding neighboring trees on the growth of replanted trees, the replanting foreground index (RFI) (Xuan et al., 2024) is incorporated to identify planting locations that have a greater impact on stand structure optimization. This method comprehensively considers various factors affecting the growth of replanted trees, thereby improving the accuracy and effectiveness of planting location determination and enhancing the efficiency of stand structure adjustment and optimization.

The calculation formula for the RFI is as follows (Equation 16):

RFI_i represents the replanting foreground index of the replanted tree i, DAA_i refers to the area of the Delaunay triangulation element in which tree i is located. Mc_i corresponds to the complete mingling of tree i, while U_i denotes its neighborhood comparison, and CI_i indicates its crown competition index. Moreover, δ_DAA , δ_Mc , δ_U and δ_CI represent the standard deviations of their respective structural parameters.

Relevant studies have shown that the optimal planting density range for Pinus yunnanensis is 1667–3333 trees per hectare. In this study, 3333 trees per hectare was selected as the upper limit for stand density after replanting. Based on the plot area, the maximum number of trees was determined for each plot: P1, P2, P3, P4, and P5 had upper limits of 1284, 1076, 418, 377, and 941 trees, respectively.

In mixed forests, the proportion and spatial arrangement of replanted tree species directly influence the degree of species mingling, while variations in DBH, H, and CW among different species substantially affect interspecific competition. In this study, replanting optimization simultaneously considered the spatial distribution and size structure of trees to achieve a more balanced species composition and stand structure. To promote the positive succession of secondary forests and enhance the mingling index, the main tree species widely distributed in the Cangshan region were selected. The seven chosen species—Pinus yunnanensis, Pinus armandii, Quercus acutissima, Vaccinium bracteatum, Camellia sinensis, Betula alnoides, and Ternstroemia gymnanthera—are either dominant or native tree species in the study area and represent the typical species composition and ecological characteristics of local secondary forests. Equal proportions were applied in the simulation to simplify the modeling framework and highlight species diversity and ecological complementarity. The replanted trees were initialized with an average DBH of 5 cm and an age of 5 years. Due to the limited availability of growth models for some species, models of ecologically similar species (Xi et al., 2015; Luo, 2021; Jiazheng et al., 2021) were adopted to estimate tree height, crown width, and crown length, and the results were subsequently validated using existing datasets. The configuration of replanted trees is shown in Table 2 .

Based on the data characteristics of different tree factors and the requirements of the prediction tasks, this study selected four different machine learning models to ensure predictive accuracy and reliability. The CNN-PSO (Ma et al., 2022) model is suitable for handling continuous variables with complex nonlinear relationships; its convolutional layers can effectively capture latent spatial patterns among features, and thus it was used to predict age and crown length. The Random Forest (RF) model (Xiaonan et al., 2024), with its strong anti-overfitting capability and ability to handle high-dimensional data, can robustly manage complex interactions among tree growth factors, making it suitable for predicting tree height and crown width. The Multilayer Perceptron (MLP) (Kim et al., 2025), due to its powerful function approximation capability, was employed to model the growth of diameter at breast height (DBH), a complex nonlinear regression problem. The Gradient Boosting Decision Tree (GBDT) (Nhat-Duc and Van-Duc, 2023), which excels in classification tasks—particularly with imbalanced data and combining multiple weak classifiers—was used to predict the probability of tree mortality, a binary classification problem. This targeted model selection strategy aims to maximize the accuracy and robustness of each prediction task.

Stand structure prediction involves forecasting parameters such as diameter at breast height, tree height, age, crown width, crown length, and mortality rate to understand the dynamic changes in stand structure. This enables dynamic optimization of stand structure and provides decision-making support for forest management, resource conservation, and ecological restoration.

We selected various indicators to construct the model, including age (AGE), diameter at breast height (DBH,1/DBH,DBH ²), tree height (H), crown width (CW), crown length (CL), number of trees per hectare (NT), stand density index (SDI), slope (SLO), aspect (ASP), height-to-diameter ratio (HDR), the sum of basal area larger than the target tree (BAL), and Hegyi competition index (HCI). To eliminate multicollinearity when building different predictive models, independent variables with a variance inflation factor (VIF) greater than 10 were excluded.

The diameter growth model was developed using the Stacked-MLP algorithm, with AGE, ASP, NT, SLO, SDI, BAL, and DBH ² selected as feature variables. The model’s target variable was set as the natural logarithm of diameter growth squared plus one, ln(DGI + 1). Since tree diameter growth occurs over long periods, shorter prediction intervals may fail to capture significant growth changes. Therefore, a five-year prediction cycle was used in the model.

The RF algorithm was used for tree height and crown width prediction models. The selected feature variables for tree height prediction were AGE, ASP, DBH, and SLOPE, with H as the target variable. For crown width prediction, the feature variables included AGE, 1/DBH, CL, ASP, SDI, HDR, and NT, with CW as the target variable.

The CNN-PSO was used for age and crown length prediction models. The feature variables for age prediction were 1/DBH, H, ASP, SLO, NT, and SDI, with AGE as the target variable. For crown length prediction, the selected feature variables were AGE, 1/DBH, ASP, SDI, HDR, and NT, with CL as the target variable.

Mortality prediction is a binary classification problem, and the GBDT model was chosen for this task. To improve prediction accuracy, the classification threshold was initially set at 0.5; however, this value is only applicable when the number of surviving and dead trees in the stand is approximately equal. In reality, mortality is a low-probability event, and cases where mortality and survival are equal are rare. Therefore, this study incorporated the HCI as an additional criterion. Trees with a mortality probability greater than 0.5 and experiencing high competition pressure (HCI > 0.75) were classified as dead. The selected feature variables included 1/DBH, H, BAL, and HCI, with tree mortality status (STATE) used as the target label.

Since the stand structure prediction model simultaneously predicts multiple stand factors, redundancy in feature variable selection is inevitable, which may lead to internal inconsistencies within the model. To address this, predictions were conducted sequentially in the following order: age, diameter at breast height, tree height, mortality, crown width, and crown length. This ensures that the feature variables used in all prediction models are updated accordingly during the prediction process.

The VIF values of the feature variables and the prediction accuracy of each model are presented in Tables 3 – 5 .

After optimization adjustments, the values of each sub-objective must not fall below their pre-optimization levels to ensure that the diversity of stand spatial structure does not decline. This means that the horizontal distribution pattern of the stand should become closer to a random distribution, the degree of species mingling should increase, the richness of vertical structure should be enhanced, and both size differentiation and competition pressure should be correspondingly reduced. During the optimization process, the number of tree diameter classes and species should not decrease, harvesting intensity should be controlled within 35%, and canopy density should not fall below 0.7. The number of replanted trees should remain within the stand density range, and after replanting, both the degree of mingling and the horizontal distribution pattern should be improved compared to prereplanting conditions (Equation 17).

M c 1 ¯ , S 1 ¯ , U 1 ¯ , C I 1 ¯ , W 1 ¯ , D 1 , T 1 , C d , N 1 respectively represent the values of complete mingling, stratification index, neighborhood comparison, crown competition index, uniform angle index, tree diameter classes, number of species, canopy density, and harvesting intensity after selective harvesting.   M c 0 ¯ , S 0 ¯ , U 0 ¯ , C I 0 ¯ , W 0 ¯ , D 0 , T 0 , C d , N 0 represent the values of the initial stand for the following indicators. M c 2 ¯ , W 2 ¯ represent the values of after complete mingling and uniform angle index after replanting. P D represents the stand density after selective redharvesting and replanting optimization.

Stand structure optimization is a multi-objective problem that typically involves multiple interrelated goals. Since these objectives are often mutually constraining, it is difficult to achieve their individual optima simultaneously. Therefore, it is necessary to adopt an integrated approach and balance multiple objectives from the perspective of the overall stand. In the optimization process, the complete mingling and stratification index should be increased to enhance stand diversity and stability; the crown competition index and neighborhood comparison should be reduced to alleviate competition among individuals and prevent excessive dominance by a few large trees; and the uniform angle index should be maintained within the range close to random distribution to ensure a reasonable horizontal pattern. Accordingly, in this study, five spatial structure parameters—uniform angle index, neighborhood comparison, complete mingling, stratification index, and crown competition index—were incorporated into a multi-objective optimization model using a multiplicative–divisive approach to construct the objective function (Equation 18):

W_i , Mc_i , S_i , U_i and CI_i represent the uniform angle index, complete mingling, stratification index, neighborhood comparison, and crown competition index of the reference tree, respectively. δ_W , δ_Mc, δ_S , δ_Wand δ_CI are the standard deviations of these structural indexes. N represents the total number of trees in the forest stand.

This study selected multi-agent deep Q-network (MADQN) as the solution algorithm. The MADQN model incorporates deep Q-network (DQN) and multi-agent Q-learning (MAQL). Each agent utilizes a deep neural network to approximate the Q-value function, enabling decision-making in large-scale and complex state spaces. By leveraging experience replay and target networks, the model achieves faster convergence, avoids overfitting, and ensures policy stability and generalization capability.

In the application of MADQN for stand structure optimization, selective harvesting and replanting serve as two key regulatory measures. Two agents, Agent1 and Agent2, were designed to achieve optimization. Each agent has its own tasks and objectives, interacting with the environment while also collaborating with each other. Through continuous learning, they work together to optimize the objective function of stand structure.Agent1 selects trees for selective harvesting based on the random selection method and adjusts its strategy according to the impact of tree removal on the objective function value. After harvesting, if the selected trees result in an increase in the objective function value, Agent1 receives a reward; otherwise, it is penalized. Through this process, Agent1 gradually learns how to select trees for harvesting, aiming to minimize unnecessary losses while maximizing the improvement of the stand’s objective function during the harvesting process.

Agent2 determines the number of replanted trees and their distribution based on the stand condition after Agent1 has completed selective harvesting and received its reward. Its goal is to compensate for gaps created by harvesting and introduce appropriate tree species to optimize stand structure, enhancing growth potential and ecological benefits. After harvesting, Agent2 utilizes the RFI to identify suitable planting locations. It then applies a curve trend-based approach (Xuan et al., 2024) to optimize the number of replanted trees. Specifically, it selects three evenly spaced replanting densities and evaluates their corresponding objective function values to establish a trend curve. Based on this trend, Agent2 adjusts the number of replanted trees by selecting new values with the same spacing before and after the current replanting density. This approach leverages curve monotonicity and extremum detection to assign rewards or penalties, refining the sampling strategy to gradually determine the optimal number of replanted trees. This ensures that under the given harvesting conditions, the replanting effect is maximized for optimal stand structure recovery.

Agent1 and Agent2 are interdependent and interact throughout the optimization process. Agent1’s selective harvesting decisions directly influence Agent2’s replanting strategy, while Agent2’s replanting results, in turn, affect the post-harvesting stand condition and consequently impact the objective function value. Through a reward and penalty mechanism, both agents continuously adjust their strategies, enabling a coordinated optimization process that dynamically balances selective harvesting and replanting to achieve the best possible stand structure.

The process of solving stand structure optimization using multi-agent deep reinforcement learning is illustrated in Figure 6 ( Algorithms 1 – 3 ).

During the subsequent optimization process, the selective harvesting and replanting agents continuously collaborate and interact, iteratively optimizing the current stand condition to achieve the most optimal stand structure. In each optimization cycle, the agents learn from environmental feedback, gradually approaching the structural characteristics of an ideal stand. However, real-world stand structures are influenced by multiple factors, and relying solely on selective harvesting and replanting may not directly achieve the desired stand structure. Therefore, a stand structure prediction model is incorporated as a crucial step to further enhance the accuracy and effectiveness of the optimization process.

The stand structure prediction model is used to forecast tree growth trends over a future period based on the current stand condition, as well as the potential changes after selective harvesting and replanting. After providing predictions on the current stand status, the agents further optimize the coordination between harvesting and replanting. At this stage, the selective harvesting and replanting agents not only rely on their original decision framework but also adjust the stand structure based on the prediction results. This process operates as a cyclic feedback mechanism, where each optimization generates a new stand condition, which in turn serves as the starting point for the next round of optimization. After each selective harvesting and replanting step, the updated stand condition is fed into the stand structure prediction model, and the model’s prediction results influence the next round of harvesting and replanting decisions. Through this approach, harvesting and replanting decisions are dynamically adjusted, ensuring that each step moves closer to an ideal stand structure. This iterative optimization process allows the model to identify the best regulatory strategies in a complex and dynamically changing stand environment. By avoiding a fixed strategy from the outset, the optimization process becomes more flexible and adaptive to stand dynamics.

The solution process for dynamic stand structure optimization is illustrated in Figure 7 ( Algorithm 1 ).

The parameter settings for the solution algorithm used in the experiment are shown in Table 6 . The initial iteration count and maximum iteration count for the solution algorithm were set to 0 and 1000, respectively. During the optimization process, the stand structure at different time periods is abstracted into a sequence. At the beginning of each iteration, the selective harvesting agent starts at the initial state ( s t a t e 1 = 0 ), while the replanting agent starts at the final state ( s t a t e 2 = 100 ). Both agents interact with the environment and collaborate with each other to decide whether to move forward ( s t a t e 1 = s t a t e 1 + 1 , s t a t e 2 = s t a t e 2 − 1 ) or move backward ( s t a t e 1 = s t a t e 1 − 1 , s t a t e 2 = s t a t e 2 + 1 ) . The iteration ends when the selective redharvesting agent and the replanting agent meet ( s t a t e 1 = s t a t e 2 ). To compare the performance of MADQN and MAQL, both algorithms were set with the same hyperparameters ( γ = 0.9 , l r = 0.01 , ϵ = 0.9 ). Additionally, for MADQN, the experience replay buffer size was set to 10000, with a batch size of 32. A three-layer fully connected network was used, with each hidden layer containing 24 neurons. These parameter settings were obtained through multiple experiments and fine-tuning to achieve optimal results.

In structured forest management, neighborhood comparison, representing size differentiation and competition intensity, complete mingling, indicating species segregation, and uniform angle index, describing horizontal distribution patterns, are the three most important spatial structure indexes. Considering the limitations of Pinus yunnanensis secondary forest plots, ( U ≤ 0.5 , M c ≥ 0.75 , 0.475 ≤ W ≤ 0.517 ) is selected as the ideal stand structure characteristic for dynamic stand structure optimization (Gangying et al., 2005, 2018).

As shown in Figure 8 , after implementing the two optimization schemes for coordinated selective harvesting and replanting, the stand structure indexes in each plot improved to varying degrees while meeting the constraint conditions, effectively enhancing stand structure. The average uniform angle index in each plot slightly decreased its deviation from 0.496, indicating that the horizontal distribution pattern of the stands became more randomly distributed. The complete mingling index significantly increased across all plots, particularly because the initial mingling degree in each plot was extremely low, leaving ample room for improvement. Notably, in Plot P4, the increase reached 16602.09%. Additionally, the crown competition index decreased substantially in all plots, indicating that tree competition pressure was alleviated after optimization. The stratification index showed a moderate increase, suggesting an improvement in vertical structural complexity and a more diverse vertical distribution pattern. However, the neighborhood comparison showed minimal changes across all plots. This is likely due to the fact that the initial average neighborhood comparison values were already in a moderate growth state, limiting the potential for significant improvement.

As shown in Table 7 , in the current stand structure optimization, both MADQN and MAQL significantly improved the objective function values. However, in terms of overall improvement, MADQN consistently outperformed MAQL. The objective function values for plots P1 to P5 under MADQN optimization increased from 0.3501, 0.3799, 0.3982, 0.3344, and 0.4294 to 0.5378, 0.5861, 0.5860, 0.5130, and 0.6034, respectively—higher than the values achieved by MAQL (0.5302, 0.5369, 0.5766, 0.5014, and 0.5906). The improvement rate under MADQN reached 49.40%, exceeding the 44.58% achieved by MAQL. These results indicate that MADQN is more effective in optimizing stand structure, guiding it toward a more optimal target state.

As shown in Figure 9 , MADQN outperformed MAQL in terms of the number of iterations required for optimization. After different numbers of iterations, MADQN exhibited a faster increase in objective function values, especially in Plots P2 and P4. This indicates that MADQN, which utilizes deep neural networks to approximate the Q-function, can find optimal strategies more quickly during the optimization process compared to MAQL, which relies on table-based Q-learning. As a result, MADQN requires less extensive exploration and achieves higher learning efficiency. It is worth noting that in Plots P1 and P3, although MADQN achieved a higher objective function value, MAQL required slightly fewer iterations to converge. This suggests that stand density and site conditions influence convergence performance, and in certain cases, MAQL’s convergence efficiency is not necessarily inferior to MADQN. As shown in Figure 10 , MADQN also demonstrated better overall performance in terms of computational efficiency. Across different plots, MADQN consistently maintained a lower or comparable runtime curve compared to MAQL, indicating superior time efficiency.

In the dynamic stand structure optimization of five plots using multi-agent deep reinforcement learning combined with structure prediction, most stand structure indexes showed significant improvements after optimization. The uniform angle index for all plots fell within the ideal range of [0.475, 0.517], indicating that the horizontal distribution pattern had reached a random distribution state. The complete mingling index increased substantially across all plots, shifting from a very low mingling state to a high mingling state. This demonstrates that dynamic optimization not only adjusts the spatial relationships among trees but also enhances stand stability and biodiversity at the species level. Notably, in Plot P4, the mingling index increased from 0.0028 to 0.7505, highlighting the strong adaptability of multi-agent deep reinforcement learning in adjusting tree species composition. Additionally, the crown competition index significantly decreased across all plots, indicating a considerable reduction in competition pressure. The stratification index also improved effectively, enhancing the vertical distribution pattern. In contrast, the neighborhood comparison showed minimal decline across all plots, suggesting that the pre-optimization stand already exhibited a relatively stable size differentiation state. As a result, despite some adjustments during optimization, fluctuations in this index remained small.

As shown in Table 8 , after incorporating structure prediction for dynamic optimization, the objective function values for all plots experienced significant improvements. Notably, in Plot P4, the objective function value increased from 0.3344 to 0.5863, achieving a remarkable 75.33% increase. Even in Plot P5, where the improvement was relatively smaller, the increase still reached 44.62%. These results indicate that dynamic optimization using multi-agent deep reinforcement learning combined with structure prediction effectively enhances stand structure stability and balance, making it more aligned with an ideal management state.