Constrained Multi-Objective Optimization of Simulated Tree Pruning with Heterogeneous Criteria

Virtual pruning of simulated fruit tree models is a useful functionality provided by software tools for computer-aided horticultural education and research. It also enables algorithmic pruning optimization with respect to a set of quantitative objectives, which is important for analytical purposes and potential applications in automated pruning. However, the existing studies in pruning optimization focus on a single type of objective, such as light distribution within the crown. In this paper, we propose the use of heterogeneous objectives for discrete multi-objective optimization of simulated tree pruning. In particular, the average light intake, crown shape, and tree balance are used to observe the emergence of different pruning patterns in the non-dominated solution sets. We also propose the use of independent constraint objectives as a new mechanism to confine overfitting of solutions to individual pruning criteria. Finally, we perform the comparison of NSGA-II, SPEA2, and MOEA/D-EAM on this task. The results demonstrate that SPEA2 and MOEA/D-EAM, which use external solution archives, can produce better sets of non-dominated solutions than NSGA-II.


Introduction
Pruning is one of the most important horticultural intervention techniques with which the vegetative and reproductive growth of fruit trees can be balanced. This balance is important for producing a stable and high quality product in consecutive yield seasons [1,2]. Proper tree pruning also helps to prevent tree degradation, and reduces the chance of disease by providing sufficient light and air flow through the tree canopy [3]. The knowledge of tree pruning is, thus, an essential skill for both professional and amateur fruit growers. In recent years, a number of software tools have emerged for computer-aided horticultural education [4][5][6][7]. They complement the field training by allowing the user to perform various actions on simulated tree models within a 3D virtual environment. The effects of such interactive tree manipulation on its growth are simulated to provide an informative feedback to the user. The use of software tools can be extended beyond education, since they allow analysis and comparison of various tree training techniques.
Another functionality enabled by software is algorithmic determination and evaluation of pruning performed on virtual trees with respect to specific goals. Existing research in this direction includes goal-oriented assessment of pruning effects [8], selection of pruning points in 3D tree reconstructions from images or point clouds as part of automated pruning [9][10][11], and algorithmic pruning optimization for support in computer-aided education [12]. In a previous study, the immediate and delayed effects of selective pruning on fruit tree development were assessed using the estimated amount of light received by the flower buds in the current and the next season [13]. In this case, the light intake was a proxy measure for the achieved ratio of vegetative and reproductive growth resulting from pruning, which was the object of optimization. It was shown that pre-growth and 1.
Integration of multiple heterogeneous objectives into a multi-objective tree pruning optimization framework.

2.
Introduction of secondary objectives, which are not the targets of optimization, but serve as constraints in the objective space of pruning solutions. 3.
Comparison and analysis of Pareto front approximations obtained with three popular advanced MO methods, namely NSGA-II, SPEA2, and MOEA/D-EAM.

Materials and Methods
In this section, we briefly review the tree growth simulation model, and describe the way pruning is formulated as an MO problem. We then define the new set of heterogeneous quantitative objectives for the pruning optimization problem, and describe the proposed constrained MO procedure.

EduAPPLE Tree Growth Simulator
Tree growth simulator EduAPPLE [7,15] was developed as a software tool that allows the user to observe the effects of various tree training techniques, such as pruning, weighing, and spreading, on the formation of an apple tree crown. To this end, EduAPPLE implements a parameterized development model that allows simulation of different growth behaviors of apple trees. EduAPPLE uses a simplified tree structure representation, where the tree skeleton is composed of short linear segments called metamers, which are linked to form branches. Metamers also present the basic unit of increment in growth simulation. Each Appl. Sci. 2021, 11, 10781 3 of 18 metamer consists of an internode and two buds, a terminal and a lateral one, which are used to extend or fork new branches upon growth. Internodes have assigned unique identifiers, which are used to indicate the cut locations for pruning. Figure 1 shows an example of a simple tree structure with depth-order internode labeling. Simulation of tree growth in EduAPPLE is based on a source-sink model, and performed in discrete seasonal steps. The yearly increments are determined by first calculating the growth resources accumulated by a tree, and then redistributing them to different parts of the tree. The growth is realized by replacing the shooting buds with new metamers. The main sources of growth material are the tree's own food reserves and photosynthesis. The amount of food reserve is modeled in EduAPPLE as a linear dependency on tree age A, while the quantity of photosynthetic product is estimated from total bud illuminance q i : Here, B denotes the set of all buds, while C 1 and C 2 are adjustable model parameters. The hyperbolic tangent non-linearity is used to limit the tree growth with age. Bud illuminance is computed by estimating the amount of shadow the bud receives from other tree elements. The details of this calculation can be found in the original EduAPPLE paper [15].
Following the allocation theory [28], the accumulated resources R are split into two exclusive portions, which are dedicated separately to vegetative and reproductive growth. The estimated reproductive requirements of a tree are calculated using the number N f of flower buds as [13]: Here, C 3 and C 4 are the allocation model parameters that can be adjusted to obtain specific behavior. The rest of the resources are allocated to vegetative growth. The reproductive pressure on resources can be reduced by pruning, which controls the balance of reproductive and vegetative growth across seasons.
Once determined, the vegetative growth resources are redistributed back to the buds in proportion to their light exposure, orientation, and distance from the roots. The probability of shooting is higher for well lit and vertically oriented buds, but reduced for older buds. Given the final amount r i of growth resources, bud i can sprout a sequence of r i new metamers. Upon creation, new buds are differentiated into vegetative and flowering ones, where the probability p f of a bud becoming a flowering one is another model parameter.

