Morphing Wing Structural Optimization Using Opposite-Based Population-Based Incremental Learning and Multigrid Ground Elements

This paper has twin aims. Firstly, a multigrid design approach for optimization of an unconventional morphing wing is proposed. The structural design problem is assigned to optimize wing mass, lift effectiveness, and buckling factor subject to structural safety requirements. Design variables consist of partial topology, nodal positions, and component sizes of a wing internal structure. Such a design process can be accomplished by using multiple resolutions of ground elements, which is called a multigrid approach. Secondly, an opposite-based multiobjective population-based incremental learning (OMPBIL) is proposed for comparison with the original multiobjective population-based incremental learning (MPBIL). Multiobjective design problems with single-grid and multigrid design variables are then posed and tackled byOMPBIL andMPBIL.The results show that usingOMPBIL in combination with a multigrid design approach is the best design strategy. OMPBIL is superior to MPBIL since the former provides better population diversity. Aeroelastic trim for an elastic morphing wing is also presented.


Introduction
Design/optimization of a morphing wing aircraft structure is a research field investigated throughout the world.It has been found that the shape and structural flexibility of a morphing aircraft greatly affect its flight performance [1].This leads to continuous shape change and is called a morphing concept, which is an attempt to avoid using conventional hinged control surfaces.From the literature, external plane transformation of the wing is a recent popular morphing concept [1].The idea is to change wing aerodynamic performance by deforming the initial wing plane, which can be categorized as wing camber variation [1][2][3][4][5], lateral wing bending [1], and wing twisting [5][6][7][8][9].The variable camber of an aircraft wing is the most popular method [1], which can be implemented via a conventional hinged mechanism, a smart material, or a compliant mechanism.The compliant mechanism is the best method for such a morphing aircraft wing concept [4].In the lateral wing bending concept, the wing can bend up and down in a vertical direction so as to produce a wing span camber.Using this concept, wing aerodynamic performance can be varied and it is carried out by using an internal mechanism.The last concept, the wing twisting, can adjust a wing angle of attack by using structural flexibility and applying external torque.Later, there has been some work attempting to exploit adaptive internal structures [6].The concept can be thought of as a combined lateral wing bending and variable camber concept.From the literature [1][2][3][4], it is revealed that the morphing wing aircraft structure is actuated directly by an external force to carry out such aircraft control.Furthermore, the morphing wing structures are elastic rather than perfectly rigid; thus, consideration of their aeroelastic behavior is essential in a design process.This leads to recent work in which the effects of the external force on aeroelastic characteristic of a morphing aircraft wing are studied [10].The result shows that the actuating force has a significant impact on the aeroelastic and mechanical characteristics and these effects should be taken into account during the design process of a morphing aircraft structure.Also, the aircraft internal structure including the structural layout and sizes significantly affects the aeroelastic characteristics.
In our recent work, it has been proposed to synthesize new internal structural layouts of an aircraft wing [8] and an aircraft morphing wing [9] using multiobjective optimization.The first study leads to an unconventional wing internal layout, which is aeroelastically superior to a traditional ribspar layout.Design variables include partial topology and sizes of segments on a structure, while design objectives are wing weight, lift effectiveness, and buckling factor.The second study found the unconventional wing structures subject to external actuating torques which are practicable to apply for the morphing aircraft.Design variables include partial topology and sizes of segments on a structure like the first study, but the second study included nodal displacements.Design objectives are wing weight, percentage of change in lift effectiveness, and buckling factor.In both study, the combination of the two types of design variables is carried out by using a ground element approach.The method used to tackle the design problems is multiobjective populationbased incremental learning (MPBIL).The optimizer used herein is one of the most powerful evolutionary methods for solving structural topology optimization [11].The method is category as nongradient optimization methods, which does not require gradient information to converge to a solution.The design strategy proposed in these works can be used as a numerical tool for synthesizing an internal structural layout of a morphing wing, which usually requires high structural and aeroelastic performance for flight operation [8,9].
The main motivation of this research work is to initiate an idea to synthesize unconventional aircraft wing structures that are expected to outperform their conventional counterparts.This paper is intended as an extension of the above literature.The work has two aims.Firstly, the use of several resolutions of ground elements at the same time is proposed instead of using one ground element resolution for a design problem and this is termed a multigrid ground element approach.Design variables are simultaneous topology, shape, and sizing optimization.The second part is performance enhancement of a multiobjective evolutionary optimizer MPBIL.The opposite-based evolutionary optimization is used to improve the search performance of the original MPBIL and it is named opposite-based multiobjective PBIL (abbreviated as OMPBIL).MPBIL and OMPBIL are employed to solve a morphing wing synthesizing design problem with the use of single-grid and multigrid ground elements.Optimum results reveal that OMPBIL is superior to MPBIL.The multigrid approach leads to better design results than the single-grid one, while the set of design variables that include structural layout, nodal positions, and sizes gives the best structural layout.
The rest of this paper is organized as follows.Section 2 details single-and multiground element design approaches for synthesizing a wing's internal structure.The oppositebased multiobjective population-based incremental learning is proposed in Section 3. A design problem and its conditions as well as a numerical experiment are given in Section 4, while the design results are in Section 5.The conclusions and discussion of the study are finally drawn in Section 6.

