Figures 9 and 10 show the Pareto-optimal fronts found in each one of the cases (Case 1 and Case 2) using the three sets of objectives: (i) all objectives; (ii) objectives Q, Power, WATS, and T; (iii) objectives Q, WATS, and T. The results are, apparently, very similar when comparing the cases with six and four objectives. In the other situation, with three objectives, the multi-objective optimization algorithm is clearly lost, since the final solution found alternates in the different runs between one type of screw and the other (i.e., between the CS and the MBS). The results for Cases 2 and 3 are very similar to those presented here and, thus, no specific discussion is made here.

By using the 11 runs performed for each case studied, the Hypervolume (HV) and the Inverted Generational Distance (IGD) were applied and the results are presented in Table 7, where it is possible to see the average and the percentage of losses when the number of objectives is reduced [42–44]. To calculate IGD, all Pareto-optimal solutions found in each run were put together in a pool and the non-dominated solutions of this pool were used for comparison.

As shown in Table 7, it is possible to conclude that the use of four objectives (Q, Power, WATS, and T) does not significantly deteriorate the final solutions found, the maximum difference found is 11.6%, which, for a process like the extrusion process, and taking into account a final population of 100 solutions, is not expressive. Additionally, the differences in the IGD value are too small, indicating that at least the solutions found for the case of four objectives are near the best solution found in the 11 runs. The results found for the situation with three objectives corroborate the results shown in Figures 9 and 10.

**Figure 9.** Pareto-optimal fronts after 100 generations for the pair of objectives identified in Figure 8 for Case 1.

**Figure 10.** Pareto-optimal fronts after 100 generations for the pair of objectives identified in Figure 8 for Case 7.

**Table 7.** Performance comparison using Hypervolume and IGD for the total number of objectives and the automatic reduction to four and three objectives (between brackets the standard deviation, and loss percentage relative to six objectives) for the four cases studied.

Finally, it is important to point out that during this process, the DM(s) play an important role in the procedure. Indeed, they have some intervention when selecting the objectives. For example, it is necessary to opt for Q or L, two objectives from the same cluster having apparently the same importance in the process. In this case, an informed DM will make the option for Q because this objective is the output of the machine and is directly linked with the economic issue of the problem, while L is the length for melting that is related to the quality of the product obtained; however, this quality is also quantified by WATS, which was already selected by the algorithm. This example shows the importance of the DM(s) that simultaneously help the optimization process and are informed about the process of obtaining the results.

### **5. Conclusions**

A methodology for reducing the number of objectives for many-objective optimization problems using population-based algorithms is proposed. This approach, based on machine learning, is an improvement over similar state-of-the-art methodologies; namely, it allows analysis of the relations variable–variable and variable–objective relations (and not only objective–objective), does not need kernel function choice and parameters optimization, allows for obtaining explainable solutions to assist the decision maker with interpreting the results, its time complexity is also low, and it supports theoretical and empirical sample sizes.

The approach showed its potential to reduce the number of objectives by capturing the complex relations between the different objectives with an additional possibility, which is to capture the objective-variable relations. This is done by applying the methodology to a set of benchmark and real-world problems. The comparison of the Pareto-optimal fronts obtained with another machine learning approach in the literature allows for the conclusion that its performance is very competitive, but with the great advantage of being much easier to use. Additionally, there is the possibility of strong interaction with the DM(s).

The application of the proposed approach to a difficult real-world problem has proven that it is automatically possible to reduce the number of objectives by losing only around ten percent of the Pareto-optimal frontier obtained, for the case of 100 individuals in the population. The use of a second possibility, which is to require the intervention of the decision maker during the process, e.g., when selecting the objectives to be considered in the optimization, can be very useful because the person interested can see how the process works and interpret the results obtained. Finally, an important characteristic of the method proposed is the capacity to explain the solutions found.

