Multiobjective Wind Driven Optimization Approach ... - IEEE Xplore

Multiobjective Wind Driven Optimization Approach Applied to Transformer Design Helon Vicente Hultmann Ayala

Emerson Hochsteiner de Vasconcelos Segundo

Pós-Graduação em Engenharia de Produção e Sistemas, PPGEPS, Pontifícia Universidade Católica do Paraná, PUCPR, Curitiba, Brazil [email protected]

Pós-Graduação em Engenharia Mecânica, PPGEM, Pontifícia Universidade Católica do Paraná, PUCPR, Curitiba, Brazil [email protected]

Luiz Lebensztajn

Viviana Cocco Mariani

Laboratório de Eletromagnetismo Aplicado, LMAG-PEA, Universidade de São Paulo, USP, São Paulo, Brazil [email protected]

Pós-Graduação em Engenharia Mecânica, PPGEM, Pontifícia Universidade Católica do Paraná Departamento de Engenharia Elétrica, Universidade Federal do Paraná, UFPR, Curitiba, Brazil [email protected]

Leandro dos Santos Coelho Pós-Graduação em Engenharia de Produção e Sistemas, PPGEPS, Pontifícia Universidade Católica do Paraná Departamento de Engenharia Elétrica, Universidade Federal do Paraná, UFPR, Curitiba, Brazil [email protected] Abstract— Metaheuristics of the natural computing field have been proposed as an alternative to mathematical optimization approaches to address non convex problems involving large search spaces. In recent years a new optimization metaheuristic algorithm was proposed called Wind Driven Optimization (WDO). WDO is a stochastic nature-inspired paradigm based on atmospheric motion. In this paper, a modified version of WDO is proposed and evaluated, based on Lévy flights (or Lévy motions) to tune its control parameters, called Lévy WDO (LWDO). Lévy flight or anomalous diffusion process is a random walk characterized by Markov chain in which the step-lengths have a probability distribution that is heavy-tailed. To evaluate the multiobjective optimization performance of the WDO and the proposed LWDO, a benchmark for optimizing of a safety isolating transformer is adopted. In this paper, the transformer design optimization is treated as a multiobjective problem, with the aim to minimize both the total mass (iron and copper materials) and losses taking into consideration design constraints. Simulation results testify that the multiobjective LWDO is a promising approach for multiobjective optimization as it outperforms the WDO in multiobjective version and the classical NSGA-II (Non-dominated Sorting Genetic Algorithm, version II). Keywords—optimization; transformer design; metaheuristics; wind driven optimization

I. INTRODUCTION Various fields of science and engineering often have multiple conflicting objectives. In a single-objective optimization problem, the optimal solution is clearly defined whereas in multiobjective optimization problem there exists a

c 978-1-5090-0623-6/16/$31.00 2016 IEEE

set of trade-offs giving rise to numerous solutions involving the simultaneous optimization of competing objectives. In the absence of any preference information, a non-dominated set of solutions is obtained by multiobjective optimization, instead of a single optimal solution. These solutions, called Pareto optimal solutions, are optimal in the wider sense that no other solutions in the search space are superior to them when all objectives are considered. There are two distinct goals, common but often conflicting, in multiobjective optimization: convergence and diversity. Achieving a balance between convergence and diversity is relevant, but usually far from trivial, in multiobjective optimization. Classical methods do convert the multiobjective optimization problem (MOP) to a single objective optimization problem by a suitable scaling/weighting factor method. On the other hand, optimization metaheuristic algorithms have been proposed as an alternative to classical optimization methods in various applications. They are very popular, well-known global optimization schemes and attempt to reproduce social behavior or natural phenomena. Furthermore, many metaheuristic algorithms have been proposed for dealing with real-world and large scale problems which are difficult to solve by traditional methods. Evolutionary algorithms and swarm intelligence approaches have been recognized to be well-suited metaheuristics to solve MOPs [1]-[3], because they deal simultaneously with a set of possible solutions (Pareto front) and are less susceptible to the shape or continuity of Pareto front.

4642