Pruning as an Optimization Problem
Algorithmic pruning optimization was introduced as an analytic tool by Strnad and Kohek [12]. In this context, a pruning is represented as a vector x = x 1 , x 2 , . . . , x n of cut locations, which correspond to individual internode identifiers. An example of this is shown in Figure 2, which illustrates the realization of a pruning vector for a given tree model. Figure 2 also shows that a cut can have no effect if another cut is present higher up in the tree hierarchy. This behavior is actually helpful during optimization, where changes to the pruning vector can activate and deactivate individual cuts. Realization of a pruning vector (i.e., genotype) x = 3, 5, 8 is a pruned tree model (i.e., phenotype). In this example, the cut on internode 5 is deactivated by the cut on internode 3 up the tree hierarchy.
The task of pruning optimization is to find a pruning solution that maximizes the values of selected quantitative objectives. Maximization of one objective, however, may result in a reduction of another objective. For instance, optimizing the crown shape towards some desired form may reduce the total resource accumulation of the tree, which affects the balance of vegetative and reproductive growth. Such conflict of goals is at the heart of multi-objective optimization proposed in this paper.
Let us formalize the pruning optimization problem with respect to a given set of objective functions f = ( f 1 , . . . , f m ). The objective value of a pruning solution vector x is, in this case, given by f (x) = ( f 1 (x), . . . , f m (x)). A solution vector x i is said to be dominated by another solution vector x j , which is denoted by x j x i , if the following conditions hold: In other words, pruning x j dominates pruning x i if it is not worse in any objective, and is strictly better in at least one. Note that the above formulation of dominance corresponds to the case of objective maximization used in this paper. For the minimization case, one needs to invert the inequalities.
The goal of MO is to construct a set X = (x 1 , x 2 , . . . , x q ) of non-dominated pruning solutions, known as the Pareto front. Only an approximation of the Pareto front can usually be obtained in practical problems. For this task, one can choose between several MO algorithms [29]. In this paper, the adaptations of three popular and efficient evolutionary MO methods were employed, known as NSGA-II [25], SPEA2 [26], and MOEA/D-EAM [27]. The latter is an extension of the original MOEA/D [30] method with a modified selection scheme that uses an external archive. In the terminology of evolutionary algorithms, the pruning solution vectors correspond to the genotypes, and their realizations on a tree model correspond to phenotypes. During the pruning optimization, the search is conducted in the genotype space, while the fitness evaluation is performed in the phenotype space.

Objectives
In a previous study by Strnad et al. [13], the bi-objective value of a pruning solution x for a given tree model T was defined using the pre-growth objective value f pre and the post-growth objective value f post : Appl. Sci. 2021, 11, 10781