Partial Topology, Shape, and Sizing Optimization
Shape and sizing optimization [12] and simultaneous partial topology integrate three different types of structural design variables in the same design problem.This is somewhat complicated but found to be an efficient design strategy [8,9,13].Such a design strategy can be achieved by using the ground finite element technique.The ground element technique is commonly used in a topology optimization process.The basic idea is to generate ground finite elements throughout a design domain being considered.Design variables determine pseudo-densities (such as element thickness and modulus) of those ground elements.Having obtained an optimum solution (e.g., from compliance minimization or dynamic stiffness maximization), ground elements with low pseudo-densities are assigned as holes on a structure whereas elements with high pseudo-density represent a structure.This concept has been proven effective for many design cases [11,14,15].For simultaneous topology, shape, and sizing optimization of an aircraft wing, ground elements include all possible combinations of wing internal segments and wing skins.Figures 1(a) and 1(b) display the ground elements or ground segments for the wing internal structure used for design demonstration in this work.The details of this simple wingbox structure can be found in our previous work [8,9].Note that diagonal segments of a structure are used because it has been found that such ground segments result in aeroelastically superior structures compared to using only chordwise and spanwise segments as with a conventional wing structure [8].The upper and lower wing skins shown in Figure 2 are also set as design variables.It has been found that the variable-thickness of element panels can enhance its aeroelastic performance [5,16].The combination of topology and sizing variables of the wing is carried out by assigning values to the thicknesses of those segments in meters (m).For example, let   ∈ [0.0005, 0.001] m be design variables for segment thickness where the lower and upper bounds are based on manufacturing tolerance.The th ground segment will be removed from the main structure if the value of   is close to its lower bound at the optimum point.Other thickness values higher than the lower bounds are used for structural sizing.It should be noted that the ideal lower bounds for topological design variables are zero segment thickness.However, in order to prevent singularity in a global stiffness matrix of a wing, we have to use a small positive value as lower bounds of the topological design variables.Figure 1 illustrates how to transform from particular design variables (Figure 1(a)) to be a wing structure (Figure 1(b)).Shape design variables, on the other hand, can be added to the partial topology and sizing optimization by assigning positions of internal nodes of those ground segments as design variables as illustrated in Figure 3.The shape design variables determine the positions of the external nodes in and -directions.Figure 4 shows the segments on lower and upper wing skins as the internal nodes positions are varied.The skin segments are also set as design variables.
The design problem is posed to minimize three design objectives simultaneously, while fulfilling safety constraints.A morphing wing in this study is based on the wing twisting concept.To achieve flight control, the wing is subject to a twisting actuating force on the wing internal structure as shown in Figure 5.This means that the wing generates lift variation through change of lift effectiveness due to actuation.Thus, the first objective function is the difference percentage between wing lift effectiveness with and wing lift effectiveness without actuating force.This function can be expressed as where  ,Fe is the lift effectiveness considering the effect of an external force.The second objective function is set to minimize structural mass, which is rather an unavoidable requirement for aircraft structural design.The third objective, buckling factor maximization, is set to increase the strength of the structure.The wing is considered as a thin-wall structure in which one of the most critical design criteria is buckling instability.The multiobjective design problem can be expressed as where  is the structural mass,   is the divergence speed,   is the flutter speed of the wing,  is the safety factor,  is the maximum equivalent stress on the wing due to aerodynamic and actuating forces,   is the yield stress,  is the buckling factor,   is the thickness of a wing segment ,   is the nodal position in the chordwise direction of node , and   is the nodal position in the spanwise direction of node .The last two constraints are assigned so as to prevent node positions and wing segments overlapping each other in cases of performing simultaneous partial topology, shape, and sizing design.For details of structural and aeroelastic models as well as other function evaluation procedures, see Sleesongsom and Bureerat [8,9].A question always arises at the preprocess stage [17], when using a ground element approach: what is the best ground element resolution for a design problem?As a result, using several sets of ground segments at the same time when performing optimization is investigated and it is termed a multigrid design approach.The multigrid approach is said to be an extension of the ground segment strategy, which has been proposed to solve a truss structure [18].Figure 6 shows three sets of ground segments with different grid resolutions used in this study.The first set of ground elements has 99 segments, which was used in the previous study [8].The second set has 82 segments, while the third set has 53 segments.Using one set of ground elements is simple to perform since we can assign thickness values to all segments (and to nodal positions in cases of shape design variables).Nevertheless, when using many ground segment sets, a special encoding/decoding scheme for multiobjective evolutionary optimization is needed [8].From Figure 6 , the variables encoding/decoding scheme for the multigrid design approach can be detailed as in Algorithm 1.As MPBIL and OMPBIL use binary design variables, the conversion of a binary string to become a real design vector x needs to be performed before entering Algorithm 1.