In this paper, a multiobjective version and an improved approach of the wind driven optimization (WDO) [4] is proposed using Lévy flights (LWDO). The original design of the WDO was inspired by earth’s atmosphere where the wind flow participates of an attempt to balance horizontally the air pressure and it can be effective in solving multidimensional numerical optimization problems. However, the WDO as other population-based algorithms can present premature convergence and stagnate at suboptimal solutions. The control parameters are generally the key factors affecting the WDO’s convergence. Nevertheless, the best parameter setting can be different for different problems. As a consequence, in order to have a proper parameter setting, the tedious and highly time-consuming trial-and-error method is usually required. The utilization of different distributions in tuning of WDO’s control parameters can be useful to escape more easily from local minima and to avoid the premature convergence than with the classical WDO approach. In this context, Lévy flights can be useful in the searching process. Lévy flights are random walks determined by step lengths that are drawn from power low distributions with heavy power law tails. The focus of application of the proposed LWDO is the multiobjective transformer design optimization (TDO). Transformers are regarded as crucial components in electrical power systems. The industry and energy market mandates the adoption of transformers design strategies yielding better performance at lower costs. Several optimization techniques have addressed TDO problems using single-objective approaches (see details in [5],[6]) including optimization metaheuristics such as genetic algorithms [7], particle swarm optimization [8], covariance matrix adaptation evolution strategy [9], and bacterial foraging algorithm [10]. Through the analysis of the transformers’ characteristics and the relationship between design parameters, a TDO can be established as a constrained multiobjective optimization problem (MOP). Examples of TDO can be viewed in [11]. The purpose of this paper is to present a TDO methodology using a multiobjective WDO approach. A benchmark for optimizing a safety isolating transformer is adopted to validate the optimization using the WDO and the proposed LWDO in multiobjetive versions. In this paper, the TDO is approached as a multiobjective problem, with the aim to minimize both the total mass (f1) and the losses (f2), taking into account the design constraints. Simulation results demonstrate the effectiveness and competitiveness of the proposed LWDO algorithm. The remainder of the paper is organized as follows. Section II gives some fundamentals of the TDO. In Section III, the fundamentals of the WDO and LWDO are detailed. In Section IV, the results are presented and discussed. The paper is concluded in Section V.

II. TRANSFORMER DESIGN OPTIMIZATION PROBLEM Transformers are widely used in electric power system to perform the primary functions, such as voltage transformation and isolation. In this new and challenging environment, there is an urgent need for the transformer manufacturing industry to improve transformer efficiency and to reduce costs, since high quality, low cost products have become the key to survival [10]. The aim of TDO is to obtain the dimensions of parts of the transformer in order to supply these data to the manufacturer. The transformer design should be carried out based on the specification given, using available materials economically in order to achieve low cost, low weight, small size and good operating performance [12]. In general, the problem of TDO can be based on minimization or maximization of one or more objective functions subjected to several constraints. As mentioned in [13], the TDO problem adopted in this paper contains seven discrete design variables. For this TDO problem, the core of the transformer has a particular shape, entitled EI, which can be observed at [13]. This kind of shape is commonly used for monophasic transformers. There are three parameters {a, b, c} for the shape of the lamination, one for the frame {d}, two for the section of enamelled wires {S1,S2}, and one for the number of primary turn {n1}. There are 24 types of lamination, 62 combinations between the laminations EI and the frames, and 62 types of enamelled wires. The number of primary turn n1 is integer but only 1,000 values are allowed, leading to 246,078,000 possible combinations. There are six inequality constraints applied on this problem. The copper and iron temperatures should be less than 120°C and 100°C, respectively. The magnetizing current Iȝ/I1 and drop voltage DV2/V2 are less than 10%. Finally, the filling factor of both coils is lower than 0.5. The objective functions are to minimize the total mass Mtot of iron and copper materials and the losses, i.e., maximize the efficiency (Ș). III. DESCRIPTION OF THE WDO AND LWDO In the following sub-sections, the procedures of the WDO and the LWDO are described. First, the formulation of the classical WDO design is provided, and then the proposed LWDO is detailed. A. The Classical WDO The WDO is a recently proposed nature inspired metaheuristic technique. It is a population based iterative heuristic process inspired on atmospheric motion. The WDO was proposed by Bayraktar et al. in 2010 [4]. The WDO algorithm uses the wind flow from earth’s atmosphere as inspiration. The wind flow is generated due to the radiation of the sun over the surface and the heating emerged from that phenomenon, which occurs irregularly, generating fluctuations of temperature that generates differences of density on the air [14]. This dynamic movement comes from the pressure gradient. The representation of the atmospheric motion in the WDO uses the Lagrangian