of 18
The quantity f (T /x ) is the estimated light intake of a tree T after being pruned according to x. The pre-growth objective f pre reflects the short-term (i.e., immediate) effects of pruning, which can be assessed directly on the pruned tree model. On the other hand, the post-growth objective f post serves as a proxy measure for long-term (i.e., delayed) pruning effects. These are estimated by averaging the light intake over the results of multiple stochastic growth simulations of the pruned tree.
Improving the light conditions within the crown is one of the principal goals of tree pruning. However, there are other important aspects that practitioners need to consider, especially when performing the pruning as a corrective measure in neglected or damaged trees. For example, the amount of removed biomass may need to be constrained in order to prevent inflicting too much stress on the tree, or the disrupted tree balance needs to be restored by improving biomass distribution. Such heterogeneous criteria can constitute conflicting objectives for the pruning optimization task. In this paper, we analyze the results of such multi-objective optimization by introducing the following pruning objectives: • The average light exposure of flower buds after pruning, denoted by f 1 , and calculated as: Here, B f (T /x ) ⊂ B(T /x ) denotes the set of the pruned tree's flower buds. Note that this objective promotes aggressive pruning, since the maximum value can be achieved by a few fully exposed buds. • Conformance of crown shape to the desired training form, which is denoted by f 2 . Its value is estimated using the inverse Hausdorff distance between the convex hull of the pruned tree and the hull's bounding volume of target symmetric shape: The convex hull H of the pruned tree is constructed using its branch tips. Once the hull is obtained, its bounding volume C with target shape is determined ( Figure 3). In our experiments, a cylindrical target shape was pursued, but other common training forms can be used (e.g., conical). The Hausdorff distance between two shapes is defined as: where d denotes the standard Euclidean distance between two points. • Tree balance, denoted by f 3 , and represented by the inverse horizontal distance d c between the above ground biomass center of gravity c and the vertical axis going through the stem origin: The center of gravity location c for a tree T /x with internodes I(T /x ) is computed as: Here, s i is the midpoint of internode i, and m i is its mass. The computation of the latter is simplified by assuming homogeneous wood density, so the mass of internode i with radius r i and length l i is proportional to its cylindrical volume: Appl. Sci. 2021, 11, 10781 6 of 18 • The proportion of remaining tree biomass, denoted by f 4 and calculated as: Here, m(T ) is the above ground biomass of tree T : Figure 3. Convex hull H, its bounding cylinder C, and the center of gravity c (yellow sphere) for a tree.

Optimization Method
In this paper, pruning optimization is treated as a discrete combinatorial optimization problem, where the task is to find the set X of non-dominated pruning solution vectors with respect to a set of objectives f = ( f 1 , f 2 , f 3 ). Objective f 4 , defined in Section 2.3, is itself not an interesting target of optimization, since it can be maximized trivially by no pruning at all. Therefore, we propose a special role for objective f 4 in order to constrain the search to certain regions of the solution space. For instance, we may prescribe that no more than 10% of the tree's biomass should be removed, which could be imposed by the constraint f 4 > 0.9. This is especially important in order to regulate aggressive pruning stimulated by objective f 1 . Adherence to constraint objective bounds is implemented in a soft manner by marking the violating solutions as dominated. Such approach does not completely reject borderline solutions, making it possible to use them for the derivation of better successors. We use only f 4 as a constraint objective in this study, but other metrics could be used. The distinctive property of a constraint objective is that it evaluates the pruning solution globally, and not with respect to individual cuts. In general, both the lower and the upper bounds for a constraint objective can be specified by the user.
The use of constraint objectives to specify feasible regions in an auxiliary objective space is a new addition to the multi-objective pruning optimization framework, which complements heuristic constraints implemented in a decision space. The latter have been introduced by Strnad et al. [13] in order to reduce the size of the solution space and make the search more efficient. This is achieved by two mechanisms: • By limiting the cut locations to those with certain local properties. In our tests, the following restrictions were enforced: -Cutting was limited to the wood at most A max years old, in order to prevent pruning of scaffold limbs.