MPBIL and OMPBIL
This section details the two algorithms MPBIL and the new optimizer OMPBIL, which will be used to tackle the optimization presented in Section 2. Evolutionary algorithms (EAs) are alternative optimizers to classical gradient-based optimizers such as a steepest descent method and sequential linear programming.Unlike the gradient-based optimizers, EAs search for optima based upon evolutionary concepts in combination with randomization [11].The methods, for example, genetic algorithms, are simple to use, robust, and capable of tackling global optimization.They can be used for almost any kinds of optimum design problems (such as discrete [11], continuous [14], and combinatorial optimization [19]) since they can search for optimum solutions without using function derivatives.In our design problem, there is no accuracy guarantee for function evaluation especially the flutter speed [10].Apart from that, our new design variables are mixed integer/continuous where it is difficult to apply any gradient-based optimizer to solve the design problem.The successful use of EAs for complex engineering design has been reported, for example, in Lencus et al. [20] and G. Nowak and I. Nowak [21].EAs for multiobjective optimization, traditionally termed multiobjective evolutionary algorithms (MOEAs), have one outstanding feature; that is, they can explore a nondominated front within one simulation run.Although having several advantages, two unavoidable drawbacks of EAs and MOEAs are the lack of convergence rate and consistency [18].This makes EAs ineffective when solving topology optimization, which is said to be a large-scale design problem with a large number of design variables.Over the past two decades, many attempts to improve search performance of EAs for solving structural topology optimization have been made, for example, a morphological geometric representation approach [22], a chromosome repairing technique [23], a ground element filtering technique [11], and antioptimization [24].The work in Kunakote and Bureerat [11] presents the comparative performance of MOEAs for solving structural topology optimization test problems based on the ground element filtering technique.The established MOEAs included MPBIL, Pareto archive evolution strategy (PAES) [25], nondominated sorting genetic algorithm (NSGA) [26], strength Pareto evolutionary algorithm (SPEA) [27], and multiobjective particle swarm optimization (MPSO) [28].The results show that MPBIL is the unanimous best optimizer [29].The simple concept of estimation of distribution algorithm embedded in MPBIL gives an advantage when dealing with multiobjective topology optimization.As a result, MPBIL is the only MOEA  chosen to solve this design problem.This section briefly details the concept of MPBIL and then it is extended to the opposition-based multiobjective population-based incremental learning.