**Author Contributions:** Conceptualization, A.G.-C., F.M. and A.D.; methodology, A.G.-C., F.M. and A.D.; software, P.C.; investigation, A.G.-C. and A.D.; resources, A.G.-C.; data curation, A.G.-C.; writing—original draft preparation, A.G.-C. and A.D.; writing—review and editing, A.G.-C.; visualization; supervision, A.G.-C. and A.D.; project administration, A.G.-C. and A.D.; funding acquisition, A.G.-C., F.M. and A.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by POR Norte under the PhD Grant PRT/BD/152192/2021. The authors also acknowledge the funding by FEDER funds through the COMPETE 2020 Programme and National Funds through FCT (Portuguese Foundation for Science and Technology) under the projects UID-B/05256/2020, and UID-P/05256/2020, the Center for Mathematical Sciences Applied to Industry (CeMEAI) and the support from the São Paulo Research Foundation (FAPESP grant No 2013/07375-0, the Center for Artificial Intelligence (C4AI-USP), the support from the São Paulo Research Foundation (FAPESP grant No 2019/07665-4) and the IBM Corporation.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Single-Loop Multi-Objective Reliability-Based Design Optimization Using Chaos Control Theory and Shifting Vector with Differential Evolution**

**Raktim Biswas † and Deepak Sharma \*,†**

Department of Mechanical Engineering, Indian Institute of Technology Guwahati, Guwahati 781039, Assam, India

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +91-361-2582661

† These authors contributed equally to this work.

**Abstract:** Multi-objective reliability-based design optimization (MORBDO) is an efficient tool for generating reliable Pareto-optimal (PO) solutions. However, generating such PO solutions requires many function evaluations for reliability analysis, thereby increasing the computational cost. In this paper, a single-loop multi-objective reliability-based design optimization formulation is proposed that approximates reliability analysis using Karush-Kuhn Tucker (KKT) optimality conditions. Further, chaos control theory is used for updating the point that is estimated through KKT conditions for avoiding any convergence issues. In order to generate the reliable point in the feasible region, the proposed formulation also incorporates the shifting vector approach. The proposed MORBDO formulation is solved using differential evolution (DE) that uses a heuristic convergence parameter based on hypervolume indicator for performing different mutation operators. DE incorporating the proposed formulation is tested on two mathematical and one engineering examples. The results demonstrate the generation of a better set of reliable PO solutions using the proposed method over the double-loop variant of multi-objective DE. Moreover, the proposed method requires 6×–377× less functional evaluations than the double-loop-based DE.

**Keywords:** multi-objective reliability-based design optimization; shifting vector approach; reliability analysis; chaos control theory; differential evolution

### **1. Introduction**

The design optimization mostly keeps design variables and parameters deterministic. It ignores the fact that uncertainties can arise owing to manufacturing variations, dimensional inaccuracy, boundary conditions, material properties, and improper loading conditions, which can lead to the infeasibility of the solution obtained through deterministic optimization. Therefore, it is necessary to consider these uncertainties in designing the process to maintain safety and the quality of the solution. Reliability-based design optimization (RBDO) [1,2] is a mathematical tool that is used for obtaining such reliable optimal solutions for problems involving uncertainties. It also enables engineers to identify solutions effectively for complex applications in the fields of the automotive, civil, mechanical, and aerospace industries [3,4]. In RBDO, the uncertainties are manifested by converting the deterministic constraints to probabilistic constraints. This is accomplished by applying a probability operator to performance functions or to limit-state functions in the literature. A generalized single-objective RBDO formulation is given in Equation (1).

> Minimize *f*(*μ***X**), subject to *P*[*Gi*(**X**) ≥ 0] ≤ *PT fi* = Φ(−*βT i* ), *i* = 1, . . . , *I*, *μ*(L) **X** ≤ *μ***X** ≤ *μ*(U) **X** , (1)

**Citation:** Biswas, R.; Sharma, D. Single-Loop Multi-Objective Reliability-Based Design Optimization Using Chaos Control Theory and Shifting Vector with Differential Evolution. *Math. Comput. Appl.* **2023**, *28*, 26. https://doi.org/ 10.3390/mca28010026

Academic Editors: Carlos Coello, Erik Goodman, Kaisa Miettinen, Dhish Saxena, Oliver Schütze and Lothar Thiele

Received: 3 July 2022 Revised: 30 January 2023 Accepted: 15 February 2023 Published: 17 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *f*(*μ***X**) is the objective function, *Gi*(**X**) is the *i*-th performance/constraint function, and *μ***X** is the mean value vector of random variable vector **X** ∈ R*n*, where *n* is the number of random design variables. *L* and *U* in the superscript of *μ***X** represent the lower and upper limits of the vector. Φ(·) represents the standard normal cumulative distribution function, *βT i* is the target reliability index of the *i*-th performance function, and *P*[·] is the probability operator that represents the failure probability of performance function (*Gi*(**X**) ≥ 0) that should be less than the target failure probability (*PT fi* ).

Equation (1) demonstrates that solving a single-objective RBDO requires a nested-loop procedure [2], where the outer optimization loop involves the inner-loop for reliability analysis. The reliability analysis can be performed using simulation-based methods [5] and analytical methods [6] on probabilistic performance function to obtain its failure probability. The simulation-based methods show better accuracy with an expense of computational cost [7], such as Monte Carlo simulation (MCS) [5], subset simulation [8], importance sampling [9], and Latin-hypercube sampling [9]. On the other hand, analytical methods are known for their computational efficiency, such as most-probable point (MPP)-based methods, in which the sub-optimization problem is solved for each performance function to obtain their respective MPP. The MPP-based methods can be broadly divided into the performance measurement approach (PMA) [10] and the reliability index approach (RIA) [6]. The optimum solution obtained using PMA and RIA is known as the most probable target point (MPTP) and the most probable failure point (MPFP), respectively. Many advanced methods have been developed to estimate the MPTP and MPFP of performance functions, and they are categorized as double-loop methods, decoupled-loop methods, and singleloop methods.

The classical double-loop methods [11,12] involve a nested optimization loop, where the inner-loop performs reliability analysis and the outer-loop is used for obtaining design solutions. All the random variables are transformed to standard normal variables [13] for performing reliability analysis. Since the nested optimization loop is computationally expensive, the reliability analysis loop (inner-loop) is decoupled and performed separately in decoupled-loop methods [14–17]. Some advanced and efficient reliability-based frameworks were also proposed based on isogeometric analysis [18,19]. The reliability analysis itself is considered as an computationally expensive procedure. Therefore, single-loop methods [20] have been proposed, in which approximate reliability analysis is performed. Different concepts such as Karush-Kuhn Tucker (KKT) conditions and quantile approximation are used to approximate MPTP that can eliminate the reliability analysis loop. The adaptive conjugate single-loop approach (AC-SLA) [21], the enhanced single-loop method (ESM) [22], the chaotic single-loop approach (CSLA) [23], the single-loop shifting vector method (SLShV-CG) [24], the sequential single-loop reliability optimization and confidence analysis method (SROCA) [25], and the approximate single-loop chaos control method (ASLCC) [26] are a few recently developed single-loop methods. Recently, some efficient evolutionary RBDO methods are also proposed to obtain the global reliable solution [27,28].

It has been found that many real-world engineering problems consist of more than one objective, which are conflicting in nature [29], and can also have uncertainties. Evolutionary algorithms are found to be promising for solving deterministic multi-objective optimization problems (MOOPs) because they can generate Pareto-optimal (PO) solutions in one run. However, these evolutionary algorithms need to be modified for generating reliable PO solutions for multi-implemented as a design optimization algorithm, and inverse reliability was performed. objective reliability-based design optimization (MORBDO) problems. To address uncertainty in MORBDO, Deb et al. [3] used a non-dominated sorting genetic algorithm (NSGA-II) [30] for design optimization, and Fast RIA for reliability analysis. A multi-objective differential evolution (MODE) [31] was also Simulation-based techniques are also used for reliability analysis and are coupled with double-loop methods. For example, a radial basis function was used for approximating the responses of the

performance function and was coupled with MCS to implement reliability analysis. NSGA-II was used to obtain PO solutions for solving the multi-objective and multi-case [32] RBDO problem. In another study, MCS and NSGA-II were coupled with entropy weighted grey relational analysis for design optimization [33] to solve the control arm problem. The multi-objective optimization design of the control arm was carried out using the Kriging surrogate model. Sun et al. [34] proposed a radial basis function-based surrogate modeling that was implemented with Latin-hypercube sampling for sensitivity analysis. MCS and multi-objective particle swarm optimization (PSO) were coupled for obtaining the reliable PO solutions. In another study, a multiple response surface method-based artificial neural network was implemented for reliability analysis [35], and a dynamic multi-objective particle swarm optimization algorithm was proposed for obtaining PO solutions. A worst-case scenario was used with fuzzy sets for reliability analysis, and a realcoded population-based incremental learning [36] was implemented with DE for obtaining the PO solutions. A multi-objective robust optimization [37] was proposed, in which the design problems consisted of parametric uncertainties involving both random and interval variables. NSGA-II was implemented to generate robust PO solutions, and MCS was performed to evaluate the impact responses of the mixed uncertainties. Constrained NSGA-II was also implemented to solve the MORBDO problem [38]. It was coupled with the hybrid method using the Kriging surrogate metamodel for reliability analysis.

A time-dependent reliability-based robust design optimization (TRBRDO) problem [39] was solved using NSGA-III [40] and the dimension reduction method. It was developed by constructing an extreme value model using the sparse grid-based stochastic collocation method for time-dependent reliability analysis. A Bayesian multi-objective RBDO [41] was proposed to solve problems involving aleatory and epistemic uncertainties. Multi-objective PSO was implemented for obtaining PO solutions, and Bayesian interference was used for reliability analysis. Another method using nested loop was proposed to solve RBDO problems [42], in which the outer-loop was performed using multi-objective PSO, and the inner-loop was solved using surrogate modeling with MCS sampling. A two-layer nested optimization problem was proposed based on a decoupling strategy. The inter-generation projection genetic algorithm was employed in the inner-loop, and the multi-objective genetic algorithm [43] was implemented at the outer-loop for solving the MORBDO problem. Another multi-objective RBDO [44] was solved by converting it into a single-objective RBDO problem. This was achieved by assigning weights to the objectives based on quantitative analysis and evidence theory. The reliability analysis was estimated using the PMA method.

From the literature, it can be seen that most of the MORBDO methods focus on PMA, RIA, MCS, or surrogate modeling for reliability analysis, and they are based on double-loop or decoupled-loop methods, which make them computationally expensive. Since evolutionary algorithms are population-based methods and require many functional evaluations, a single-loop method for solving MORBDO can improve the computational efficiency. Moreover, single-loop methods that are solved using steepest descent search to estimate MPTP are often stuck with periodic oscillation [26,45] for highly nonlinear functions. This leads to the motivation of this paper, in which a new MORBDO formulation is proposed, based on adaptive multi-objective DE. An adaptive mutation scheme is used for selecting different variants of mutations for exploration in the search space. Both trial and target vectors take part in the MORBDO formulation to estimate the reliable PO solutions. The following are the contributions of the paper.

• The formulation is further developed by incorporating target and trial vectors of differential evolution for better exploration of the search space.

The proposed method is tested on three benchmark examples from the literature. The results are compared with a double-loop variant of multi-objective differential evolution using PMA for reliability analysis.

The organization of the paper is as follows. In Section 2, a brief discussion on multiobjective RBDO, PMA, chaos control method, single-loop method, and shifting vector approach are presented. The proposed single-loop multi-objective reliability-based design optimization method is discussed in Section 3, along with its implementation. The adaptive mutation scheme and the detailed steps of multi-objective differential evolution are also discussed in this section. Numerical examples are solved and discussed in Section 4. Finally, the paper is concluded in Section 5 with a note on future work.

### **2. Preliminaries**

*2.1. Multi-Objective Reliability-Based Design Optimization*

A generalized MORBDO formulation can be written as

$$\begin{array}{ll}\text{Minimize} & f\_{\mathfrak{m}}(\boldsymbol{\mu}\_{\mathbf{X}}), & m = 1, \ldots, M, \\ \text{subject to} & P[G\_i(\mathbf{X}) \ge 0] \le P\_{f\_i}^T = \Phi(-\boldsymbol{\beta}\_i^T), & i = 1, \ldots, I, \\ & \boldsymbol{\mu}\_{\mathbf{X}}^{(\mathcal{L})} \le \boldsymbol{\mu}\_{\mathbf{X}} \le \boldsymbol{\mu}\_{\mathbf{X}}^{(\mathcal{U})}, \mathbf{X}^{(\mathcal{L})} \le \mathbf{X} \le \mathbf{X}^{(\mathcal{U})}, \end{array} \tag{2}$$

where *fm*(·) is the *m*-th conflicting objective function that is written using the mean value (*μ***X**) of the random variable (**X**). **X**(L) and **X**(U) are the upper and lower limits on **X**. Solving Equation (2) generates a set of reliable PO solutions in the design space. The reliability analysis is performed on the probabilistic performance function to estimate the failure probability by solving a multidimensional integral, as given in Equation (3).

$$P\_{f\_i} = P[\mathbb{G}\_i(\mathbf{X}) \ge 0] = \int \cdots \int\_{G\_i(\mathbf{X}) \ge 0} f\_{\mathbf{X}}(\mathbf{X}) d\mathbf{X} \tag{3}$$

where *f***X**(**X**) is the joint probability density function of **X**. Solving this multidimensional integral is difficult, and therefore, it is approximated with reliability analysis [7]. The firstorder reliability method (FORM) [6] and second-order reliability method (SORM) [46] are analytical methods for reliability analysis. Both FORM and SORM estimate the reliability index *β* that represents the minimum distance from the origin to the performance function in the standard normal space. The reliability index *β* can be obtained by solving a suboptimization problem, and the reliability (R) can be estimated using Φ(*β*) (R = 1 − *Pf* = 1 − Φ(−*β*) = Φ(*β*)). Due to its computational efficiency and stability in generating a reliable solution, PMA is widely used to solve the sub-optimization problem [47].

### *2.2. Performance Measure Approach (PMA)*

PMA estimates the failure probability of performance function *G*(**X**) by finding MPTP in the standard normal space (*U*-space). After transforming *G*(**X**) to the *U*-space using the Rosenblatt transformation [13], the MPTP can be estimated using the steepest descent direction. When all the random variables are independent, the joint cumulative distribution function (CDF) is calculated via the product of the marginal CDFs. The Rosenblatt transformation is given as

$$\Phi(\boldsymbol{u}\_i) = \mathcal{F}\_{\boldsymbol{X}\_i}(\mathbf{x}\_i) \implies \boldsymbol{u}\_i = \Phi^{-1}(\mathcal{F}\_{\boldsymbol{X}\_i}(\mathbf{x}\_i)),\tag{4}$$

where *FXi* (*xi*) is the marginal CDF of *Xi* and Φ(·) is the CDF of the standard normal random variable. After transforming variables to the standard normal space by using Equation (4), MPTP is calculated by performing the following sub-optimization problem.

$$\begin{array}{ll}\text{Minimize} & G(\mathbf{U}),\\\text{subject to} & ||\mathbf{U}|| = \boldsymbol{\beta}^T.\end{array} \tag{5}$$

where **U** is the random variable in the standard normal space, and *βT* is the target reliability index for the performance function *G*(**U**). To efficiently obtain the optimum solution of Equation (5), the advanced mean value algorithm is used and the expression is presented in Equation (6).

$$\mathbf{U}^{(k+1)} = \boldsymbol{\beta}^T \frac{\nabla G(\mathbf{U})}{||\nabla G(\mathbf{U})||}. \tag{6}$$

If the performance function value at MPTP is less than or equal to zero, it is satisfied for the given target reliability, as presented in Equation (2).

### *2.3. The Chaos Control Method*

It has been observed that PMA performs well for simple nonlinear performance functions, but it fails to converge for highly nonlinear performance functions. To overcome this issue, chaos control theory [45] was proposed based on a stability transformation method [48]. The modification is achieved while updating the iterative point **U**(*k*+1) of Equation (5). The formulation for estimating the iterative point via the chaos control (CC) method is as follows.

$$\begin{aligned} \mathbf{U}\_{\mathbf{CC}}^{(k+1)} &= \mathbf{U}\_{\mathbf{CC}}^{(k)} + \lambda \mathbf{C} [\mathbf{F}(\mathbf{u}^{(k)}) - \mathbf{U}\_{\mathbf{CC}}^{(k)}], \\ \mathbf{F}(\mathbf{u}^{(k)}) &= \mathbf{U}^{(k+1)} = \boldsymbol{\beta}^T \frac{\nabla G(\mathbf{U})}{||\nabla G(\mathbf{U})||}, \end{aligned} \tag{7}$$

where U(*k*) *CC* is the MPTP calculated using CC method in the *k*-th iteration; **C** is the involutory matrix with only one element in each row and is assumed as identity matrix **I** for simplicity. The matrix **C** is usually selected to stabilize the unstable fixed point of the chaotic dynamical system in Equation (7). The chaos control factor *λ* is determined according to the eigenvalues of the original system's Jacobian matrix, and the value is considered within interval [0, 1]. When *λ* is considered as one, the formulation of the CC method is similar to Equation (5) and can have the same issue as discussed earlier. Therefore, a small value of *λ* is considered for stable convergence. **F** is the vector of the response function that is estimated via nonlinear mapping with respect to the iterative values of **U**(*k*+1) , as shown in Equation (7). Although the CC method eliminates the issue of oscillation in the convergence of MPTP, it is considered to be an inefficient process. Therefore, a modified chaos control (MCC) [12] was proposed. The modification is achieved by extending the iterative search to the *β*-hypersphere that is at the constraint boundary in the standard normal space. Thus, MPTP is located on the constraint boundary, and convergence is improved by controlling the tangential step size instead of the radial step size, which was the case for the CC method. The formulation of MCC is given as

$$\begin{aligned} \boldsymbol{\tilde{\mathfrak{n}}}^{(k+1)} &= \mathbf{U}\_{\text{CC}}^{(k)} + \lambda \mathbf{C} [\mathbf{F}(\mathbf{u}^{(k)}) - \mathbf{U}\_{\text{CC}}^{(k)}], \\ \mathbf{U}\_{\text{MCC}}^{(k+1)} &= \boldsymbol{\beta}^{T} \frac{\boldsymbol{\tilde{\mathfrak{n}}}^{(k+1)}}{||\boldsymbol{\tilde{\mathfrak{n}}}^{(k+1)}||}. \end{aligned} \tag{8}$$

where **n**˜ *k* is the modified search direction updated using **U**(*k*+1) *CC* of Equation (7). **U**(*k*+1) *MCC* is the MPTP evaluated using the MCC method.

### *2.4. Single-Loop Method*

The single-loop method (SLM) [20] has been proposed to approximate the reliability analysis of the double-loop method, and establish an equivalent deterministic performance function that is computationally efficient. The approximate MPTP is estimated by using the KKT optimality conditions of Equation (5), and is given in Equation (9).

$$
\nabla G(\mathbf{U}) - \hat{\lambda}\nabla H(\mathbf{U}) = 0,\tag{9}
$$

where *λ*ˆ is the Lagrange multiplier, and *H*(**U**) = **U** 2 − *βT i* 2 after squaring both sides of the equality constraint of Equation (5). Using Equation (9) and ∇*H*(**U**) = 2**U** yields ∇*G*(**U**) − 2**U***λ*ˆ = 0. After simplification, **U** can be written as ∇*G*(**U**) 2*λ*ˆ , and multiplying it with ∇*G*(**U**) in the numerator and denominator, and further simplifying, we obtain

$$\mathbf{U} = \frac{||\nabla G(\mathbf{U})||}{2\lambda} \frac{\nabla G(\mathbf{U})}{||\nabla G(\mathbf{U})||} = \boldsymbol{\beta}^T \mathbf{a}\_\prime \tag{10}$$

where *α* = ∇*G*(**U**) ∇*G*(**U**) is the unit gradient direction, and *βT* = ∇*G*(**U**) 2*λ*ˆ is a constant at the optimal solution **U**∗. The gradient is calculated in *U*-space and the random design variables lie in the *X*-space. Therefore, the transformation from *X*-space to *U*-space is used for the evaluation of approximate MPTP, using the following relationship.

$$
\mathbf{X} = \mu\_{\mathbf{X}} + \sigma\_{\mathbf{X}} \mathbf{U},
\tag{11}
$$

where *σ***X** is the standard deviation of **X**. Substituting **U** from Equation (10) in Equation (11) and using the chain rule, we obtain MPTP in the *X*-space as

$$\mathbf{X}\_{\rm MPTP} = \mu\_{\mathbf{X}} + \sigma\_{\mathbf{X}}\beta\mathbf{a} = \mu\_{\mathbf{X}} + \sigma\_{\mathbf{X}}\beta^T \frac{\sigma\_{\mathbf{X}}\nabla\chi\mathbf{G}(\mathbf{X})}{\|\sigma\_{\mathbf{X}}\nabla\chi\mathbf{G}(\mathbf{X})\|} \tag{12}$$

where **X***MPTP* is the MPTP of the performance function *G*(**X**).

### *2.5. Shifting Vector Approach*

The concept of the shifting vector (**S**(*k*) *i* ) has been proposed [14] to decouple the doubleloop structure of the RBDO problem. It separates the optimization and reliability analysis loop and performs it sequentially in the sequential optimization and reliability assessment (SORA) [14] method. Using this process, the computational efficiency of SORA has been improved as compared to the double-loop method. The concept of the shifting vector is used to shift the violated performance function towards the feasible direction. It is given as

$$\mathbf{S}\_{i}^{(k)} = \mu\_{\mathbf{X}}^{(k-1)} - \mathbf{X}\_{i, \text{MPTP}}^{(k-1)} \tag{13}$$

where (**S**(*k*) *i* ) is the shifting vector at the *k*-th iteration, **X**(*k*−1) *i*,*MPTP* is the MPTP for the *i*-th constraint, and *μ*(*k*−1) **X** is the mean of the random variable **X** in the (*k* − 1)-th iteration. Figure 1 shows the schematic diagram of the shifted constraint based on the MPTP. It can be seen that (**S**(1) *i* ) is estimated based on **X**(1) *i*,*MPTP* and *μ*(1) **X** , and the shifted constraint is evaluated at *μ*(1) **X** − (**S**(1) *i* ) until the reliability of the constraint is achieved. Here, the shifting vector (**S**(*k*) *i* ) is generated via an iterative process that helps to estimate the feasibility of the performance function until its reliability is satisfied.

**Figure 1.** Shifting vector approach.

### **3. The Proposed Method and Its Implementation**

*3.1. Single-Loop MORBDO Formulation Using Chaos Control and the Shifting Vector Approach*

The single-loop MORBDO formulation can be written using the approximate MPTP given in Equation (12) as

$$\begin{aligned} \text{Min.} \, f\_{\mathbf{m}}(\mu\_{\mathbf{X}}), \qquad & m = 1, \dots, M, \\ \text{s.t.:} \, \mathbb{G}\_{i}(\mathbf{X}\_{i, \text{MPTP}}^{(k)}) &\le 0, \qquad & i = 1, \dots, I, \\ \text{where } \mathbb{X}\_{i, \text{MPTP}}^{(k)} &= \mu\_{\mathbf{X}}^{(k)} + \boldsymbol{\beta}\_{i}^{T} \sigma\_{\mathbf{X}} \mathbf{a}\_{i, \mathbf{X}'}^{(k)} \\ \text{a}\_{i, \mathbf{X}}^{(k)} &= \frac{\sigma\_{\mathbf{X}} \, \nabla G\_{i, \mathbf{X}}(\mathbf{X}\_{i, \text{MPTP}}^{(k-1)})}{||\sigma\_{\mathbf{X}} \, \nabla G\_{i, \mathbf{X}}(\mathbf{X}\_{i, \text{MPTP}}^{(k-1)})||}, \\ \mu\_{\mathbf{X}}^{(\text{L})} &\le \mu\_{\mathbf{X}} \le \mu\_{\mathbf{X}}^{(\text{U})}. \end{aligned} \tag{14}$$

where **X**(*k*) *i*,*MPTP* is the approximate MPTP of the '*i*' performance function at the *k*-th iteration, and *α*(*k*) *i*,**X** is the unit gradient vector of the performance function '*i*' with respect to random variable (**X**). In Equation (14), the probabilistic performance functions of Equation (2) are converted into deterministic performance functions, which eliminate the MPTP search of the inner-loop at every iteration. Thus, the computational efficiency can be improved significantly. It is to be noted that the steepest descent search is used to evaluate the approximate MPTP, which has a tendency to oscillate during convergence [45].

In the proposed formulation, chaos control theory replaces the steepest descent search for approximating MPTP. The concept of the shifting vector approach is incorporated to formulate a novel single-loop MORBDO formulation, as shown in Equation (15).

$$\begin{aligned} \text{Min.} \, f\_{\mathfrak{m}}(\mu\_{\mathbf{X}}), \quad & i = 1, \ldots, M, \\ \text{s.t.:} \, \mathbf{G}\_{i}(\mathbf{F}^{(k)}) \le 0, \quad & i = 1, \ldots, I, \\ \text{where } \mathbf{F}^{(k)} = \begin{cases} \mathbf{X}\_{i, MTPP}^{(k)}, & \forall \text{ target vectors,} \\ \mathbf{\mu}\_{\mathbf{U}}^{(k+1)} - \mathbf{S}\_{i}^{(k+1)}, & \forall \text{ trial vectors,} \\ \mathbf{\mu}\_{\mathbf{X}}^{(k+1)} - \mathbf{X}\_{i, MTPP}^{(k)}, & \\ \mathbf{\mathcal{X}}\_{i, MTPP}^{(k)} = \mathbf{T}^{-1}(\mathbf{U}) = \boldsymbol{\mu}\_{\mathbf{X}}^{(k)} + \sigma\_{\mathbf{X}} \mathbf{U}\_{i, SLCC}^{(k)} \\ \boldsymbol{\mu}\_{\mathbf{X}}^{(L)} \le \boldsymbol{\mu}\_{\mathbf{X}} \le \boldsymbol{\mu}\_{\mathbf{X}}^{(L)}. \end{aligned} \tag{15}$$

where **U**(*k*) *i*,*SLCC* is the approximate MPTP in the *U*-space that is estimated using the MCC method. *μ*(*k*+1) **U** is the trial vector of differential evolution in the *U*-space in the (*k* + 1)-th iteration. In the proposed formulation, the performance function *Gi*(**Ψ**(*k*) ) includes both **X**(*k*) *i*,*MPTP* and (*μ*(*k*+1) **U** − **S**(*k*+1) *i* ), which are used for evaluating the performance function for each target vector and trial vector, respectively. The vector (*μ*(*k*+1) **U** − **S**(*k*+1) *i* ) shifts the violated performance function towards a feasible direction for the population of trial vectors. **U**(*k*) *i*,*SLCC* in the standard normal space is given in Equation (16).

$$\mathbf{U}\_{i,SLCC}^{(k)} = \beta\_i^T \frac{\mathbf{U}\_i^{(k-1)} + \lambda\_i^{(k)} \mathbf{C}[\mathbf{U}\_i^{(k)} - \mathbf{U}\_i^{(k-1)}]}{||\mathbf{U}\_i^{(k-1)} + \lambda\_i^{(k)} \mathbf{C}[\mathbf{U}\_i^{(k)} - \mathbf{U}\_i^{(k-1)}]||} \tag{16}$$

where **U**(*k*) *i* and **U**(*k*−1) *i* are the MPTPs estimated for the *i*-th constraint in the *k*-th and (*k* − 1) th generations, respectively. The value of **U**(*k*) *i*,*SLCC* is calculated after the transformation, as given in Equation (17).

$$\mathbf{U}\_{i}^{(k)} = \mathbf{T}(\mathbf{X}\_{i}^{(k)}) = (\mathbf{X}\_{i}^{(k)} - \mu\_{\mathbf{X}}) / \sigma\_{\mathbf{X}}.\tag{17}$$

The proposed single-loop MORBDO formulation given in Equation (15) is developed based on a single-loop methodology that eliminates the integrated reliability analysis involved in double-loop formulation, as given in Equation (2). The approximated formulation for reliability analysis is established through KKT optimality conditions, where the search direction is calculated by using modified chaos control theory. Furthermore, the shifting vector is integrated with the single-loop MORBDO formulation that uniquely involves the target and trial vectors of differential evolution.

### *3.2. Multi-Objective Differential Evolution with Adaptive Mutation Scheme*

Differential evolution (DE) [49] is a population-based meta-heuristic algorithm that works with a set of vectors and optimizes an optimization problem by iteratively improving each vector based on an evolutionary process. It explores the design space by maintaining a population of vectors and creating new vectors by combining existing ones. It starts with a random generation of vectors, which are referred to as target vectors, *μ*(*k*) **X** (*t*), in which *t* represents the *t*-th target vector, and *k* represents the *k*-th generation counter. Since DE is used for solving the MORBDO problem, the notation for vector is kept the same as the mean value of the random variable. Each target vector (*μ*(*k*) **X** (*t*)) is transformed to the mutant vector (*μ*(*k*+1) **V** (*t*)) using the randomly chosen vectors (*μ*(*k*) **r**1 (*t*)), (*μ*(*k*) **r**2 (*t*)) and (*μ*(*k*) **r**3 (*t*)). In this paper, an adaptive mutation scheme is used, in which the mutation vector (*μ*(*k*+1) **V** (*t*)) is generated, either by using a random vector or the best vector. The scheme for generating (*μ*(*k*+1) **V** (*t*)) is given in Equation (18).

$$\mu\_{\mathbf{V}}^{(k+1)}(t) = \begin{cases} \mu\_{\mathbf{r1}}^{(k)}(t) + \hat{\mathbf{r}} \times (\mu\_{\mathbf{r2}}^{(k)}(t) - \mu\_{\mathbf{r3}}^{(k)}(t)), & \zeta > \epsilon, \\\\ \mu\_{\mathbf{best}}^{(k)}(t) + \hat{\mathbf{r}} \times (\mu\_{\mathbf{r2}}^{(k)}(t) - \mu\_{\mathbf{r3}}^{(k)}(t)), & \text{otherwise,} \end{cases} \tag{18}$$

where **r1** = **r2** = **r3** are the three randomly chosen vectors from the current population, and *F*ˆ is the scaling factor. The variant "DE/rand/bin/1" is found to be effective in exploring the search space during the initial generations because the mutant vector is generated a using random vector. When DE starts converging towards the Pareto-optimal front, the "DE/best/bin/1" variant replacing *μ*(*k*) **r**1 (*t*) to *μ*(*k*) **best**(*t*) can improve the convergence. The *μ*(*k*) **best**(*t*) vector for each target vector is found by calculating the Euclidean distance of the *t*-th target vector with respect to all non-dominated target vectors in the objective

space. The closest non-dominated target vector is selected as *μ*(*k*) **best**(*t*) for the *t*-th target vector. Since both the variants have their own merits, a heuristic convergence parameter (*ζ*) is proposed that can help DE to use either of these variants, depending on the user-defined parameter *-*. The parameter *ζ* is calculated using the hypervolume (HV) performance indicator [50] that is given as

$$\mathcal{L} = I\_H^{(k)} - I\_H^{(k-1)}\, \text{.} \tag{19}$$

where *I* (*k*) *H* and *I* (*k*−1) *H* are the hypervolume calculated with respect to the non-dominated target vectors in the (*k*) and (*k* − 1) generations. It is noted that the non-dominated target vectors in the (*k* − 1) and (*k*) generations are normalized together for estimating the hypervolume with respect to the dominated point. Thereafter, the trial vector (*μ*(*k*+1) **U** (*t*)) is created for each target vector (*μ*(*k*) **X** (*t*)), which is given as

$$\mu\_{\mathbf{U}}^{(k+1)}(t\_{\bar{j}}) = \begin{cases} \mu\_{\mathbf{V}}^{(k+1)}(t\_{\bar{j}}) & \text{if } r \le p\_{\varepsilon} \text{ or } j = rnbr(i), \\\mu\_{\mathbf{X}}^{(k)}(t\_{\bar{j}}) & \text{if } r > p\_{\varepsilon} \text{ and } j \ne rnbr(i), \end{cases} \tag{20}$$

where subscript *j* with *t* in *μ*(*k*) **X** (*tj*), *μ*(*k*+1) **V** (*tj*), and *μ*(*k*+1) **U** (*tj*) represent the *j*-th component of the target, mutant, and trial vectors, respectively. *r* is a random number between 0 and 1, *pc* is the crossover rate, and *rnbr*(*i*) is a randomly chosen index ∈ {1, 2, ... , *n*}, which ensures that *μ*(*k*+1) **U** (*tj*) obtains at least one component from *μ*(*k*+1) **V** (*tj*). Thereafter, all target vectors and trial vectors are combined (*μ*(*k*) **X** 8 *μ*(*k*+1) **U** ) to find the rank of the combined population using the non-dominated sorting [30] of NSGA-II. The crowding distance is also calculated for maintaining the diversity for the selection of the next generation of target vectors. The best *N* target vectors for the next generation are selected by using the environmental selection scheme of NSGA-II [30]. Multi-objective DE is terminated if the generation counter (*k*) is more than the total number of generations (*K*). Otherwise, the generation loop continues till the termination condition becomes satisfied.

### *3.3. Steps for Implementation*

In this section, the steps for implementing DE with an adaptive mutation scheme for the proposed MORBDO formulation are presented, which are as follows.

	- 3.1 Calculate the objective function values, *fm*(*μ*(*k*) **X** (*t*)).
	- 3.2 Calculate MPTP for each performance function (*i*) using Equations (15) and (16), and estimate shifting vector **S**(*k*+1) *i*,*μ***X** = *μ*(*k*) **X** − **X**(*k*) *i*,*MPTP*.
	- 3.3 Calculate the constraint violation of each performance function using the MPTP that is estimated through the chaos control theory given in Equation (15).
	- 7.1 Calculate the objective function, *fm*(*μ*(*k*+1) **U** (*t*)).

### **4. Numerical Examples**

In this section, three mathematical examples and one engineering example are solved to demonstrate the performance of the proposed method. All the examples consist of two objective functions, along with the nonlinear performance functions. The proposed method is abbreviated as SLMDE since it is developed via a single-loop method using multi-objective DE. The results of SLMDE are compared with double-loop multi-objective differential evolution (DLMDE). It is noted that PMA is used with DLMDE for reliability analysis. The reliable PO solutions are generated via both methods for different values of the target reliability index (*βT*). HV performance indicator values and number of function evaluations are used to compare the outcome. Both the methods are run 30 times with different initial populations. The standard deviation (SD) is also evaluated to see the dispersion of HV values. The Wilcoxon signed-rank test at a 5% significance level is also used to determine the difference for the statistical significance between SLMDE and DLMDE. The parameters of SLMDE and DLMDE are as follows: the scaling factor (*F*ˆ) is taken as 0.3, the crossover probability (*pc*) is 0.9, the population size (*N*) is 200, and the total number of generations (*K*) is 100 for the first example, 250 for the second example, and 200 for the car side impact example. The chaos control factor (*λ*) is considered as 0.2 [26]. The user-defined parameter (*-*) in Equation (18) is considered as 10−3. The MATLAB R2016b platform is used for developing both methods.

### *4.1. Example 1*

The first MORBDO example [3] consists of two objectives that are developed using two independent random normal variables with a standard deviation of 0.03. The example is subjected to two linear performance functions that are shown in Equation (21).

$$\begin{aligned} \min: f\_1(\mu\_\mathbf{X}) &= \mu\_{x\_1}, \\ \min: f\_2(\mu\_\mathbf{X}) &= \frac{1 + \mu\_{x\_2}}{\mu\_{x\_1}}, \\ \text{s.t.:} &P[G\_i(\mathbf{X}) > 0] \le \phi(-\beta\_i^T), \quad i = 1, 2, \\ &G\_1(\mathbf{X}) = x\_2 + 9x\_1 - 6, \\ &G\_2(\mathbf{X}) = -x\_2 + 9x\_1 - 1, \\ &0.1 \le \mu\_{x\_1} \le 1, \; 0 \le \mu\_{x\_1} \le 5. \end{aligned} \tag{21}$$

Table 1 presents the best, median, and worst values of HV obtained via SLMDE and DLMDE. SLMDE has converged to better values of HV for different *βT* values. This indicates that SLMDE generates a better set of PO solutions for the given example. It is to be noted that for a larger value of *βT*, the HV value becomes reduced, as compared to the lower *βT* value. This is because a larger value of *βT* signifies a high degree of reliability that makes the obtained PO solutions more conservative and pushes them away from the deterministic PO front inside the feasible region.

**Table 1.** Best, median, and worst HV values obtained by both methods are presented for Example 1 for different values of *βT*. The best performances are highlighted in bold font.

Figure 2 demonstrates the PO solutions obtained by both methods for different values of *βT*. The reliable PO solutions shown in Figure 2a,b correspond to the median HV values from Table 1. It can be seen that for larger values of *βT*, the PO solutions become conservative and move inside the feasible region. The same figure also demonstrates that some solutions coincide with the deterministic PO front that is located at the bottom right. This is because for those solutions, the target reliability is satisfied for the performance function *G*1(**x**).

**Figure 2.** The obtained PO solutions by both methods for example 1 for different *βT* values.

The computational efficiencies of both methods are measured with the help of a number of function evaluations that are presented in Table 2. It can be seen that the proposed method requires 202,000 function evaluations, which is only 14.85% of DLMDE. This is because SLMDE is based on a single-loop method, where the reliability of the performance function is estimated using KKT optimality conditions. On the other hand, DLMDE performs PMA for reliability estimation, which requires many function evaluations. Since the number of iterations for PMA is kept fixed, the number of function evaluations is the same for DLMDE with different values of *βT*.

**Table 2.** Number of function evaluations required by both methods for example 1.

The Wilcoxon test results are shown in the same table with symbols (+, =, −). The symbol '+' suggests a significantly better performance of SLMDE over DLMDE. Other symbols '−' and '=' suggest a significantly bad performance and an equivalent performance of SLMDE over DLMDE, respectively. It can be seen from the table that SLMDE shows a significantly better performance over DLMDE.

The progress of HV and heuristic convergence parameter *ζ* with respect to iterations are shown in Figure 3. It can be seen that there are some initial fluctuations in both HV and *ζ*, which subsidise after 10 generations and stabilize after 50 generations.

**Figure 3.** Progress of hypervolume and *ζ* of SLMDE with respect to number of generations for example 1.

### *4.2. Example 2*

The second example [51] consists of two objective functions which are highly nonlinear. It has four linear and two nonlinear performance functions that are developed using two independent random normal variables, each with a standard deviation of 0.3. The RBDO formulation of this example is given in Equation (22).

$$\begin{aligned} \text{minim: } f\_1(\mu\_\mathbf{X}) &= -[25(\mu\_{\mathbf{x}1}-2)^2 + (\mu\_{\mathbf{x}2}-2)^2 + (\mu\_{\mathbf{x}3}-1)^2 + (\mu\_{\mathbf{x}4}-4)^2 + (\mu\_{\mathbf{x}5}-1)^2], \\ \text{minim: } f\_2(\mu\_\mathbf{x}) &= [\mu\_{\mathbf{x}1}^2 + \mu\_{\mathbf{x}2}^2 + \mu\_{\mathbf{x}3}^2 + \mu\_{\mathbf{x}4}^2 + \mu\_{\mathbf{x}5}^2 + \mu\_{\mathbf{x}6}^2], \\ \text{s.t.: } P[G\_i(\mathbf{X}) > 0] &\le \phi(-\theta\_i^T), \quad i = 1, \dots, 6 \\ G\_1(\mathbf{X}) &= \mathbf{x}\_1 + \mathbf{x}\_2 - 2, \\ G\_2(\mathbf{X}) &= 6 - \mathbf{x}\_1 - \mathbf{x}\_2, \\ G\_3(\mathbf{X}) &= 2 - \mathbf{x}\_2 + \mathbf{x}\_1, \\ G\_4(\mathbf{X}) &= 2 - \mathbf{x}\_1 + 3\mathbf{x}\_2, \\ G\_5(\mathbf{X}) &= 4 - (\mathbf{x}\_3 - 3)^2 - \mathbf{x}\_4, \\ G\_6(\mathbf{X}) &= (\mathbf{x}\_5 - 3)^2 + \mathbf{x}\_6 - 4, \\ 0 &\le \mu\_{\mathbf{x}1}, \mu\_{\mathbf{x}2}, \mu\_{\mathbf{x}3} \le 10, \; 1 \le \mu\_{\mathbf{x}3}, \mu\_{\mathbf{x}5} \le 5, \ 0 \le \mu\_{\mathbf{x}4} \le 6. \end{aligned}$$

Table 3 presents the statistical values of HV obtained via both methods. In can be seen that SLMDE has converged to better values of HV for different *βT* values. This indicates that SLMDE generates a better set of PO solutions for this given example. In this case, a similar observation can also be made where for larger values of *βT*, the HV values becomes reduced. This is due to the fact that larger values of *βT* signify a larger degree of reliability, which leads to the generation of conservative PO solutions. The Wilcoxon test results are shown in the same table with symbols (+, =, −). It can be seen from the table that SLMDE shows a significantly better performance over DLMDE.

Figure 4 shows the reliable PO solutions generated in the run, corresponding to a median HV value from Table 3. It can be seen that for larger *βT*, PO solutions move inside the feasible region and away from the deterministic PO front. The spread of solutions is less in the case of SLMDE for *βT* = 1.0. The solutions are nicely distributed in the case of DLMDE. The shift of the solutions is more for larger values of *βT*, which leads to smaller values of HV that can be seen from Table 3.

**Table 3.** Best, median, and worst HV values obtained by both methods are presented for Example 2 for different values of *βT*. The best performances are highlighted in bold font.

**Figure 4.** The PO solutions obtained via both methods for example 2, for different *βT* values.

Table 4 presents the computational efficiency of both methods. The proposed method only requires 3,915,600 function evaluations, which is only 3.5–2.7% of DLMDE. It suggests that DLMDE needs many function evaluations because PMA is performed for reliability estimation.

**Table 4.** Number of function evaluations required by both methods for Example 2.

Figure 5 shows the progress of HV and *ζ* with respect to the number of generations. It can be seen that there are fluctuations for all values of *βT* until the termination criterion is achieved. The initial fluctuations can also be observed for *ζ*, which subsidize after 150 generations.

**Figure 5.** Progress of hypervolume and *ζ* of SLMDE with respect to the number of generations for example 2.

### *4.3. Car Side Impact Example*

The car side impact [3] example is considered as an engineering RBDO example, which is formulated by using 2 objectives and 10 performance functions. It consists of 11 random design variables that are normally distributed and that are grouped into random variables (*x*1, ... , *x*7) and random parameters (*x*8, ... , *x*11). The details of the variables with their standard deviation values are given in Table 5. The RBDO formulation is presented in Equation (23). The mathematical expressions for each function are given in Table 6.

$$\begin{aligned} \text{minim: } &f\_1(\mathbf{x}\_{\text{|}}) \equiv \text{Structure risk}, \\ \text{min: } &f\_2(\mathbf{x}\_{\text{|}}) \equiv \text{Average risk deflection}, \\ \text{s.t.: } &P[G\_1(\mathbf{X}) > 0] \le \phi(-\beta\_1^T), \quad i = 1, \ldots, 10 \\ &G\_1(\mathbf{X}) = \text{Abdomen load} \le 1 \text{KN}, \\ &G\_2(\mathbf{X}) = V \times \text{Cupper} \le 0.32 \text{m/s}, \\ &G\_3(\mathbf{X}) = V \times \text{Cupner} \le 0.32 \text{m/s}, \\ &G\_4(\mathbf{X}) = V \times \text{Cuper} \le 0.32 \text{m/s}, \\ &G\_5(\mathbf{X}) = \text{Upper risk deflection} \le 32 \text{mm}, \\ &G\_6(\mathbf{X}) = \text{Midderib deflection} \le 32 \text{mm}, \\ &G\_7(\mathbf{X}) = \text{Lower risk deflection} \le 32 \text{mm}, \\ &G\_8(\mathbf{X}) = \text{Public force} \le 48 \text{K}, \\ &G\_9(\mathbf{X}) = \text{Velocity of V-Pillar} \le 9.9 \text{mm/ms}, \\ &G\_{10}(\mathbf{X}) = \text{Front door velocity of V-Pillar} \le 15.7 \text{mm/s}, \\ &0.5 \le \mu\_{x\_1}, \mu\_{y\_1}, \mu\_{x\_4} \le 1.5, 0.45 \le \mu\_{x\_2} \le 1.35, 0.875 \le \mu\_{x\_3} \le 2.625, \\ &0.4 \le \mu\_{y\_2} \le 1.2, 0.4 \le \mu\_{y\_3} \le 1.2, 0.192 \le \mu\_{y\_4}, \mu\_{y\_5} \le 0.75. \end{aligned}$$

The statistical values of the HV values obtained from both methods with respect to different *βT* are presented in Table 7. The proposed method has converged to larger values of HV for all *βT* values. This signifies a better distribution of PO solutions of SLMDE as compared to DLMDE. The observation of reducing HV values with larger *βT* values remains the same. The Wilcoxon test results are shown in the same table with symbols (+, =, −). It can be seen from the table that SLMDE shows significantly better, bad, and equivalent performances over DLMDE for *βT* = 1, 2 and 3, respectively.

The obtained reliable PO solutions for both methods are shown in Figure 6. As observed with previous examples, for larger values of *βT*, the PO solutions start moving away from the deterministic PO front inside the feasible region.

Table 8 presents the computational efficiency of both methods. In this example, SLMDE requires only 4,623,000, the number of function evaluations, which is only 0.3–0.26% that of DLMDE. Since SLMDE performs an approximate reliability estimation by using the KKT optimality conditions, it saves many function evaluations compared to DLMDE. Figure 7 shows a similar progress for HV and *ζ* with respect to the number of generations. There are initial fluctuations for all values of *βT*, which subside after 80 generations.

**Table 5.** Details of design variables and their standard deviation values.

**Table 6.** The objectives and performance functions of Example 3.

**Table 7.** Best, median, and worst HV values obtained via both methods, presented for car side impact example for different values of *βT*. The best performances are highlighted in bold font.

**Table 8.** Number of function evaluations required by both methods for car side impact example.

**Figure 6.** The obtained PO solutions by both methods for car side impact example for different *βT* values.

**Figure 7.** Progress of hypervolume and *ζ* of SLMDE with respect to number of generations for car side impact problem.

### *4.4. Example 4*

The fourth example [44] consists of two objective functions and both of them are quadratic functions. The example has a linear performance function developed with three independent random normal variables. The RBDO formulation is given in Equation (24).

$$\begin{aligned} \text{min:} & f\_1(\mu\_\mathbf{X}) = (\mu\_{x\_1} - 1)^2 + (\mu\_{x\_2} - 2)^2 + (\mu\_{x\_3} - 3)^2, \\ \text{min:} & f\_2(\mu\_\mathbf{X}) = \mu\_{x\_1}^2 + 2\mu\_{x\_2}^2 + 3\mu\_{x\_3}^2, \\ \text{s.t.:} & P[G\_i(\mathbf{X}) > 0] \le \phi(-\beta\_i^T), \quad i = 1, \\ & G\_1(\mathbf{X}) = \mathbf{x}\_1 + \mathbf{x}\_2 + \mathbf{x}\_3 - 1, \\ & 0.1 \le \mu\_{x\_i} \le 6, \ i = 1, 2, 3. \end{aligned} \tag{24}$$

where *x*1 ∼ *N*(1, 0.05), *x*2 ∼ *N*(2, 0.1), and *x*3 ∼ *N*(3, 0.15).

Table 9 presents the statistical values of HV obtained via SLMDE and DLMDE. It can be observed that in most of the cases, SLMDE converged to better values of HV for different *βT*. The HV values become reduced with larger values of *βT*. The Wilcoxon test results are shown in the same table with symbols (+, =, −). It can be seen from the table that SLMDE shows an equivalent performance with DLMDE for *βT* = 2 and 3, and a bad performance for *βT* = 1.

**Table 9.** Best, median, and worst HV values obtained by both methods are presented for example 4 for different values of *βT*. The best performances are highlighted in bold font.

Figure 8 shows the reliable PO solutions generated in the run corresponding to the median HV value obtained via both methods for different values of *βT*. As observed in the previous examples, for larger *βT*, PO solutions move inside the feasible region, away from the deterministic PO front.

**Figure 8.** The PO solutions obtained via both methods for example 4 for different *βT* values.

Table 10 presents the computational efficiencies of both methods. The proposed method requires only 50% of function evaluations as that of DLMDE. Figure 9 also shows similar observations for HV and *ζ* during the progress of the generations. There is an initial fluctuation which reduces after 20 generations.

**Table 10.** Number of function evaluations required by both methods for Example 4.

**Figure 9.** Progress of hypervolume and *ζ* of SLMDE with respect to number of generations for Example 4.

### **5. Conclusions**

A single-loop multi-objective reliability-based design optimization has been proposed for generating reliable PO solutions quickly. It was developed by applying KKT optimality conditions to PMA for generating an approximate expression of MPTP. The search direction of approximate MPTP was modified via chaos control theory. The concept of the shifting vector approach was implemented with the novel formulation to include both target and trial vectors. DE was made adaptive, using the heuristic parameter that helped DE to perform different mutation operators. The proposed SLMDE was tested on three mathematical and one engineering bi-objective RBDO examples. It was found that SLMDE generated more reliable PO solutions for all examples compared to DLMDE. The results demonstrate that the convergence of SLMDE takes less function evaluations than DLMDE. For all four examples, the SLMDE was able to generate better HV values. For example 2, a lot of fluctuations during the progress of hypervolume can be observed, which stabilize gradually. The user-defined parameter *ζ* shows stable progress for all the examples. In the future, the proposed method can be modified for quick convergence by incorporating quantile approximation for reliability analysis. The proposed method can also be tested on other real-world examples having many nonlinear functions.

**Author Contributions:** Conceptualization: R.B. and D.S.; methodology: R.B.; validation: R.B.; formal analysis: R.B. and D.S.; investigation: R.B.; writing—original draft preparation: R.B.; writing—review and editing: R.B. and D.S.; visualization: R.B.; supervision: D.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:

### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Mathematical and Computational Applications* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/mca

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-6981-9


**3. Findings and discussions**

#### **3.1 Findings concerning the first sub-problem**

The first sub-problem of the study is "What are the pre-service science teachers' attitudes towards educational technologies?" The findings concerning this sub-problem reveal that the mean score for the pre-service science teachers' attitudes towards educational technologies is 169.66, standard deviation is 14.19; the lowest score taken from the attitude scale is 125 and the highest score is 207. According to these scores, the pre-service science teachers' general attitude is in the category of "Agree". This finding shows that the students in general have positive attitudes towards educational technologies. There are similar findings reported in the literature. Gunter, Gunter & Wiens (1998) found that pre-service teachers have more positive attitudes towards working on computer and learning through computer and technology in general. Ylmaz (2005), in his thesis study, investigated the effects of technology on students' achievement and attitude and found positive impacts on achievement and attitude. In another study, Sevindik (2006) found positive effects of using smart classes in higher education on students' academic achievement and attitudes. Yavuz and Coşkun (2008) found that pre-service elementary school teachers have positive attitudes towards and opinions about the use of technological tools and equipments.

#### **3.2 Findings concerning the second sub-problem**

The second sub-problem of the study is "Do the pre-service science teachers' attitudes towards educational technologies vary significantly depending on the environment where they have been brought up?" T-test was conducted to test whether there is a statistically significant

The Effects of Environment and Family Factors on Pre-Service Science Teachers'

education or are illiterate.

mothers' educational status

of the test are presented in Table 5.

their mothers' educational status

teachers' fathers' educational status

Attitudes Towards Educational Technologies (The Case of Muğla University-Turkey) 449

According to Table 3, 63.4% of the pre-service science teachers' mothers are secondary or elementary school graduates, 24.8% are high school graduates and 6.9% are graduates from graduate or undergraduate programs. 5.0% of them have no education or are illiterate. When we look at the educational status of the fathers, we can see that 49.5% are secondary school or elementary school graduates, 28.27% are high school graduates and 21.8% are graduates of undergraduate or graduate programs. There are no fathers not having

**Educational status N** X **S**  Undergraduate/graduate 7 3.86 .25

Secondary school/elementary school 64 3.96 .34

Table 4. Arithmetic means and standard deviations for the pre-service science teachers'

In Table 4, it is seen that there are differences among the arithmetic means. ANOVA test was carried out to determine whether these differences are statistically significant and the results

**Source of the variance Sum of squares df Mean of squares F p**  Between-groups .14 3 .05 .42 .74

Table 5. Anova results for the pre-service science teachers' attitude scale scores in relation to

The results of Table 5 show that there is no significant difference based on the mothers' educational status among the pre-service science teachers' attitudes towards educational technologies [F(3-97) = .42 , p>.05]. That is, there is no correlation between the pre-service science teachers' attitudes towards educational technologies and their mothers' educational

**Educational status N** X **S**  Undergraduate/graduate 22 3.93 .45

Secondary school/ elementary school 50 3.94 .30

Table 6. Arithmetic means and standard deviations concerning the pre-service science

High school 29 3.97 .29

Others - - - **Total** 101 3.95 .33

status. This finding is supported by the findings reported by Erol and Gezer (2006).

Within-groups 10.75 97 .11

**Total** 10.89 100

High school 25 3.91 .32

Others 5 4.03 .29 **Total** 101 3.95 .33

difference among the pre-service science teachers' attitudes towards technology based on the environment they have been brought up and the results are presented in Table 1 and Table 2.

Table 1. Distribution of the pre-service science teachers according to environment where they have been brought up

Of the pre-service science teachers, 15.8% come from rural area and 84.2% from urban area. As the pre-service science teachers coming from urban area probably have had better technological opportunities, they are expected to exhibit more positive attitudes.

Table 2. T-test results for the science pre-service teachers' attitude scale scores in relation to environment where they have been brought up

According to the t-test results presented in Table 2, there is no significant difference among the attitudes based on the environment where they have been brought up [t (99) = .22, p > .05]. This finding indicates that there is no significant relationship between the environment where the pre-service science teachers have been brought up and their attitudes towards educational technologies. This finding concurs with the findings of Can (2010); Erol and Gezer (2006).

#### **3.3 Findings concerning the third sub-problem**

The third sub-question of the study is "Is there a significant relationship between the preservice science teachers' attitudes towards educational technologies and their parents' educational status?" The findings concerning this problem are related to the relationship between the pre-service science teachers' attitudes and their parents' educational status. First, the distribution of the pre-service science teachers according to their parents' educational status is given in Table 3 and then the ANOVA test was carried out to determine whether there is a significant correlation between the pre-service science teachers' attitudes and their parents' educational status and then findings are presented in Table 4, Table 5, Table 6 and Table 7.

Table 3. Distribution of the pre-service science teachers' parents' educational status

difference among the pre-service science teachers' attitudes towards technology based on the environment they have been brought up and the results are presented in Table 1 and Table 2.

**Environment Frequency** % Urban area 85 84.2 Rural area 16 15.8 **Total** 101 100.0 Table 1. Distribution of the pre-service science teachers according to environment where

Of the pre-service science teachers, 15.8% come from rural area and 84.2% from urban area. As the pre-service science teachers coming from urban area probably have had better

**Environment N** X **S df t p**  Urban area 85 169.80 14.40 99 .22 .82

Table 2. T-test results for the science pre-service teachers' attitude scale scores in relation to

According to the t-test results presented in Table 2, there is no significant difference among the attitudes based on the environment where they have been brought up [t (99) = .22, p > .05]. This finding indicates that there is no significant relationship between the environment where the pre-service science teachers have been brought up and their attitudes towards educational technologies. This finding concurs with the findings of Can (2010); Erol and Gezer (2006).

The third sub-question of the study is "Is there a significant relationship between the preservice science teachers' attitudes towards educational technologies and their parents' educational status?" The findings concerning this problem are related to the relationship between the pre-service science teachers' attitudes and their parents' educational status. First, the distribution of the pre-service science teachers according to their parents' educational status is given in Table 3 and then the ANOVA test was carried out to determine whether there is a significant correlation between the pre-service science teachers' attitudes and their parents' educational status and then findings are presented in Table 4,

**Educational status Mother Father**

Table 3. Distribution of the pre-service science teachers' parents' educational status

Graduate/Undergraduate 7 6.9 22 21.8 High school 25 24.8 29 28.7 Secondary school/Elementary school 64 63.4 50 49.5 Other 5 5.0 - - **Total** 101 100.0 101 100.0

**f % f %**

technological opportunities, they are expected to exhibit more positive attitudes.

Rural area 16 168.93 13.44

environment where they have been brought up

**3.3 Findings concerning the third sub-problem**

Table 5, Table 6 and Table 7.

they have been brought up

According to Table 3, 63.4% of the pre-service science teachers' mothers are secondary or elementary school graduates, 24.8% are high school graduates and 6.9% are graduates from graduate or undergraduate programs. 5.0% of them have no education or are illiterate. When we look at the educational status of the fathers, we can see that 49.5% are secondary school or elementary school graduates, 28.27% are high school graduates and 21.8% are graduates of undergraduate or graduate programs. There are no fathers not having education or are illiterate.

Table 4. Arithmetic means and standard deviations for the pre-service science teachers' mothers' educational status

In Table 4, it is seen that there are differences among the arithmetic means. ANOVA test was carried out to determine whether these differences are statistically significant and the results of the test are presented in Table 5.

Table 5. Anova results for the pre-service science teachers' attitude scale scores in relation to their mothers' educational status

The results of Table 5 show that there is no significant difference based on the mothers' educational status among the pre-service science teachers' attitudes towards educational technologies [F(3-97) = .42 , p>.05]. That is, there is no correlation between the pre-service science teachers' attitudes towards educational technologies and their mothers' educational status. This finding is supported by the findings reported by Erol and Gezer (2006).

Table 6. Arithmetic means and standard deviations concerning the pre-service science teachers' fathers' educational status

The Effects of Environment and Family Factors on Pre-Service Science Teachers'

service science teachers' attitudes towards educational technologies.

universities.

**5. References**

234

*– TOJET,* 2 (2)

*Technology – TOJET, 4* (1), Article 12

*Bilim Toplum,* Cilt 1, Say 1, Kş Says , 51-60

Doğu Teknik Üniversitesi, 6-8 Eylül, Ankara

Technology. *Journal of Technology Education*, *13 (*2), 6-19

of Recent Research in Tojet. *TOJET*, 8(2).

Thousand Oaks, pp. 1-23, California

of such technologies in their teaching.

Attitudes Towards Educational Technologies (The Case of Muğla University-Turkey) 451

given in their former education and in this way they may have developed more positive attitudes towards educational technologies. Another finding of the present study is that there is no significant correlation between the parents' educational status and the pre-

But before making some generalizations in light of the findings of the present study, the limitations of the study should be mentioned. First, the present study is limited to its study group and data collection tools used in the present study. Therefore, further research may look at students from different departments, different faculties or different

Utilization of educational technologies in the field of education can enrich education and enhance students' motivation, in this way; students are promoted to develop positive attitudes towards educational technologies. Positive attitudes developed by pre-service science teachers towards educational technologies may help them to make more efficient use

Akdeniz, A.R. & Yiğit, N. (2001). Fen Bilimleri Öğretiminde Bilgisayar Destekli

Akpnar, Y. (2003). Öğretmenlerin Yeni Bilgi Teknolojileri Kullanmnda Yükseköğretimin

Akpnar, Y. (2004). Eğitim Teknolojisiyle İlgili Öğrenmeyi Etkileyebilecek Baz Etmenlere

Akpnar, E., Aktamş, H., & Ergin, Ö. (2005). Fen Bilgisi Dersinde Eğitim Teknolojisi

Aksoy, H. H. (2003). Uluslararas Karşlaştrma Ölçütlerinin Kullanm ve Türkiye. *Eğitim*

Aktamş, H., Akpnar, E., & Ergin, Ö. (2002). *Yapsalc Kurama Örnek Bir Uygulama*. V. Ulusal

Alper, A. & Gülbahar, Y. (2009). Trends and Issues in Educational Technologies: A Review

Becker, K. H. & Maunsaiyat, S. (2002). Thai Students' Attitudes And Concepts of

Berne, R. (1994). Educational Input and Outcome Inequities in New York State. *Outcome*

Binbaşoğlu, C. (1994). *Genel Öğretim Bilgisi*. Yedinci Basm, Kadoğlu Matbaas, Ankara.

*TOJET* July ISSN: 1303-6521 Volume 3, Issue 3, Article 15. Sakarya

Materyallerin Öğrenci Başars Üzerine Etkisi. Yeni Bin yln Başnda Türkiye'de Fen Bilimleri Eğitimi Sempozyumu Bildirileri, Maltepe Üniversitesi, İstanbul, 229-

Etkisi: İstanbul Okullar Örneği. *The Turkish Online Journal of Educational Technology*

Karş Öğretmen Yaklaşmlar. *The Turkish Online Journal of Educational Technology –*

Kullanlmasna İlişkin Öğrenci Görüşleri. *The Turkish Online Journal of Educational*

Fen Bilimleri ve Matematik Eğitim Kongresi'nde bildiri olarak sunulmuştur. Orta

*Equity in Education*, 15th Annual Yearbook of the American Education Finance Association. Editors: Robert Berne ve Lawrance O. Picus. Corwin Press, Inc.

Variance analysis was conducted to see whether the differences seen among arithmetic means in Table 6 are significant, and the results are presented in Table 7.

Table 7. Anova results for the pre-service science teachers' attitude scale scores in relation to their fathers'' educational status

Anova results presented in Table 7 show that there is no significant difference based on the fathers' educational status among the pre-service science teachers' attitudes towards educational technologies [F(2-98) = .07, p>.05]. That is, there is no correlation between the preservice science teachers' attitudes towards educational technologies and their fathers' educational status. This finding is in compliance with the findings of Erol and Gezer (2006).


Sustainable Development Authoritative and Leading Edge Content for Environmental Management

*Edited by Sime Curkovic*

**SUSTAINABLE**

**DEVELOPMENT –**

**MANAGEMENT**

Edited by **Sime Curkovic**

**AUTHORITATIVE AND**

**FOR ENVIRONMENTAL**

**LEADING EDGE CONTENT**


*3.3. Monte Carlo Simulation and Visual Minteq*

Due to individual differences and sampling limitations, the results of the health risk calculation were highly uncertain. In order to overcome this defect, this study used the Monte Carlo Simulation (MCS) method, which is a probabilistic statistical mathematical

method to evaluate uncertainty by random sampling each variable value [40]. This simulation method was combined with the U.S. Environmental Protection Agency risk assessment technology to quantify and minimize the uncertainty, and analyze the sensitivity of pollutant exposure in drinking water wells. Different numbers related to each variable were sampled for a large number of digital repetition operations. The main steps are as follows: (1) The respective distribution functions were fitted to the data of each variable, and the hypothetical variable was defined; (2) The calculation formula was inputted and the decision variable was defined; (3) A random sampling was conducted from the distribution of the hypothetical variables; (4) The randomly selected parameter sequence was used for a large number of repeated operations, to obtain as output the probability distribution of the operation results (Figure 3). In this study, the Monte Carlo Model with 10,000 iterations was used to evaluate the carcinogenic risk of exposure to As and Cr6+ for children and adults during the rainy and dry seasons from 2017 to 2018. Then, sensitivity analysis was carried out using the MCS method with 10,000 repetitions using Oracle Crystal Ball®. Finally, the frequency diagram and sensitivity analysis diagram drawn by it were used to identify the input parameters with the greatest impact on the output of the risk assessment model. After using the sensitivity analysis to obtain the main influencing factors of health risk assessment results, Visual Minteq was used to analyze the possible ion forms of the main factors in the hydrogeological environment of the study area. *Water* **2022**, *14*, x FOR PEER REVIEW 10 of 23 analyze the possible ion forms of the main factors in the hydrogeological environment of the study area.

**Figure 3.** Monte Carlo model operation flowchart.

**Figure 3.** Monte Carlo model operation flowchart.

4.1.2. Non‐Carcinogenic Risk Assessment

**4. Results and Discussion** *4.1. Health Risk Assessment* 4.1.1. Health Risk Assessment

risk of eight ions and the non‐carcinogenic risk of two heavy metal ions [42].

The concentrations of ions in water were used to assess the possible health risks from human exposure to drinking water through both oral intake and skin contact. This study considered two groups, namely adults and children [41], and evaluated the carcinogenic

The Hazard Index (HI) was used to measure non‐carcinogenic health risk as the sum of the hazard quotients (HQs) from both oral intake and dermal contact of several parameters including F−, NH4+, NO2−, NO3−, Cr6+, Mn, As, Fe (including Fe2+, Fe3+), and


**6. Summary**

VFTO [p.u.]

In this work, a higher damping effect was investigated by designing a spiral litz coil connected in series to the damping busbar, which increased the total inductance value. The simulation results show that there was a considerable increase in the damping effect, as shown in Figure 10. Furthermore, improving the design and material type could be useful to achieving an optimal suppressing effect. The main advantages of this improved design are the higher damping effect, accuracy, simplicity, and ability to be used in any existing GIS without any change in the GIS design or construction. The transient characteristics of the VFTO after adding a 0.1 mH litz coil are illustrated in Table 6.

**Table 6.** The transient characteristics of VFTO under an opening operation at load side with and without installing the spiral litz coil (0.1 mH).

Another obvious advantage of this improved design is no magnetic saturation problem and mitigation of the overvoltage magnitude due to dissipation of the transient energy in the damping resistance. On the other hand, litz wire has been classified as a cost-effective choice to reduce eddy current loss in high-frequency transformer and inductor windings. Some solutions could be successful in increasing the dielectric strength for a high voltage level or for safety requirements, like adding thermoplastic insulation such as PTFE, which could be used for bundle insulation on litz wire in some applications [36].

The proposed improved damping busbar has limitations when thermoplastic insulation is added due to the extra cost. So, in future work, it will be fundamental to develop the design to improve the estimation quality to meet the industry's requirements for low costs. However, compared with other existing expensive VFTO suppression methods that require changes in GIS design and construction, the proposed method can be easily implemented in any existing GIS without any change in its structure, making it very cost-effective. Therefore, it is necessary to use simulation calculations to optimize the insulation reliability of the uniform electric field and increase the insulation of the busbar. Table 7 shows the mitigation effect of different damping methods, including the new method.

**Table 7.** Comparison of different mitigation methods with the new damping busbar.


The EI has a close relationship with *mp*, i.e., EI = *mp*/(π/4). A linear regression method was used to estimate the *mp* value based on the number of leaves. However, EI can also be calculated directly using *A*, *L* and *W*. The EI quantifies the extent of a leaf shape deviating from an ellipse regardless of the eccentricity of the ellipse. Thus, we used the *W/L* ratio as another leaf-shape index. ANOVA with Tukey's honestly significant difference (HSD) test [31] and a 0.05 significance level were used to test the significance of the differences in EI as well as the *W/L* ratio between the two taxa at the combined data level and at the intraspecific level across different dates of collection. For any two groups, we calculated the 95% confidence interval (CI) of the Tukey's test confidence interval, which equals the difference in the means between the two groups ± HSD, and we determined its significance by checking whether the 95% CI included 0. We also used the linear and parabolic equations

to fit the EI (and the *W/L* ratio) vs. collection month data, where month was set to a numeric variable to test whether leaf shape has a monotonic or parabolic trend across different investigation dates (i.e., leaf ages). to be fitted (referred to as regression coefficients). For the simple linear regression, 2 β was set to be zero. Each regression coefficient's statistic in Equation (5) follows a *t* distribution with *n* − *p* degrees of freedom, where *n* is the sample size and *p* is the number of

*Plants* **2022**, *11*, x FOR PEER REVIEW 4 of 9

different investigation dates (i.e., leaf ages).

$$\text{LS} = \beta\_0 + \beta\_1 \text{Month} + \beta\_2 \text{Month}^2,\tag{5}$$

where LS is the leaf-shape index of interest (EI or the *W/L* ratio); Month is the month of sampling (ranging from 3 to 8 [i.e., March to August 2021] for *P. × fraseri*, and from 4 to 11 [i.e., April to November 2021] for *P. serratifolia*); and 0 β , 1 β and 2 β are the parameters

as well as the *W/L* ratio between the two taxa at the combined data level and at the intraspecific level across different dates of collection. For any two groups, we calculated the 95% confidence interval (CI) of the Tukey's test confidence interval, which equals the difference in the means between the two groups ± HSD, and we determined its significance by checking whether the 95% CI included 0. We also used the linear and parabolic equations to fit the EI (and the *W/L* ratio) vs. collection month data, where month was set to a numeric variable to test whether leaf shape has a monotonic or parabolic trend across

> 2 01 2 LS =+ + β β Month β Month *,* (5)

where LS is the leaf-shape index of interest (EI or the *W/L* ratio); Month is the month of sampling (ranging from 3 to 8 [i.e., March to August 2021] for *P.* × *fraseri*, and from 4 to 11 [i.e., April to November 2021] for *P. serratifolia*); and β0 , β1 and β2 are the parameters to be fitted (referred to as regression coefficients). For the simple linear regression, β2 was set to be zero. Each regression coefficient's statistic in Equation (5) follows a *t* distribution with *n* − *p* degrees of freedom, where *n* is the sample size and *p* is the number of parameters [32]. The statistical software R (version 4.2.0) [28] was used to carry out calculations and to draw figures. The 'agricolae' package (version 1.3-5) was used to implement the Tukey's HSD test. **3. Results**

The statistical software R (version 4.2.0) [28] was used to carry out calculations and to draw figures. The 'agricolae' package (version 1.3-5) was used to implement the Tukey's HSD test. The Montgomery equation (ME) was valid in calculating leaf area. For the pooled intraspecific data, the correlation coefficient was greater than 0.99 (*p* < 0.001). The esti-

### **3. Results** mated values of the Montgomery parameters of *P. × fraseri* and *P. serratifolia* were 0.6494

coefficient.

parameters [32].

The Montgomery equation (ME) was valid in calculating leaf area. For the pooled intraspecific data, the correlation coefficient was greater than 0.99 (*p* < 0.001). The estimated values of the Montgomery parameters of *P.* × *fraseri* and *P. serratifolia* were 0.6494 (95% CI = 0.6482, 0.6506) and 0.6981 (95% CI = 0.6971, 0.6990),respectively (Figure 2). There was a significant difference in leaf shape between the two taxa, i.e., the leaves of *P. serratifolia* were on average wider and more elliptical in shape (Figure 3). The means of the *W/L* ratios for the two taxa were 0.38 (*P.* × *fraseri*) and 0.40 (*P. serratifolia*); the means of the leaf elliptical indices (EIs) for the two taxa were 0.8276 (*P.* × *fraseri*) and 0.8894 (*P. serratifolia*). (95% CI = 0.6482, 0.6506) and 0.6981 (95% CI = 0.6971, 0.6990), respectively (Figure 2). There was a significant difference in leaf shape between the two taxa, i.e., the leaves of *P. serratifolia* were on average wider and more elliptical in shape (Figure 3). The means of the *W/L* ratios for the two taxa were 0.38 (*P. × fraseri*) and 0.40 (*P. serratifolia*); the means of the leaf elliptical indices (EIs) for the two taxa were 0.8276 (*P. × fraseri*) and 0.8894 (*P. serratifolia*).

**Figure 2.** Log-log bivariate plot of leaf area (*A*) vs. the product leaf length (*L*) and width (*W*). The symbol ෝ represents the estimated value of the Montgomery parameter, i.e., the estimated proportionality coefficient of the Montgomery equation. The subscripts Pf and Ps represent the hybrid **Figure 2.** Log-log bivariate plot of leaf area (*A*) vs. the product leaf length (*L*) and width (*W*). The symbol *m*ˆ *p* represents the estimated value of the Montgomery parameter, i.e., the estimated proportionality coefficient of the Montgomery equation. The subscripts Pf and Ps represent the hybrid *P.* × *fraseri*, and one of its parents *P. serratifolia*, respectively. \*\*\* denotes *p* < 0.001 for a correlation coefficient.

*P. × fraseri*, and one of its parents *P. serratifolia*, respectively. \*\*\* denotes *p* < 0.001 for a correlation

**Figure 3.** Comparison of leaf shape between the two *Photinia* taxa for the pooled data. (**A**) Ratio of leaf width to length; (**B**) Leaf ellipticalness index. In each panel, the letters on the whiskers show differences between the two taxa (i.e., taxa with different letters are significantly different at the α = 0.05 significance level using Tukey's HSD); the values at the top of whiskers represent the coefficients of variation (%) for each taxon; the horizontal solid line represents the median; and the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. The difference in mean leaf *W/L* ratios is equal to 0.0193 with HSD = 0.0017, and the corresponding 95% Tukey's test confidence interval is 0.0176, 0.0210, which does not include 0, indicating a significant difference; the difference in mean leaf ellipticalness indices is equal to 0.0617 with HSD = 0.0019, and the corresponding 95% Tukey's test **Figure 3.** Comparison of leaf shape between the two *Photinia* taxa for the pooled data. (**A**) Ratio of leaf width to length; (**B**) Leaf ellipticalness index. In each panel, the letters on the whiskers show differences between the two taxa (i.e., taxa with different letters are significantly different at the α = 0.05 significance level using Tukey's HSD); the values at the top of whiskers represent the coefficients of variation (%) for each taxon; the horizontal solid line represents the median; and the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. The difference in mean leaf *W/L* ratios is equal to 0.0193 with HSD = 0.0017, and the corresponding 95% Tukey's test confidence interval is 0.0176, 0.0210, which does not include 0, indicating a significant difference; the difference in mean leaf ellipticalness indices is equal to 0.0617 with HSD = 0.0019, and the corresponding 95% Tukey's test confidence interval is 0.0598, 0.0636, which also does not include 0, indicating a significant difference.

confidence interval is 0.0598, 0.0636, which also does not include 0, indicating a significant difference. There were significant differences in leaf shape (measured as either the *W/L* ratio or the EI) among the different months of collecting (Figure 4). Three of the four slopes of the linear model of leaf shape vs. month were statistically significant (with three *p* values < 0.05) with the exception of the EIs of *P. serratifolia*. However, the four coefficients of deter-There were significant differences in leaf shape (measured as either the *W/L* ratio or the EI) among the different months of collecting (Figure 4). Three of the four slopes of the linear model of leaf shape vs. month were statistically significant (with three *p* values < 0.05) with the exception of the EIs of *P. serratifolia*. However, the four coefficients of determination (i.e., *r* 2 values) were all less than 0.06, i.e., a very weak linear relationship was observed between leaf shape and sampling month. Thus, there is no robust evidence for a monotonic increase or decrease in leaf shape across leaf ages. The linear and quadratic coefficients in the parabolic model were statistically significant (*p* values < 0.05) for the leaf *W/L*

mination (i.e., *r*2 values) were all less than 0.06, i.e., a very weak linear relationship was observed between leaf shape and sampling month. Thus, there is no robust evidence for a monotonic increase or decrease in leaf shape across leaf ages. The linear and quadratic

leaf *W/L* ratio vs. month relationship and the EI vs. month relationship for *P. × fraseri*. However, the two coefficients were not statistically significant (*p* values > 0.05) for any leaf-shape indices for *P. serratifolia*. Three of the four coefficients of determination of the parabolic regression were smaller than 0.03, with the exception of the EI vs. month data of *P. × fraseri* (*r*2 = 0.1452). Thus, there was also no strong evidence for a parabolic trend in

leaf shape across the different leaf ages (see Table 1 for details).

*Photinia × fraseri*  Leaf width/length ratio

*Photinia × fraseri*  Leaf ellipticalness index

ratio vs. month relationship and the EI vs. month relationship for *P.* × *fraseri*. However, the two coefficients were not statistically significant (*p* values > 0.05) for any leaf-shape indices for *P. serratifolia*. Three of the four coefficients of determination of the parabolic regression were smaller than 0.03, with the exception of the EI vs. month data of *P.* × *fraseri* (*r* 2 = 0.1452). Thus, there was also no strong evidence for a parabolic trend in leaf shape across the different leaf ages (see Table 1 for details). *Plants* **2022**, *11*, x FOR PEER REVIEW 6 of 9

**Figure 4.** Comparison of leaf shapes (reflected by the ratio of leaf width to length and leaf ellipticalness index) across different sampling months for *P. × fraseri* (**A**,**B**) and *P. serratifolia* (**C**,**D**). In each panel, the letters on the whiskers represent the significance of the differences in leaf shape between any two sampling months among which the letters a, b, c, and d are used to represent the significance of differences in leaf shape among different sampling months; groups sharing a common letter are not significantly different at the 0.05 significance level; the values at the top of whiskers represent the coefficients of variation (%) for each sampling month; the horizontal solid line represents the median; the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. For the label of the *x*-axis, the months 3 to 11 represent March to November 2021, respectively. The differences in the means of leaf-shape indices, HSD values, and the corresponding 95% Tukey's test confidence intervals between any two groups for *P. × fraseri* and *P. serratifolia* are provided in Tables S3 and S4, respectively. **Figure 4.** Comparison of leaf shapes (reflected by the ratio of leaf width to length and leaf ellipticalness index) across different sampling months for *P.* × *fraseri* (**A**,**B**) and *P. serratifolia* (**C**,**D**). In each panel, the letters on the whiskers represent the significance of the differences in leaf shape between any two sampling months among which the letters a, b, c, and d are used to represent the significance of differences in leaf shape among different sampling months; groups sharing a common letter are not significantly different at the 0.05 significance level; the values at the top of whiskers represent the coefficients of variation (%) for each sampling month; the horizontal solid line represents the median; the red asterisk represents the mean. The whiskers extend to the most extreme data point, which is no more than 1.5 times the interquartile range from the box. For the label of the *x*-axis, the months 3 to 11 represent March to November 2021, respectively. The differences in the means of leaf-shape indices, HSD values, and the corresponding 95% Tukey's test confidence intervals between any two groups for *P.* × *fraseri* and *P. serratifolia* are provided in Tables S3 and S4, respectively.

*a* 0.3733 <0.001

*c* −0.0006 0.0031

*y = a + bx a*0.8549 <0.001 0.0591 *b* −0.0050 <0.001

*y = a + bx + cx*2 *a*0.9668 <0.001 0.1452 *b* 0.0504 <0.001

*b* 0.0054 0.0365 0.0233

**Taxon/Leaf-Shape Index Equation † Parameter Estimate Significance** *r***2**

*y = a + bx + cx*2

**Table 1.** Regression statistics of the leaf-shape index vs. sampling month.

**Table 1.** Regression statistics of the leaf-shape index vs. sampling month.

† Here, *y* represents a specific leaf-shape index, and *x* represents the sampling month.

### **4. Discussion and Conclusions**

Shi et al. [33] used the Montgomery parameter (*mp*) to fit leaf *A* vs. *LW* of 101 bamboo taxa using > 10,000 bamboo leaves. In their study, the estimated value of *mp* was 0.6959 (95% CI = 0.6952, 0.6966). Li et al. [23] estimated a *mp* value of 0.6840 (95% CI = 0.6827, 0.9855) using > 2200 leaves from nine Magnoliaceae species, whereas Ma et al. [13] estimated a *mp* value of 0.7710 (95% CI = 0.7696, 0.7714) using > 6500 leaves of an alpine oak species (*Quercus pannosa*) from different tree size groups. The present study shows that the estimated value of *mp* of the two *Photinia* taxa also approximates 0.7, i.e., the leaf area approximately accounts for 70% of the area of a hypothetical rectangle with the leaf length and width as its two side lengths. This result is significantly different from the numerical value of the alpine oak provided in ref. [13], because the leaf shape of *Q. pannosa* is a special superellipse with a mean leaf ellipticalness index greater than unity [34–36], which is larger than those of the two *Photinia* taxa examined in this study. Since leaf shape is found to be closely related to climatic factors [18,37], it would be worth studying the link between the *mp* and climate for some widely distributed plants in the future.

The leaf shapes of two *Photinia* taxa were quantified using different metrics, and found to significantly differ in spite of the fact that the two types of leaves appear to be superficially very similar in appearance. The leaves of *P. serratifolia* are wider and more elliptical than those of its hybrid. In addition, leaf shape significantly varied across the different leaf ages, but manifested no monotonic or parabolic tendency across the different leaf ages. The variations in leaf shape among the different months of sampling may reflect variations resulting from random samplings of individual plants. This variation, however, does not mask the fact that the leaf shape of the *Photinia* taxa is relatively stable, and that the Montgomery equation, which assumes a proportional relationship between leaf area and the product of leaf length and width, is valid at different leaf growth stages. These results highlight an approach that permits the non-destructive estimation of leaf area, i.e., direct field measurements of leaf length and width are used to estimate lamina area by multiplying their product by a proportionality coefficient (*mp*).

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants11182370/s1, Table S1: The raw leaf data of *Photinia* × *fraseri*; Table S2: The raw leaf data of *P. serratifolia*; Table S3: Tukey multiple comparisons of the means of leafshape indices between any sampling months for *P.* × *fraseri*; Table S4: Tukey multiple comparisons of the means of leaf-shape indices between any sampling months for *P. serratifolia*.

**Author Contributions:** Methodology, K.J.N., D.A.R. and P.S.; formal analysis, K.J.N., D.A.R. and P.S.; investigation, X.Z., K.J.N., D.A.R., Y.J., H.D. and P.S.; writing—original draft preparation, X.Z. and P.S.; writing—review and editing, K.J.N. and D.A.R.; funding acquisition, X.Z. and H.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** X.Z. and H.D. were funded by the Project of Biological Resources Survey in Wuyishan National Park (grant number: HXQT2020120701).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in the present work have been listed in the online Supplementary Materials.

**Acknowledgments:** We thank Kexin Yu for her valuable help during the preparation of this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**