-
Pruning was restricted to locations following a branch fork that result in removal of at least m min metamers.
• By constraining the dimension of the solution vectors (i.e., the number of cuts) to a range [d min , d max ].
The high-level concept of the proposed multi-objective pruning optimization is outlined in Figure 4 and Algorithm 1. It is an adaptation of a general evolutionary process Appl. Sci. 2021, 11, 10781 7 of 18 for the construction of Pareto front approximation. Given a tree model T as input, a list I cut (T ) ⊂ I(T ) of potential cut locations is first built according to the heuristic constraints described above. From this list, the initial population  for k = 0 . . . M/Ps do 6: for ∀x 10: P (k+1) ← apply evolutionary operators on P (k) using Cr, Mr and P 13: return X The procedure then enters the optimization loop. The number of iterations N is either specified directly, or calculated from the maximum number of objective evaluations M as N = M Ps . Within the loop, three main steps are performed: 1.
The current population P The set of non-dominated solutions is updated, using the constraint objective f 4 to demote unsuitable solutions. 3.
The next generation of solutions is produced by applying evolutionary operators to the current population.
Evaluation of individual candidate solution x i is performed by pruning the tree model according to x i , and computing the values of objectives on the resulting pruned tree. This part of the computation is performed by the simulation model, and is independent of the used optimization algorithm.
The set X of non-dominated pruning solutions is updated next, which is a step that depends on the evolutionary method used. In NSGA-II, X is updated with non-dominated solutions from the current population P (k) . In SPEA2 and MOEA/D-EAM, the current archive P (k) is used to update X in each iteration. In both cases, solutions with objective value f 4 below the prescribed limit are not included in X even if non-dominated.
The last step is the derivation of a new population P (k+1) from the current population P (k) using a method-dependent procedure. In all of our experiments, the following problem-specific implementations of evolutionary operators are used: • Selection of two parent vectors from the population (NSGA-II) or the mating pool (SPEA2 and MOEA/D-EAM). • Uniform crossover of parent vectors with probability Cr. • Mutation of child vectors with change probabilities P = P move , P add , P remove and mutation rate Mr.
The selection step is method-dependent. In NSGA-II [25] and SPEA2 [26], binary tournaments are used to select the parents from the population and the mating pool, respectively. The mating pool is also used by MOEA/D-EAM [27], but its selection scheme combines random selection with a neighborhood-based one. During crossover, special consideration is necessary when recombining parent vectors of different sizes. In such cases, the surplus genes of the longer parent vector are distributed uniformly between the two child vectors ( Figure 5). Mutation of a child vector is performed in two steps. First, a randomly selected cut is relocated to a different valid position from I cut (T ) with probability P move , removed from pruning with probability P remove , or added to pruning with probability P add . Cut addition and removal are possible only if the size of the resulting pruning vector remains within the range [d min , d max ]. In the second phase, all of the other cuts are replaced by randomly selected alternatives with small probability Mr. The complete implementation of crossover and mutation is shown in Algorithm 2. The exploration/exploitation behavior of the search is regulated through meta-parameters Cr, Mr, and P. Sufficient exploration can usually be achieved even with small values of crossover and mutation rates, because local changes of a genotype (i.e., pruning vector) can result in significantly different phenotypes (i.e., pruned tree models).

Results
The goals of the experiments presented in this section are to: 1.
Determine the properties of Pareto front approximations obtained from pruning optimization with multiple heterogeneous objectives. In particular, we are interested in the relation of generated pruning solutions to three reference proposals, corresponding to non-pruning, non-selective pruning to cylindrical shape, and distance-based pruning, where the secondary branches are removed if their distance to the primary ones is below the threshold.

2.
Investigate the pruning patterns produced by MO, and how they reflect the trade-offs between conflicting pruning goals.

3.
Evaluate the effect of the proposed constraint objectives on the resulting Pareto front approximations and solution characteristics.