2016 IEEE Congress on Evolutionary Computation (CEC)

4643

mathematical representation due to its infinitesimal fluid parcels collection that can be governed by Newton’s second law of motion. In the case of wind flow, it is assumed that the atmosphere is a homogeneous fluid and that hydrostatic balance is present. Considering the Cartesian coordinate description (x, y, z) and that the horizontal movement is stronger than the vertical movement the flow can be treated as horizontal motion only [15]. The movement comes from the pressure gradient ∇P that expressed in rectangular coordinates system have the form given by § ∂P ∂P ∂P · ∇P = ¨¨ , , ¸¸ © ∂x ∂y ∂z ¹

(1)

where P is the pressure. Considering the fact that the air has finite mass and finite volume ( δV ), the wind blows in the direction of high pressure area to a low pressure area with a velocity (V) is proportional to the pressure gradient force (FPG), given by FPG = −∇P ⋅ δV

(2)

The representation of the atmospheric motion in this model uses the Lagrangian description [15] due to its infinitesimal fluid parcels collection that can be over-ruled by Newton’s second law of motion. This representation of the WDO considers each air parcel dimensionless, what do not complicate the implementation by separate coordinates for each corner of the cuboid. In addition, each face of the cuboid can receive different pressure values causing its deformation by several forms. In the case of wind flow, it is assumed that the atmosphere is a homogeneous fluid and that hydrostatic balance is present. Considering the fact of the rectangular coordinate decryption and that the horizontal movement is stronger than the vertical movement, the wind blows can be treated as horizontal motion only [15]. However, the WDO algorithm can operate under Ndimensional search space considering certain assumptions and simplifications. The process of the air parcel trajectory calculation starts from the Newton’s second law of motion, which states that the total force applied on an air parcel causes its acceleration according to

ρa = ¦ Fi

(3)

where ȡ is the air density, a is the acceleration and Fi is the external force. The relation of the air parcel pressure P and its density is given by P = ρRT

(4)

where R is the ideal gas constant and T is the temperature. In this model, the pressure gradient can be considered the fundamental force that initiates the motion of the air parcel,

4644

although, that are four different forces that can cause influence on the movement such as pressure gradient force defined previously, friction force (α), gravitational constant (g) and Coriolis effect (c), given respectively by equations given by FF = −ραv

(5)

FG = ρ ⋅ δV ⋅ g

(6)

FC = −2Ωxv

(7)

where Ω represents the rotation of the earth and v is the velocity. These four forces are the major contributors to the motion of air parcel, but there are other forces not considered, for example, advection and turbulent drag force. Assembling all the forces described in the equation and using the ideal gas law the density can be written in terms of the pressure, where the index c corresponds to the current value of the parameter. In this case, we have the velocity (v) and the position of the air parcel (x) updated by § − 2ΩxvRT · § RT · ¸ ¸ − αu + ¨ Δv = g + ¨¨ − ∇P ¨ ¸ ¸ Pc ¹ Pc © ¹ ©

(8)

Considering that Δv = vnew − vc , that the acceleration of

(

)

gravity g =| g | (0 − xc ) and ∇Pres = Presopt − Presc xopt − xc ,

where the index opt is the optimum, we rewrite the equation (8) assuming Pc as index i to obtain the equations (9) and (10) where Δt = 1 . In this case, we have

(

)

· §1 v new = (1 − α )v c − gx c + ¨¨ − 1 x opt − x c RT ¸¸ ¹ ©i § − 2ΩxvRT · +¨ ¸ i © ¹ xnew = xc + vnew Δt

(9)

(10)

The implementation of the WDO starts with the initialization of the parameters population size, maximum number of iterations, coefficients, boundaries and pressure function definition that compounds the objective function. Than allocate random positions and velocities for the air parcels for the evaluation of the pressure for each air parcel and the update of velocity and position and it is checking under the velocity limits and boundaries conditions until the maximum number of iterations as stopping criterion. B. The Proposed LWDO In the classical WDO, there are multiple control parameters that must be chosen prior to starting an optimization using WDO, including the R, T, g, α, and c. The stability of the WDO performance against parameter settings is generally unknown. In contrast, to improve the diversity of the multiobjective WDO the proposed LWDO approach uses Lévy flights [6] in order to tune all control parameters.