Multiobjective Population-Based Incremental Learning.
The principle of a PBIL search is governed by the learning rule, which is iteratively limiting the design space depending on current best design variables, and a random process [29].The design space is iteratively narrowed until the optimum is approached.Unlike the popular genetic algorithm and other evolutionary algorithms (EAs), a binary population used in PBIL is represented by an estimated probability distribution for each binary bit, which is called a probability vector.Table 1 shows three examples of probability vectors and their corresponding binary populations given that a design solution is a row of a population.The th element of a probability vector determines the probability of having a string "1" on the th column of the population.It can be seen that one probability vector can produce a variety of binary populations.Unlike single objective PBIL, MPBIL uses several probability vectors to represent a binary population in order  to maintain population diversity.Thus, a set of probability vectors is called a probability matrix [29].The search procedure of MPBIL starts with an initial Pareto archive for collecting nondominated solutions and a probability matrix whose elements are full of "0.5." Each row of the probability matrix is a probability vector that will be used to create a subpopulation.Let   be the number of design solutions in a population,  the number of probability vectors, and   the number of binary bits for one design vector x.The probability matrix, therefore, has the size of  ×   where each row of the matrix results in approximately   / design solutions as one subpopulation.Each design variable is represented by a binary string with   / elements.Figure 7 displays a 3 × 4 probability matrix and its corresponding population (  = 12), which is the union set of all three subpopulations according to the three rows of the probability matrix [30].
Having generated the population and evaluated their corresponding objective values, the nondominated members sorted from the combination of members in the current population and the old external Pareto archive are taken as a new external Pareto set.If the external Pareto set is full, some solutions are removed from the external Pareto set using an archiving technique [29].In MOEAs, the archiving technique is proposed to filter some nondominated solutions out of the archive matrix to prevent excessive memory usage of a computer because MOEAs can normally explore a large number of Pareto solutions within one run.The probability matrix and the nondominated solutions are improved iteratively until the termination criterion is met.
For the updating scheme for the probability matrix on each generation, the matrix is updated in such a way that  0 <   binary solutions are randomly selected from the current Pareto archive to modify the th row of the probability matrix.The th row of the probability matrix is updated using the relation where   ∈ (0, 1) is called the learning rate to be defined and   is the average value of the th bit position of those  0 randomly selected nondominated solutions.In our work, the learning rate is set as where rand ∈ [0, 1] is a uniform random number.It is also useful to apply mutation to the probability matrix at some predefined probability as where ms (default = 0.2) is the amount of shift used in the mutation.