4.
Compare the performance of NSGA-II, SPEA2, and MOEA/D-EAM on this problem, which is important from the perspective of future framework development.
In the continuation, we first describe the configuration and workflow of experiments, and then address the above research questions in order.
The hardware used for the experiments was a desktop computer with Intel i7 CPU, NVIDIA GeForce GTX 1060 GPU, and 16 GB of RAM. The software environment included the Linux operating system (kernel version 5.13.11) and the GCC compiler version 11.1.
The tree models for experiments were generated by using the simulator with growth parameters C 1 = 80, C 2 = 40, C 3 = 80, C 4 = 20, and p f = 0.08. The models are shown in Figure 6, and present different levels of complexity for pruning optimization. Their main structural properties and the corresponding combinatorial search space sizes are reported in Table 1.
In the continuation, we first describe the configuration and workflow of experiments, and then address the above research questions in order.
The hardware used for the experiments was a desktop computer with Intel i7 CPU, Nvidia Geforce GTX 1060 GPU, and 16 GB of RAM. The software environment included the Linux operating system (kernel version 5.13.11) and the gcc compiler version 11.1.
The tree models for experiments were generated by using the simulator with growth parameters C 1 = 80, C 2 = 40, C 3 = 80, C 4 = 20, and p f = 0.08. The models are shown in Figure 6, and present different levels of complexity for pruning optimization. Their main structural properties and the corresponding combinatorial search space sizes are reported in Table 1. Figure 6. Tree models used in the experiments.  Figure 6b 1,234 145 9.75 × 10 27 Figure 6c 1,919 253 2.54 × 10 34 Figure 6d 2,506 341 5.92 × 10 37 Figure 6e 2,647 354 1.55 × 10 38 Figure 6f 3,740 429 2.18 × 10 40 Tuning of methods' meta-parameters was performed with a grid search using the pools Ps ∈ {20, 30, 50, 100}, Mr ∈ {0.01, 0.03, 0.05, 1.0}, and Cr ∈ {0.7, 0.8, 0.9, 1.0}. The best out of 5 optimization runs with 5000 objective evaluations was selected for each configuration. The configurations were ordered using the hypervolume indicator I H [31]. For NSGA-II and SPEA2, the configuration Ps/Mr/Cr = 50/0.05/0.8 was the configuration with the highest I H value. Its performance was also stable across a selection of tree models, and was therefore used in the final experiments. For MOEA/D-EAM, a smaller Ps = 30 and Mr = 0.03 in combination with Cr = 1 proved to be optimal. The MOEA/D-EAM uses additional neighborhood size Ns parameter, which was set to the recommended value Ps/10. A wider neighborhood setting Ns = Ps/5 was also experimented with, but performed worse than the default.   Tuning of methods' meta-parameters was performed with a grid search using the pools Ps ∈ {20, 30, 50, 100}, Mr ∈ {0.01, 0.03, 0.05, 1.0}, and Cr ∈ {0.7, 0.8, 0.9, 1.0}. The best out of 5 optimization runs with 5000 objective evaluations was selected for each configuration. The configurations were ordered using the hypervolume indicator I H [31]. For NSGA-II and SPEA2, the configuration Ps/Mr/Cr = 50/0.05/0.8 was the configuration with the highest I H value. Its performance was also stable across a selection of tree models, and was therefore used in the final experiments. For MOEA/D-EAM, a smaller Ps = 30 and Mr = 0.03 in combination with Cr = 1 proved to be optimal. The MOEA/D-EAM uses additional neighborhood size Ns parameter, which was set to the recommended value Ps/10. A wider neighborhood setting Ns = Ps/5 was also experimented with, but performed worse than the default.
In order to analyze the effects of individual meta-parameters in more detail, we used the ratio of non-dominated individuals (RNI) [32], which is computed between solution sets of the selected reference configuration and the modified configurations. The best solutions, in terms of RNI, were selected out of 5 optimization runs and used for configuration comparison. The effects of individual meta-parameters are reported in Table 2. Table 2. Analysis of meta-parameter effects on quality of non-dominated solution sets, measured by the RNI metric [32]. For each optimization method, the parameters Ps, Mr and Cr were individually modified. The RNI for a configuration with specific value of the observed parameter is computed with respect to other configurations using different values of the same parameter. It can be observed that, for NSGA-II and MOEA/D-EAM, the selection of Cr is less sensitive than for SPEA2. For NSGA-II, larger population sizes can be used than for the other two methods, while SPEA2 performs comparatively for population sizes 30 and 50. Because the canonical NSGA-II does not use an external archive, the initial diversity with a small population can be exhausted before a large enough pool of feasible non-dominated solutions can be constructed. MOEA/D-EAM prefers smaller mutation rates than NSGA-II and SPEA2 because the changes are propagated faster across the population via the neighborhood mechanism.