Lα ,γ =

1 ∞ α· § ³ exp¨ − γq ¸ cos(qz )dq © ¹ π 0

(11)

The parameters characterizing the distribution are α, controlling the sharpness of the graph, and γ, controlling the scale unit of the distribution. For α = 2, the distribution is equivalent to the Gaussian distribution. For α = 1, the integral can in principle be carried out analytically and is equivalent to the Cauchy probability distribution. Thus, the smaller α parameter causes the distribution to make longer jumps since there will be longer tail. The sign of the skewness parameter γ indicates the skew direction, positive to the right, and negative to the left. When γ = 0, the distribution is symmetric [18]. In the focus of this paper, a truncated Lévy flight in the range [0,λ] to tune the control parameters R, T, g, α, and c with λ equal to 3, 3, 0.2, 0.4 and 0.3, respectively, during the iterations cycle. This approach is adopted in the multiobjective LWDO in order to maintain the population diversity and consequently a better distribution of Pareto-optimal solutions. C. WDO and LWDO in Multiobjective Version To deal with conflicting objectives, one should adapt the WDO and LWDO algorithms such that the goals in multiobjective optimization are attained. Namely, the solutions which approximate the true Pareto front are expected to be as close to its values and as diverse as possible, in such a way that at the end of the optimization procedure the designer has representative solutions to choose from. To this end, the present work adopts the same truncation procedure based in non dominance and crowding distance as in [20]. Thus, at the end of each iteration, the set containing the solutions formed by the current and previous iterations are sorted according to (i) non dominance and (ii) crowding distance and then truncated such that the number of solutions for the next iteration is kept constant. The best individual at each iteration is chosen randomly in the 10% first sorted individuals according to the two aforementioned criteria, as in [21]. There are inequality constraints to the approached TDO problem. So, the WDO and LWDO multiobjective versions given, respectively, by acronyms MO-WDO and MO-LWDO must be supplemented with a mechanism to efficiently handle

constraints. In this paper, a third objective function f3 to be minimized related to constraints-handling is adopted. IV. OPTIMIZATION RESULTS Both MO-WDO and MO-LWDO algorithms were run 30 times with a population size of 100 and the stopping criterion equal to 50,000 function evaluations. The control parameters adopted for MO-WDO were set according to [4], where -2ȍRT is equal to 0.4, RT is equal to 3, g is equal to 0.2, α is equal to 0.4 and maximum velocity is equal to 0.3. In multiobjective optimization, the goal of the MO-WDO or MO-LWDO is to search a solution set of which the corresponding objective vectors are closely and evenly distributed on the Pareto front. In terms of optimization results evaluation, the NSGA-II (Non-dominated Sorting Genetic Algorithm, version II) [20] was adopted to perform a comparison among the WDO based algorithms. The parameters for NSGA-II were a crossover probability equal to 0.8, a mutation probability equal to 0.25. The population size and the stopping criterion are the same of the MO-WDO and MO-LWDO. The total Pareto solutions (all feasible solutions, where f3 0) obtained in 30 runs using the data of the Pareto front filtered out of 30 runs for MO-WDO, MO-LWDO and NSGA-II were 44, 300 and 89, respectively, as shown in Table I. On the other hand, in terms of Euclidean distance of the objective function values with the origin at (0, 0), the MOLWDO presented the best (minimum) value of the tested optimization approaches. TABLE I.

PERFORMANCE INDICES FOR THE TDO (30 RUNS) Indices

MOMOWDO LWDO Pareto front of feasible solutions * 44 300 Normalized spacing (f1, f2) * 0.1988 0.0038 Normalized Euclidean distance (f1, f2) * 0.7301 0.3112 *this index uses the data of the Pareto front filtered out of 30 runs

NSGA II 89 0.0034 0.4398

The resulting algorithm exhibits promising performance in terms of the Pareto fronts as shown in Fig. 1.

MO-WDO MO-LWDO NSGA-II

0.12 0.1 f2: 1-η