Opposite-Based Multiobjective Population-Based
Incremental Learning.MPBIL is governed by the probability update rule (3), which is a variation of Hebb's rule one of the first neural network learning laws [31].From (3), if   approaches one, the probability update rule quickly forgets  old  and remembers only the most recent   .Otherwise, the probability update rule remembers only  old  and quickly forgets   .A small value of learning rate is usually recommended for the conventional PBIL [32].This   setting is not always effective for MPBIL as showed in Bureerat and Srisomporn [30], which had good results from randomly defining   value in the range of [0.4,0.6].The question arises: how can we choose a proper value of   for a general problem, which can accelerate a convergence rate or increase search performance?A new idea to deal with such a problem is motivated from Rahnamayan et al. [33] and Wang et al. [34] as they proposed to enhance EAs using oppositionbased learning (OBL).Such a concept could be used for improving population initialization [33] or a hybrid with EA main operators [34,35].Both methods can enhance EA search performance.Furthermore, it has been found that this idea was used to maintain diversity in PBIL and increase its performance [36] for single objective design cases.
Opposition is concerned with the relationship between entities, objects, or their abstractions of the same nature, which are completely different in some manner [37].For example, "remember" and "forget" are descriptions of the brain, but completely different.The opposition-based thinking is a basic element of human thinking; it has been used in many fields.In statistics, the Bernoulli distribution says if the probability of an event is , then the probability of its contrary is 1 − .Also, the pair (, −) represent opposite real numbers.It had been proved mathematically and experimentally that utilizing opposition in learning yields a more efficient algorithm than using only pure randomness [38].For this reason, EAs were enhanced by using OBL.When evaluating a solution x for a given problem, usually computing its opposite solution is performed together.This process will provide another chance for finding a candidate solution, which is closer to the global optimum.
Definition 1 (opposite points, see [34]).Let  be a real number defined on the interval [, ].The opposite number x is defined as follows: For  = 0 and  = 1, we have Similarly, the opposite number in a multidimensional search space can be defined.Let x = { 1 , . . .,   }  be a point in an -dimensional coordinate system with  1 , . . .,   ∈ R and   ∈ [  ,   ].The opposite point x is defined by its coordinates x 1 , . . ., x  where The above opposition-based concept can be embedded into MPBIL, which is governed by a probability update rule.In MPBIL, a probability matrix represents a binary population where   is a user-defined value in the range of [0, 1].When evaluating the value of   for a probability matrix P, it is possible to simultaneously compute its opposite   (or   ) for the opposite probability matrix (OP), which will provide better population variety in a PBIL search process.For OMPBIL, an initial probability matrix P and its opposition OP are sized /2 ×   .A binary population is generated using both matrices.When updating the th row of the matrix by using (3), the th row of the opposite probability matrix is also updated using   and the th row of P as where   = 1 −   .Two binary subpopulations with   / members based upon to the th row of the probability matrix and its opposite values obtained from ( 9) are then created.Other computational steps are similar to MPBIL.The outline of the OMPBIL algorithm, which includes the oppositionbased concept, is shown in Algorithm 3.

Numerical Experiment
The purpose of this research is to find a new internal structure for the morphing wings, which can adjust the aeroelastic characteristics by applying external torque actuation.This can be classified as a twisting morphing wing concept.Firstly, we start with finding the internal structural layout of a morphing aircraft wing three times a problem by employing an optimization process to synthesize the unconventional internal layout of the morphing aircraft wing.Then, we investigate aeroelastic response due to an actuating torque at a particular wing angle of attack to maintain or change the altitude of flight.

Synthesizing the Internal Structure of a Morphing Aircraft
Wing.For the first purpose of this study, an initial design domain was a simple rectangular unswept and untapered wing-box.In order to laterally verify the proposed design approach, four design problems are posed, while OMPBIL and MPBIL are used to solve the design problems three optimization runs for all test problems.The four sets of design variables or test cases are presented in The aeroelastic constraints require that the flutter and divergence speeds must be greater or equal to 150 m/s to avoid the aircraft structural failure at a desired flight speed range (100 m/s).Buckling and stress constraints are also considered.The void threshold is set to be 0.0005 m.
The parameters used for all wing-box models are given in Table 3.The wing is made of aluminum with Young modulus  = 70 GPa, Poisson's ratio ] = 0.34, density  = 2700 kg/m and yield stress   = 100 MPa.Aerodynamic analysis for static and dynamic aeroelasticity is carried out by using the vortex ring panel method [7][8][9][10].The wing lifting surface has 8 × 5 panels, while its wake has 8 × 30 panels.Flutter analysis is computed through the discrete-time aeroelastic model using a reduced-order aerodynamic model with static correction [39].Buckling factor and stress are computed considering both aerodynamic force and actuator torque.The first load (lift force) is computed at a wing angle of attack 5 ∘ , while the free stream velocity is 150 m/s.The actuator torques applied to twist a wing to increase or decrease wing lift is ±60 N-m (positive for clockwise torque).The forces are applied at the first rib as shown in Figure 5.All design constraints are computed where the prestressed effects due to external actuation are taken into account.For more details of each objective and constraint, see our previous studies [8][9][10].
The multidisciplinary optimization procedure is the combination of a multiobjective evolutionary algorithm, aerodynamic analysis, aeroelastic analysis, and finite element analysis.The vertex ring method for aerodynamic analysis is coded in MATLAB, while ANSYS software is used for FEA.Both optimization algorithms (MPBIL and OMPBIL) are also coded using the MATLAB language.Structural analysis, aerodynamic analysis, and an optimization algorithm are connected using MATLAB.
For performance assessment of the propose method, the conventional MPBIL and the new OMPBIL are employed to solve all of the design problems where the population size is set to be 50, the total number of generations is set to be 400, the external Pareto archive size is 50, and the number of binary bits for each design variable is 4. The learning rate (  ) is generated randomly in the interval [0.4,0.6].The mutation probability and mutation shift are 0.1 and 0.2, respectively.These parameters settings were successfully used with MPBIL in Bureerat and Sriworamas [29], Bureerat and Srisomporn [30], Kunakote and Bureerat [11], and Sleesongsom and Bureerat [8,9].Design constraints are handled by means of the nondominated sorting scheme [40].The procedure is terminated when reaching the maximum generation number.