Method
The probability mass function for different mutation types was taken from the study by Strnad et al. [13], where P = P move , P add , P remove = 0.3, 0.35, 0.35 performed well. Using equal probabilities for extending and shortening the solution vectors also encourages searching with mean pruning size. An empirical lower bound value 0.9 was used for the constraint objective f 4 , i.e., no more than 10% of the tree's biomass should be removed by pruning, which is a pragmatical limit. In order to reduce the search space and prevent negligible cuts, the same values for heuristic constraints A max = 4 and m min = 10 as in [13] were used. The number of cuts in a pruning solution was also constrained empirically to the range [5,25].
The number of objective evaluations in the final experiments was set to M = 10, 000. The NSGA-II, SPEA2, and MOEA/D-EAM optimization methods were executed q = 11 times for each tree model. Within each method, the obtained q non-dominated sets were compared in order to determine the number of overall non-dominated solutions in them. The fronts were then sorted by the decreasing number of non-dominated solutions. The best, median, and worst run of both optimization methods were finally used in the analysis of results and method comparison. Figure 7 shows 3D approximations of Pareto fronts, obtained by the best runs of NSGA-II and SPEA2 for the tree models from Figure 6. Figure 8 shows the same comparison for MOEA/D-EAM and SPEA2. best, median, and worst run of both optimization methods were finally used in the analysis of results and method comparison. Figure 7 shows 3D approximations of Pareto fronts, obtained by the best runs of NSGA-II and SPEA2 for the tree models from Figure 6. Figure 8 shows the same comparison for MOEA/D-EAM and SPEA2. Figure 7. Non-dominated sets obtained by NSGA-II (blue) and SPEA2 (orange) for the corresponding trees from Figure 6. The plots include reference objective vectors for no pruning (red), cylindrical pruning (green), and rule-based pruning (black). Figure 8. Non-dominated sets obtained by MOEA/D-EAM (blue) and SPEA2 (orange) for the corresponding trees from Figure 6. The plots include reference objective vectors for no pruning (red), cylindrical pruning (green), and rule-based pruning (black).
In order to convey a better sense of shape and relation to reference solutions, the 2D projections of the Pareto front approximation for the tree in Figure 6a are shown in Figure  9. The situation is similar for other tree models. It can be observed that all of the reference solutions are dominated by the majority of non-dominated individuals found by MO methods, because their positions in objective space projections are within the dominated regions behind the non-dominated fronts. As expected, the non-selective pruning to cylindrical form improves the value of objectives f 2 and f 3 with respect to no pruning, while the rule-based pruning increases the value of objective f 1 . However, the results of pruning optimization demonstrate that simultaneous improvement on all objectives is best, median, and worst run of both optimization methods were finally used in the analysis of results and method comparison. Figure 7 shows 3D approximations of Pareto fronts, obtained by the best runs of NSGA-II and SPEA2 for the tree models from Figure 6. Figure 8 shows the same comparison for MOEA/D-EAM and SPEA2. Figure 7. Non-dominated sets obtained by NSGA-II (blue) and SPEA2 (orange) for the corresponding trees from Figure 6. The plots include reference objective vectors for no pruning (red), cylindrical pruning (green), and rule-based pruning (black). Figure 8. Non-dominated sets obtained by MOEA/D-EAM (blue) and SPEA2 (orange) for the corresponding trees from Figure 6. The plots include reference objective vectors for no pruning (red), cylindrical pruning (green), and rule-based pruning (black).
In order to convey a better sense of shape and relation to reference solutions, the 2D projections of the Pareto front approximation for the tree in Figure 6a are shown in Figure  9. The situation is similar for other tree models. It can be observed that all of the reference solutions are dominated by the majority of non-dominated individuals found by MO methods, because their positions in objective space projections are within the dominated regions behind the non-dominated fronts. As expected, the non-selective pruning to cylindrical form improves the value of objectives f 2 and f 3 with respect to no pruning, while the rule-based pruning increases the value of objective f 1 . However, the results of pruning optimization demonstrate that simultaneous improvement on all objectives is In order to convey a better sense of shape and relation to reference solutions, the 2D projections of the Pareto front approximation for the tree in Figure 6a are shown in Figure 9. The situation is similar for other tree models. It can be observed that all of the reference solutions are dominated by the majority of non-dominated individuals found by MO methods, because their positions in objective space projections are within the dominated regions behind the non-dominated fronts. As expected, the non-selective pruning to cylindrical form improves the value of objectives f 2 and f 3 with respect to no pruning, while the rule-based pruning increases the value of objective f 1 . However, the results of pruning optimization demonstrate that simultaneous improvement on all objectives is possible. In relative terms, the hardest improvements to make are with respect to objective f 1 . This can be explained by analyzing the behavior of that objective, which computes the average light exposure of flower buds. By pruning away branches, the number of buds that contribute to the objective value is reduced, but at the same time the light conditions for the remaining ones are improved and the denominator in Equation (4) is lowered. This gain vs. loss ratio of light intake can be highly discontinuous even for small changes in the solution vector, which results in multiple local optima that make optimization difficult.  Figure 10 shows realizations of pruning vectors from different regions of the nondominated set on the target tree model. It is informative to compare the pruning solutions that focus on maximization of individual objectives. For example, the pruning that maximizes the light objective tends to remove branches that cast or receive a lot of shadow, whereas the balance-oriented pruning promotes one-sided removal of wood in order to restore tree equilibrium. The pruning solution focusing on shape objective f 2 advocates branch removal from multiple sides, but contains elements of the other two prunings in order to improve on objectives f 1 and f 3 as well. The next question addressed with the experiments was the evaluation of constraint objective f 4 importance for regulating the optimization. To this end, we executed the SPEA2 optimization using a relaxed lower bound 0.8 for f 4 . The comparison of non-dominated sets for the best SPEA2 run with tight and relaxed bound is shown in Figure 11a. It is evident that further numerical improvement of objective values can be achieved by relaxing the constraint. However, the solutions tend to start overfitting on one or more pruning objectives through extensive biomass removal (Figure 11b). A proper empirical selection of bounds is, thus, necessary to confine the solutions to viable regions. objectives through extensive biomass removal (Figure 11b). A proper empirical selection of bounds is, thus, necessary to confine the solutions to viable regions.
(a) (b) Figure 11. The effect of relaxing constraint f 4 on SPEA2 optimization for the tree from Figure 6a: a) The SPEA2 non-dominated sets obtained with tight (blue) and loose (orange) bound. b) An example of excessive pruning that overfits on the balance objective. Figure 11a shows that the effect of the constraint objective is in bounding the feasible region of the objective space. By relaxing the constraint, we allow a wider spread of individual objective values, which increases the possible distance of non-dominated fronts from the origin. Such indirect specification of feasible regions is in contrast to defining constraints in terms of optimization objectives themselves, which is more common in MO [33].
The final goal of experiments is comparison of method performance. Figure 7 suggests that SPEA2 outperformed NSGA-II consistently in our tests, and also performed better than MOEA/D-EAM in general. These conclusions are supported by Tables 3 and 4, where the RNI values and the I H values are reported based on the comparison of best, median and worst runs of each method. The performance comparison of NSGA-II and R1 A6 MOEA/D-EAM shows that MOEA/D-EAM is always better in terms of RNI, but in some cases NSGA-II can achieve better results with respect to the hypervolume indicator due to a better spread of solutions in the objective space. The reason for weak performance of MOEA/D-EAM in some cases is in the inherited MOEA/D replacement strategy. It was shown that the employed neighborhood replacement scheme in MOEA/D can result in poor population diversity and premature convergence [34]. While MOEA/D-EAM addresses the reproduction step of MEOA/D, it uses the same replacement strategy, so the problem can persist.   Figure 11a shows that the effect of the constraint objective is in bounding the feasible region of the objective space. By relaxing the constraint, we allow a wider spread of individual objective values, which increases the possible distance of non-dominated fronts from the origin. Such indirect specification of feasible regions is in contrast to defining constraints in terms of optimization objectives themselves, which is more common in MO [33].
The final goal of experiments is comparison of method performance. Figure 7 suggests that SPEA2 outperformed NSGA-II consistently in our tests, and also performed better than MOEA/D-EAM in general. These conclusions are supported by Tables 3 and 4, where the RNI values and the I H values are reported based on the comparison of best, median, and worst runs of each method. The performance comparison of NSGA-II and MOEA/D-EAM shows that MOEA/D-EAM is always better in terms of RNI, but in some cases, NSGA-II can achieve better results with respect to the hypervolume indicator due to a better spread of solutions in the objective space. The reason for weak performance of MOEA/D-EAM in some cases is in the inherited MOEA/D replacement strategy. It was shown that the employed neighborhood replacement scheme in MOEA/D can result in poor population diversity and premature convergence [34]. While MOEA/D-EAM addresses the reproduction step of MEOA/D, it uses the same replacement strategy, so the problem can persist.
The common feature of SPEA2 and MOEA/D-EAM, which gives them important advantage over NSGA-II, is their use of an external archive. This finding is consistent with recent studies, in which the NSGA-II equipped with external archive was used to improve optimization results of the canonical version [35,36]. In both SPEA2 and MOEA/D-EAM, the archive is used to construct the mating pool for recombination. In the case of SPEA2, the archive is also used to update the current set of non-dominated solutions, while in MOEA/D-EAM the non-dominated set is the external archive itself. Separating the two in SPEA2 may be one reason for the slight performance difference between SPEA2 and MOEA/D-EAM. Another one is that similarity-based selection in MOEA/D-EAM results in a more locally directed search, which requires more time to escape out of large basins of similar pruning solutions. It can be concluded that the exploration mechanism of SPEA2 is better suited to the discrete optimization task at hand, and allows building better Pareto front approximations within the constrained number of fitness evaluations. Table 3. The RNI values for non-dominated sets produced by best, median, and worst runs of each optimization method. The metric is computed across the combined results of all three methods for a given run.       The advantage of NSGA-II over the other two methods is faster execution, but the differences are small because the cost of objective evaluations dominates the run-time. The average optimization run-times are reported in Table 5. The run-time variances between optimization executions of the same method on a single tree model are negligible. The slight differences between methods are in line with the expected time complexities of the methods, which are discussed in more detail in Section 4. Comparison of Tables 1 and 5 reveals that the run-time grows linearly with tree model complexity, which is expected because the objective functions in Equations (4)- (7) are linear in the number of internodes.