Non-Gaussian, heavy-tailed statistics is becoming a commonly used tool in several applications, such as behavior of numerous animals and insects [16],[17]. Lévy distribution or anomalous diffusion is stable and has probability density function that is analytically expressible. Lévy flights processes describe a class of random walks whose step lengths follow a power-law tailed distribution. Furthermore, the Lévy flight random walk patterns have infinite variance (except for the Gaussian case) and possess scale-invariance and selfsimilarity properties. This flight behavior has been applied to optimization and search algorithms, and the reported results show its importance in the field of solution search algorithms [18],[19]. The Lévy distribution has the probability density:

0.08 0.06 0.04 2

4

6 8 f1: M total (kg)

10

12

Fig. 1. Pareto set points (filtered of 30 runs) using MO-WDO, MO-LWDO and NSGA-II.


4645

According to the results, the MO-LWDO presented promising results in terms of the quality of the solutions found in the Pareto front. Furthermore, MO-LWDO solution sets have better Pareto front convergence than MO-WDO and NSGA-II preserving solution set diversity. The distribution of the non-dominated solutions along the variable space is another important point of analysis. Figs. 2-4 show the histograms of the optimization variables for the nondominated solutions obtained by the MO-WDO, MOLWDO and NSGA-II, respectively. The proposed method (MO-LWDO) showed the best distribution in the Pareto front. By observing the histograms of optimization variables, one can notice that MO-LWDO has a more continuous distribution in the following variables: number of turns of the primary, wire section in primary and wire section in the secondary. The shapes and frame (a, b, c and d) have similar distribution for all the three methods. Without a good distribution on number of turns of the primary, wire section in primary and in the secondary a Pareto set is not obtained in this problem.

Fig. 3. Histogram of the optimization variables for the MO-LWDO (filtered of 30 runs).

Fig. 2. Histogram of the optimization variables for the MO-WDO (filtered of 30 runs).

Fig. 4. Histogram of the optimization variables for the NSGA-II (filtered of 30 runs).

4646


V. CONCLUDING REMARKS For MOPs, evolutionary and swarm intelligence algorithms in general have demonstrated to be effective and efficient tools for finding approximations of the Pareto front. In this context, the present paper evaluated the use of WDO in multiobjective optimization applied to electromagnetics problems. Moreover, we introduced a novel approach based on Lévy flights which improved the final results in comparison to NSGA-II. This paper proposes multiobjective MO-WDO and MOLWDO algorithms for TDO, which uses an external archiving scheme and the ranking with non-dominance and crowding distance. The mentioned optimization algorithms have been successfully applied to the TDO problem, mainly the MOLWDO with superior results when compared with NSGA-II. Further research work is being carried out to efficiently optimize other constrained TDO cases. How to improve the efficiency of the MO-LWDO remains to be studied further. Many possibilities are present in the literature with respect to the improvement of global search algorithms [22]-[26]. Other research line will focus on approaches which adapt the technique to the problem at hand. ACKNOWLEDGMENTS This work was supported by the National Council of Scientific and Technologic Development of Brazil — CNPq — under Grants 303908/2015-7/PQ and 303906/2015-4/PQ. REFERENCES [1]

[2] [3]

[4]

[5]

[6]

[7]

[8]

[9]

A. Zhou, B. -Y. Qu, H. Li, S. -Z. Zhao, P. N. Suganthan, and Q. Zhang, “Multiobjective evolutionary algorithms: a survey of the state of the art,” Swarm and Evolutionary Computation, vol. 1, no. 1, pp. 32-49, 2011. K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms, John Wiley & Sons, Baffins Lane, Chichester, 2001. C. A. Coello Coello, “Evolutionary multiobjective optimization: A historical view of the field,” IEEE Computational Intelligence Magazine, vol. 1, no. 1, pp. 28-36, 2006. Z. Bayraktar, J. P. Turpin, and D. H. Werner, “Nature-inspired optimization of high-impedance metasurfaces with ultrasmall interwoven unit cells,” IEEE Antennas and Wireless Propagation Letters, vol. 10, pp. 1563-1566, 2011. A. Khatri and O. P. Rahi, “Optimal design of transformer: a compressive bibliographical survey,” International Journal of Scientific Engineering and Technology, vol. 1, no. 2, pp. 159-167, 2012. H. D. Mehta and R. M. Patel, “A review on transformer design optimization and performance analysis using artificial intelligence techniques,” International Journal of Science and Research, vol. 3, no. 9, pp. 726-733, 2014. A. K. Yadav, O. P. Rahi, H. Malik, and A. Azeem, “Design optimization of high-frequency power transformer by genetic algorithm and simulated annealing,” International Journal of Electrical and Computer Engineering, vol. 1, no. 2, pp. 102-109, 2011. V. Rashtchi, V. Bagheri, A. Shabani, and S. Fazli, “A novel PSO-based technique for optimal design of protective current transformers,” International Journal for Computation and Mathematics in Electrical and Electronic Engineering (COMPEL), vol. 30, pp. 505-518, 2011. S. Tamilselvi and S. Baskar, “Modified parameter optimization of distribution transformer design using covariance matrix adaptation evolution strategy,” International Journal of Electrical Power & Energy Systems, vol. 61, pp. 208-218, 2014.