Aircraft
Trimming.For this study, the aeroelastic characteristics and mechanical parameters need to be studied in steady level flight to ensure that the new internal structure can be used as a morphing aircraft wing.Normally, motion of an aircraft depends upon forces acting on an aircraft such as lift (), drag (), thrust (), and weight ().Both lift and drag can be found by aerodynamic analysis in which lift is perpendicular, while drag is parallel to the free stream direction.Lift is created by the wing to carry the weight of the aircraft whilst drag is a resistance of the aircraft motion, which is against the thrust force produced by an engine or propeller.At a steady level flight, with thrust vector aligned with the direction of flight, the static equilibrium can be expressed as follows: Equation (10) shows lift is balanced by weight in steady level flight.This equation has been used to study how much the angle of attack for an aircraft is trimmed to maintain a steady level flight if an aircraft is disturbed by changes in the structural stiffness with moving of rotating spars [9].For our work, ( 10) is used to find the angle of attack to maintain the aircraft at steady level flight.At this angle of attack, if an aircraft is disturbed by input torque (control input) to change the aircraft to other steady level flight, this implies that the new internal structure can be used as the morphing aircraft wing.To find an angle of attack at a flight speed and weight of aircraft, the bisection method is employed in Algorithm 4.
The searching procedure for a required angle of attack to maintain the aircraft in steady level flight is terminated when  2 is equal to the weight of aircraft (Equation ( 10)) at the root of the angle of attack  2 .