Discussion
The experiments demonstrate that diverse sets of non-dominated solutions can be obtained by using heterogeneous pruning objectives, which is valuable for both educational and analytical purposes. However, algorithmic pruning optimization has wider potential applicability to complement rule-based solutions in the developing field of automated pruning [37]. The progress of scanning technology and computer vision allows producing increasingly faithful digital reconstructions of real trees [38][39][40]. A wider availability of relatively low-cost devices equipped with LiDAR scanners should enable further advance of the field, as already demonstrated by some recent research [11,41]. In such scenarios, the multi-objective assessment of pruning effects may need to rely on estimation of tree parameters that are not fully observable, or cannot be captured to a sufficient level of accuracy.
The main limitation of the proposed methodology is the assumption of exact knowledge of tree properties, which are provided by the simulation environment. In order to account for the possibility of missing information, fuzzy objective measures should be incorporated into the optimization process. This would also improve the alignment of the methodology with the human perception of achieved pruning goals.
Another practical drawback of current implementation is its run-time, which needs to be improved for the use of methodology in interactive sessions. The differences between the optimization methods in this respect are small, because the run-time is governed by objective evaluations. The definitions of objectives in Section 2.3 indicate that their evaluation should increase approximately linearly with problem size, i.e., the number of tree internodes. The objective f 1 in Equation (4) sums over the flower buds, whose number is linear in the number of internodes. The most time consuming part of the objective f 2 computation is the generation of a convex hull. This can be done in O(nh), where n is the number of branch tips and h is the number of points in the hull. While n increases approximately linearly with the tree size, h was determined to be within a small constant factor for all trees. Objectives f 3 and f 4 are directly related to the number of internodes, so their computation is also linear. The linear relationship of run-time and problem size has been experimentally confirmed, as presented in Table 5. The small run-time differences between methods arise from their average time complexities, which are O(N 2 ) for NSGA-II, O(N 2 log(N)) for SPEA2 and O(N 2 + N M) for MOEA/D-EAM, where N is the population size and M is the archive size. Further optimization and parallelization of objective value computation is one of priorities for future framework development.
The presented use of additional objective constraints to prevent overfitting of pruning on actual target objectives is a generally applicable concept. It can be used in other MO problems where the bounds are hard to specify on optimization objectives directly. In the proposed approach to pruning optimization, the objective constraints complement heuristic decision constraints in the search space. However, the violations of the latter can be easily detected and corrected within the optimization loop. In this way the waste of computational resources on evaluation of infeasible solutions can be avoided. Violations of objective constraints, on the other hand, can only be established after their evaluation, and it is usually not clear how to resolve them. A strategy for handling the infeasible solutions in the population is therefore required. In the currently implemented strategy, all solutions that violate the constraint objective become equally dominated due to the assignment of low score. However, a finer distinction between nearly feasible solutions and truly bad ones with adaptively increasing constraint violation penalty could improve the exploration efficiency of the method.

Conclusions
In the paper, we introduced the use of heterogeneous objectives within an MO framework for virtual tree pruning. We also proposed the use of constraint objectives, which are not the optimization targets, but are employed as heuristic constraints in the objective space of pruning solutions. Finally, a performance analysis of three popular MO methods was performed on the presented discrete optimization problem.
The experiments have shown that, by using additional objectives to place heuristic constraints on search in the fitness space, the overfitting of solutions to individual objectives can be restricted. It was demonstrated that the use of an external archive allows SPEA2 and MOEA/D-EAM to achieve better performance than NSGA-II in the studied problem, while the slight advantage of NSGA-II is the shorter run-time.