[10] S. Subramanian and S. Padma, “Optimization of transformer design using bacterial foraging algorithm,” International Journal of Computer Applications, vol. 19, no. 3, pp. 52-57, 2011. [11] E. I. Amoiralis, P. S. Georgilakis, M. A. Tsili, and A. G. Kladas, “Global transformer optimization method using evolutionary design and numerical field computation,” IEEE Transactions on Magnetics, vol. 45, no. 3, pp. 1720-1723, 2009. [12] A. K. Yadav, A. Azeem, H. Malik, and O. P. Rahi, “Optimization of power transformer design using simulated annealing technique,” International Journal of Electrical Engineering, vol. 4, no. 2, pp. 191198, 2011. [13] http://l2ep.univ-lille1.fr/come/benchmarktransformer_fichiers/Optimization_Problem.htm [14] C. D. Ahrens, Meteorology Today: An Introduction to Weather, Climate, and the Environment. 7th ed. Belmont, CA, USA: ThomsonBrook, 2003. [15] R. D. Thompson, Atmospheric Processes and System. New York, NY, USA: Routledge, 1998. [16] M. Ghaemi, Z. Zabihinpour, and Y. Asgari, “Computer simulation study of the Lévy flight process,” Physica A, vol. 388, no. 8, pp. 15091514, 2009. [17] G. Viswanathan, S. Buldyrev, S. Havlin, M. da Luz, E. Raposo, and H. Stanley, “Optimizing the success of random searches,” Nature, vol. 401, pp. 911-914, 1999. [18] H. Hakl and H. U÷uz, “A novel particle swarm optimization algorithm with Lévy flight,” Applied Soft Computing, vol. 23, pp. 333-345, 2014. [19] L. S. Coelho, T. C. Bora, and V. C. Mariani, “Differential evolution based on truncated Lévy-type flights and population diversity measure to solve economic load dispatch problems,” International Journal of Electrical Power & Energy Systems, vol. 57, pp. 178-188, 2014. [20] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multi-objective genetic algorithm: NSGA-II,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182-197, 2002. [21] C. R. Raquel and P.C. Naval, Jr., “An effective use of crowding distance in multiobjective particle swarm optimization,” in Proceedings of GECCO, Washington, DC, USA, 2005. [22] S. Mahdavi, M. E. Shiri and S. Rahnamayan, “Metaheuristics in largescale global continues optimization: A survey,” Information Sciences, vol. 295, pp 407-428, 2015. [23] G. Sun, R. Zhao, and Y. Lan, “Joint operations algorithm for large-scale global optimization,” Applied Soft Computing, vol. 38, pp. 1025-1039, 2016. [24] E. Sayed, D. Essam, R. Sarker, and S. Elsayed, “Decomposition-based evolutionary algorithm for large scale constrained problems,” Information Sciences, vol. 316, pp. 457-486, 2015. [25] Q. Fan and Y. Zhang, “Self-adaptive differential evolution algorithm with crossover strategies adaptation and its application in parameter estimation,” Chemometrics and Intelligent Laboratory Systems, vol. 151, pp. 164-171, 2016. [26] M. Mirjalili, S. Saremi, S. M. Mirjalili, and L. S. Coelho, “Multiobjective grey wolf optimizer: a novel algorithm for multi-criterion optimization,” Expert Systems with Applications, vol. 47, pp. 106-119, 2016.


4647