Unconventional Structures of a Morphing Aircraft Wing.
Nondominated fronts and the search history (front hypervolume versus generation numbers) of all design problems obtained from the best runs of both MPBIL and OMPBIL are shown in Figures 8 and 9, respectively.In Figure 8, approximate surfaces are used to represent the Pareto fronts as it is difficult to display a 3D front.The search history shows  that OMPBIL outperforms MPBIL before the halfway point of searching except for the case of F2 where the two methods interchangeably outperform each other and OMPBIL won at the end.The fronts obtained using OMPBIL are slightly changed for the last 10 iterations.Therefore, we assume that the optimizer converge to the global optimum.The comparative Pareto fronts are somewhat difficult to observe.As a result, a quantitative study, the use of a hypervolume indicator, is conducted to compare the quality of the Pareto fronts.This indicator determines an area (for two objectives) or volume (for more than two objectives) covered by a particular Pareto front with respect to a predefined reference point [41].The hypervolumes of the fronts of all test problems for all optimization runs are given in Table 4 where the reference point for computing the hypervolumes is set to be {0, 6.0 kg, −1}  .This table reveals that the fronts for all case studies obtained from using OMPBIL are better (higher average hypervolume) than those obtained from using a conventional MPBIL.The simultaneous partial topology, shape, and sizing design approach gives better results than the approach without nodal position design variables as the average hypervolume of F4 is higher than F3 and that of F2 is higher than F1.In addition, the results from using a multigrid approach are far superior to those from using single resolution of ground segments (F1 versus F3 and F2 versus F4).Some unconventional internal structures for the third and fourth design problems (first run) are selected using an even Pareto filter technique [42].The selected structures of F3 with their aerodynamic (a.a.) and elastic axes (e.a.) are shown in Figure 10 whereas their lower and upper skin thickness distributions are given in Figure 11.The elastic axis can be thought of as a beam representing a wing structure.When the wing is subject to torsion, it will rotate with respect to this axis.For all wings in Figure 10, the elastic axis is backward behind the aerodynamic axis, which means the wing under aerodynamic loads will have an increased angle of attack and consequently higher value of lift effectiveness.Moreover, some of the wings have forward swept elastic axes, which means they behave like a forward swept wing (usually have high lift effectiveness) although their wing shape are unswept.The structures have various internal layouts and skin thickness distributions.The optimum structures are mostly from the ground segments with medium (Figure 6 (2) element resolution (Figure 6(a)) is not a good choice for F1 and F2.However, in practice, a designer will never know, which resolution is most suitable; therefore, employing a multigrid approach is advantageous.Figure 12 shows the selected wing structures of F4 where their skin thickness distributions are illustrated in Figure 13.Similarly to the case of F3 in Figure 10, all the wings in Figure 12 have their elastic axes backward behind their aerodynamic axes, while some of them behave like a forward swept wing.It is shown that all the structures are from the low ground element resolution.This shows that when adding nodal positions to the design problems the proper ground elements resolution is different from the F3 case.The use of a multigrid approach automatically helps finding the suitable resolution during the optimization process.

Results of Aircraft Trimming.
To study the performance of the new morphing aircraft wing, a steady level flight at a flight speed of 100 m/s is defined and the mass of the whole aircraft is set to weigh 25 kg for simplicity.The wingbox (4) illustrated in Figure 12 is selected for this study.From the iterative procedure, it is found that an appropriate angle of attack is 5 ∘ .Next, the aeroelastic characteristics and mechanical parameters need to be examined at this angle of attack when an input torque is applied.An input torque in the range of [−60, 60] N-m is applied to the wing at the first wing rib by means of couple as shown in Figure 5.The aeroelastic characteristics and mechanical parameters of the unconventional wing due to the various input torque values are given in Table 5.The relationship between the input torques and the aeroelastic parameters of the wing is illustrated in ( (1) (4) Figure 14, while the buckling factor and the maximum von Misses stress on the wings are displayed in Figure 15.The lift effectiveness of the unconventional wing ranges between 0.967 and 1.045.This shows that the new unconventional wing can change the lift effectiveness over a wider range compared to the conventional two-spar wing in Sleesongsom et al. [9].The lift effectiveness of the unconventional wings, in cases of unloaded conditions, is 1.006.When increasing the input torque in the positive direction, the lift effectiveness of the wing is increased whereas it is decreased with negative torques.The divergence speed of the unconventional wing has maximum values when no external forces are applied and it is slightly reduced as the torque magnitudes become larger.The flutter speed of the wing has a similar trend to the divergence speed as its value decreases with increasing input torque.Buckling and stress on the wing are varied in such a way that the higher external force magnitude will cause more damage to the structure.The buckling factor is reduced, while the maximum stress is increased.

Conclusions and Discussion
Two significant contributions of this paper are the improved version of multiobjective PBIL based on an oppositionbased approach and the use of a multigrid ground element strategy to synthesize a wing structure.Firstly, the multigrid design approach for optimization of an unconventional morphing wing is proposed.The new design strategy is the   5.
procedure for simultaneous partial topology, sizing, and/or shape optimization, which uses multiple ground segment resolutions.This multigrid approach is more efficient than using single-resolution ground elements in the sense that the suitable grid resolution is automatically detected and used in one optimization run.The simultaneous topology, shape, and sizing design gives better results than design without shape design variables.The approximate Pareto front gives a variety of morphing wing internal layouts and skin thickness distributions.Rather than guessing a resolution of ground elements as in the case of using one resolution, an optimizer will automatically find the proper resolutions for design problems with the multigrid approach.The proposed design approach can be extended for practical aircraft design.Secondly, opposite-based multiobjective populationbased incremental learning is proposed to be compared with the original multiobjective population-based incremental learning.The multiobjective design problems of morphing wing structure with single-grid and multigrid design variables are then posed and tackled by both OMPBIL and MPBIL.The optimum results show that the proposed optimizer is superior to the traditional MPBIL.In addition, the results show that using OMPBIL in combination with a multigrid design approach is the best design strategy.The proposed concept is said to be acceptable as an alternative aircraft morphing wing design.
The proposed ideas can be extended real aircraft wing and then realized and used in some type of aerial vehicles.For a small unmanned aerial vehicle, the obtained complicated wing structures can be manufactured by means of threedimensional printing such as the Southampton University laser sintered aircraft [43,44].Our future work is to apply this design concept to a real wing, which is made of high-grade aluminum or carbon fiber composite.That means higher lift effectiveness and structural strength can be expected.Furthermore, rather than generating simple ground elements as used in this paper, using chaos generation [45,46] or complex networks [47,48] for generating unconventional wing internal topologies can be a potential means to enhance the wing performance.Moreover, it is also possible to use chaos for improving EA search performance.

Figure 1 :Figure 2 : 6 ×
Figure 1: Wing-box (a) ground segments without nodal position and (b) an internal structural layout of a wing-box.

Figure 3 :Figure 4 : 6 ×
Figure 3: Ground segments with nodal position variables (e) and an internal structural layout of a wing-box.
Initialization probability matrix P = [0.5]×  , external Pareto archive Pareto = {}.(1) Generate a binary population B from P (see Figure 8).(2) Decode the binary population to be x ×  and find the objective values f ×  .(3) Update Pareto by replacing it with the non-dominated solutions of a union set Pareto ∪ x. (4) If the number of members in Pareto exceeds the predefined archive size   , remove some of them by using an archiving technique.(5) If the termination criterion is fulfilled, stop the procedure.Otherwise, go to step 6: (6) Update P. (6.1)For  = 1 to . (6.1.1)Select  0 binary solutions from Pareto randomly.(6.1.2)Generate   using (4).(6.1.3)Update the th row of P by using (3).(6.1.4)Generate rand ∈ [0, 1] a uniform random number.(6.1.5)If rand < the predefined mutation probability, update the th row of P using (5).(6.2) Next .(7) Go to step 1.

Figure 10 :
Figure 10: Internal structural layouts of the selected wing-boxes in Figure 8(c).

Figure 12 :
Figure 12: Internal structural layouts of the selected wing-boxes in Figure 8(d).

Table 1 :
Probability vectors and corresponding binary populations.

Table 2 :
Initialization Probability matrix P = [0.5]/2×,oppositeprobabilitymatrixOP= [0.5]/2×,externalParetoarchivePareto= {}.(1)Generate a binary population B from P.(2) Decode the binary population to be x ×  and find the objective values f ×  .(3)UpdateParetobyreplacing it with the non-dominated solutions of a union set Pareto ∪ x. (4) If the number of members in Pareto exceeds the predefined archive size   , remove some of them by using an archiving Compute   = 1 −   and generate the th row of the opposite probability matrix OP using (9).(6.2.5) Generate rand ∈ [0, 1] a uniform random number.(6.2.6)If rand < the predefined mutation probability, update the th row of P and OP using (5).(6.2.7) Generate binary subpopulations SB and SBO from the th row of P and OP respectively.Four design cases.

Table 3 :
3 , (1) Start with the initial interval of the root angle of attack [ 1 ,  3 ].It should be noted that the lift of this interval covers the given weight of an aircraft; otherwise, find a new interval.(2) Calculate lift,   1 and   3 at  1 and  3 respectively.(3) Find  2 = ( 1 +  3 )/2, then calculate   2 .(4) If   2 −  is sufficiently small for accuracy requirement, then the procedure is stopped.Otherwise, the new interval of angle of attack is determined.(5) If   2 < , the new interval is  1 =  2 .Set   1 =   2 and go to step 3. (6) Otherwise, the new interval is  3 =  2 .Set   3 =   2 and go to step 3. Wing-box parameters.

Table 5 :
Aeroelastic characteristics and mechanical parameters due to various input torques of wing-box (4) in Figure14at trimming angle of attack = 5 ∘ .