On Structural Parameter Optimization Method for Quad Tilt-Wing UAV Based on Indirect Size Estimation of Domain of Attraction

To deal with the flexibility in the electric aircraft design, this study develops a structural parameter optimization method. The quad tilt-wing aircraft exhibits rotary-wing flight and fixed-wing flight modes by tilting its wings. This study focuses on the domain of attraction, which is a class of state vectors that converge to a static equilibrium point after infinite time. The method indirectly estimates the size of domain of attraction of an aircraft by calculating the radius of a sphere that is encompassed by the domain. In addition, after investigating the relationship between change in structural parameters and the size of the domain of attraction, the structural parameters that maximize the size of the domain are surveyed. The effectiveness of the proposed method is demonstrated for the optimization of the quad tilt-wing aircraft. The quad tilt-wing aircraft exhibits rotary-wing flight and fixed-wing flight modes by tilting its wings. As the characteristics of the dynamics of the quad tilt-wing aircraft demonstrate a significant nonlinear change depending on the tilt angle of the wing, we propose a method to optimize the structural parameters for an easy control of the aircraft. Simulation results show that the optimization scheme identifies the parameters that enable the aircraft to perform agile control to address a disturbance moment. Further, although the same controller is used in the original and optimized aircraft, the optimized aircraft can handle a disturbance in a shorter time when compared with the original aircraft.


I. INTRODUCTION
Throughout recent years, the electric aircraft technology is widely recognized as the future of the aircraft industry [1], [2]. Full or hybrid electric aircraft will achieve significant reduction in fuel costs, cut in maintenance, and quitter takeoff-and landing. To make the aviation more environment friendly, the development of the electric aircraft is inevitable. The electric-based technologies for aircraft have been employed for unmanned aerial vehicles (UAVs) prior to passenger aircraft. Many of the UAVs operate from electric energy stored in batteries and can offer more efficient evaluation of the technologies. UAVs are expected to be widely used in many fields like observations, mappings, hobbies, and missions for disaster management may common application areas. For example, natural disaster monitoring [3], [4], UAV-assisted disaster management [5], [6] have been demonstrated. Also, thanks to the separation of energy and thrust The associate editor coordinating the review of this manuscript and approving it for publication was Rosario Pecora . generation, in general, electric UAVs have a lot of flexibility in the design [7]. A quad tilt-wing (QTW) UAV is a good example. The QTW-UAV has four rotors placed on wings that operate according to tilt mechanisms. By changing the tilt angles of front and rear wings, the QTW-UAVs change their dynamics and combine the characteristics of the rotary-wing and fixed-wing flight modes. When the two tilt angles are 0 • , the QTW-UAV cruises with a fixed-wing mode. When the angles are 90 • , the QTW-UAV cruises with a rotary-wing mode (i.e., behaves as a multicopter). Rotating wing and rotor as a single piece avoids the impinging of the propeller slipstream on wings that reduces the thrust in the hover [8].
One of the most important characteristics of QTW-UAVs is the significant non-linear change in dynamics depending on the tilt angles. Recently, to address the nonlinearly changing dynamics and develop an actual QTW-UAV, the flight control for the QTW has been studied for this decade. Mathematical modeling and vertical flight modeling were reported in 2012 [9]. Cetinsoy et al. presented an actual design of the QTW-UAV and its corresponding simulation and limited real flight test results [10]. An autonomous attitude controller that can be applied around rotary-wing mode QTW-UAV was developed [11]. Sato and Muraoka have addressed the QTW problem and shown the effectiveness of the gain-scheduled controllers to address the dynamic change in the aerodynamic characteristics on an actual QTW-UAV [12]- [14]. For The QTW-UAV, the controller design method has been continuously developed [15], [16]. In this study, a design method for structural parameters of the electric QTW-UAV is developed. Design methods for aircraft have been proposed in several styles. Wang et al. composed a method that consists of the preceding design of a basic flowfield and the following generation of an airplane [17]. Their method successfully improved the volume performance of waveriders. Alvarez et al. focused on the aerodynamic interaction induced by the use of multiple rotors operating in close proximity [18]. They mentioned that the computational cost required to calculate the interaction makes conventional numerical methods prohibitive for design space exploration and explored the capability of the viscous vortex particle method for the conceptual design of electrical vertical takeoff and landing aircraft. The area of the tilt-wings including QTW, as in Ref. [19], has been a growing area of research and development. The behavior of the wings and propellers are tightly coupled and thus the design of the stable QTW-UAV system is a challenging task.
In this study, to deal with the significant change in the flight mode and increase in the flexibility in the design, we focus on the relationship between the structural parameters of the QTW-UAV and domain of attraction (DOA). The DOA is a class of state vectors that converge to a static equilibrium point after infinite time variation. In many applications, finding a static equilibrium point is not sufficient to analyze the system, because a very small neighborhood of the equilibrium point should also be considered in practice. And thus, DOA is an important tool in the stability analysis, because the size of the DOA shows that how much the initial points can be far away from the equilibrium point and trajectories can still be converged. Since finding an exact DOA is a very difficult problem, alternative methods that estimate the shape or size of the DOA have been studied [20], [21]. A DOA size estimation technique has been employed for a system analysis [22], attitude tracking technique for spacecraft [23]. If the size of the DOA of a system is expanded, the system can handle high disturbance. Hara et al. designed some controllers for a QTW-UAV based on different control theories and proposed to use one that yields the largest size of the DOA [24]. The relationship between the size of the DOA and structural parameters was discussed in Ref. [25]. In the paper, they proposed a structure of an inverted pendulum that expands the size of the DOA. Fig. 1 shows the schematic sketch of the structural parameter optimization method developed in this study. In the method, a numerical optimization algorithm and DOA size estimation technique are combined to identify optimal parameters. The DOA size is estimated for a set of structural parameters, and the parameters are changed such that the estimated size of the DOA is expanded. Since the size of the DOA directly reflects how much a system can deal with the change in state variables from an equilibrium point, expanding the size directly improves the performance of the system. The relationship between the estimated DOA size and structural parameters is used to determine the following behavior of the optimization algorithm. By conducting the scheme iteratively, we can obtain a better set of structural parameters that ensures the expansion of the DOA, i.e., a larger stable region in the variable space. In the proposed method, the expansion of the DOA size is ensured at the end of every steps. The system with a larger DOA can achieve more agile response. In this study, structural parameters like the length and positions of the wings, and wing areas are surveyed in an optimal framework to maximize the DOA size of the QTW-UAV. Optimal shapes of the QTW-UVA are surveyed to identify the structural parameters that enable the system to be easily controlled.
The remainder of this paper is organized as follows: in Sec. II, the dynamical model for the QTW-UAV and methodology for the optimization that employs the indirect estimation of the DOA are described. Using the proposed method, the size of the DOA is indirectly estimated as the radius of a sphere that is encompassed by the DOA. This method can reduce the computational resources used for the calculation. Section III discusses certain optimization results. The relationship between the DOA size and structural parameters are investigated using parametric survey, and all the parameters are optimized using a numerical method. Section IV concludes this paper and discusses the scope for future studies.

A. DYNAMICAL MODEL FOR QUAD TILT-WING UAV
The dynamical model for the QTW-UAV is presented. In this study, the model is developed based on the model that was proposed in Ref. [26] by incorporating a non-linear aerodynamic force model. To study the design problem of the QTW-UAV in an optimization framework, structural parameters are considered as variables. Fig. 2 illustrates the coordinate systems and structural parameters of the QTW-UAV. In this study, the body-fixed coordinate system (x b , x b , x b ) and inertial coordinate system (x, y, z) are introduced. The length l tx1/2 are the distances between the centroid and front/rear wings. The length l ty is the distance between x b -axis and propeller, and ξ 1/2 are front/rear wing tilt angles, respectively. The non-linear equations for the transitional and rotational motions for the QTW-UAV in the body-fixed coordinate system (x b , y b , z b ) are given bẏ where V is the transitional velocity vector, ω is the angular velocity vector, C B/I X is the transform matrix from the inertial frame to the body-fixed frame, G is the gravitational acceleration vector, m is the mass of the QTW-UAV, F a is the aerodynamic force, α b is the angle between V and x b -axis, ξ is the tilt angle vector, F u is the control force, J is the inertia tensor, M a is the aerodynamic moment, M u is the control moment, and M d is the disturbance moment. Vectors and matrices in Eqs. (1) and (2) are given by Eqs. (3)- (6), as shown at the bottom of the next page: The position, velocity, and acceleration of the QTW-UAV in the inertial coordinate system are obtained using other transform matrices [26]. The forces and moments generated by the propellers and flaperons (F u and M u ) are separated from those generated by the aero dynamical effects (F a and M a ). They are given by Eqs. (7) and (8), as shown at the bottom of the next page, where F Lw1/Lw2 are front/rear wing lifts, F Dw1/Dw2 are drags, F Lb and F Db are lift and drag of the body, T 1−4 are thrusts generated by the propellers, and F f1−f4 are flaperon forces. The moments are given by Eqs. (9) and (10), as shown at the bottom of the next page, where ρ is the air density, S v is the vertical stabilizer area, C Lv is the lift coefficient of the vertical stabilizer, and M p1-p4 are propellers anti torques.
The angles of attack are shown in Fig. 3. For a QTW aircraft, the angle of attack should be calculated to deal with the change in the tilt angle of a wing. The lift and drag acting on an i-th wing are defined as where S w i is the i-th wing area. To address the change in the dynamical mode caused by the tilt mechanism, we introduce the following C Lw and C Dw models that depend on the angle of attack α i of the i-th wing as in Eqs. (13) and (14), as shown at the bottom of the page. In Eqs. (13) and (14), C L0 is the lift coefficient at zero angle of attack, C Dp is the minimum parasitic drag coefficient, e 0 is the Oswald efficiency factor, and AR is the aspect ratio.
Using AR, C Lα is given by In the model given by Eqs. (13) and (14), two wing models are smoothly connected at a stall angle α 0 using the following blending function: where the η (= 0.8·180/π ) is the slope parameter. By assuming C L0 = 0 and C Dp = 0.00361, the C Lw and C Dw curves over the change in the angle of attack are shown in Fig. 4. As shown in Fig. 4, two different models are smoothly connected at the stall angle α 0 = 15 • . The lift and drag forces exerted on the body are given by where S b is the representative area of the body, and C Lb and C Db are the body lift and drag coefficients, respectively. We assume that the lift and drag coefficients of the body C Lb and C Db are constant. The force generated by the flaperon is given by where q i and q pi represent the dynamic pressure caused by the UAV velocity and propeller, S f is the wing area of the flaperon, C Df is the drag coefficient of the flaperon, and ζ fi is the angle of the flaperon. The dynamic pressures are VOLUME 10, 2022 The propeller slipstream U pi is given by Unlike the main wings, the flaperon force is assumed to be proportional to the angle of attack. Thus, the flaperon angles are assumed to be |ζ i | ≤ 15 • . The position of the flaperons are behind each propeller. The thrust and anti-torque of the i-th propeller are given by where k p is the lift coefficient of the propeller, pi is the angular velocity of the propeller, and b p is the drag coefficient of the propeller. For simplicity, we ignored the inertia of the propeller in this study. Furthermore, we ignored the lift and drag of the vertical stabilizer [26].
In this study, the optimal shape for the QTW-UAV is searched. The variables for the structural optimization problem are, l tx1 , l tx2 , l ty , l wf , S w1 , S w2 , b 1 , and b 2 . Using these parameters, I xx , I yy , I zz , I xz , m, and AR are indirectly optimized. In this study, the moments of inertia I xx , I yy and I zz are given by Eqs. (25)- (27), as shown at the bottom of the page, where l bx is the length of the body along the x b -axis, and l wx is the length of the wing along the x b -axis. Further, the product of inertia is assumed to be I xz = 0.1I xx . Therefore, the mass of the QTW-UAV is given by m = ρ b ∆y∆zl bx body dydz + ρ w ∆y∆zl wx wing dydz (28) The value of the body and wing density ρ b and ρ w , respectively, which are determined by the material and structure, are assumed to be 58.77 kg/m 3 to equalize the mass m to the original mass of 1.2 kg [26].

B. ESTIMATION OF DOMAIN OF ATTRACTION
For a given locally asymptotically stable equilibrium point, the DOA is a set of initial states for which all trajectories starting from this region converge to an equilibrium point (see also Fig. 5). We assume a state vector x and non-linear function f (x) as If we express the DOA as R 0 , it can be written as If an initial condition of the system is inside the DOA, the trajectory converges to the static equilibrium point x eq . If an initial condition is outside the DOA, the trajectory diverges. These two regions can be divided by the limit cycle such that the trajectories converge to a cycle in case the system is stable. By calculating the shape and size of DOA of a system, we can directly predict the stability and safety margin of a system for an initial state.
The shape of a DOA for a given non-linear system can be estimated using the following steps. If we integrate all the non-linear equations of motion for a sufficient number of initial state vectors for sufficiently longer simulation time to identify if an initial state converges to the equilibrium point, we can directly investigate the shape of the DOA. As the accuracy of the DOA estimation of this method depends on the number and range of the initial condition and simulation time, it requires a huge computational resource. The sum of squares (SOS) relaxation method indirectly estimates the size of the DOA by calculating the radius of a n dimension ellipsoid encompassed by the DOA [27]. If the system can be modeled by the SOS form, we can calculate the radius of the ellipsoid. The size of the DOA is indirectly estimated as the radius of the ellipsoid. However, the SOS relaxation method requires to solve the semidefinite programming problem, and it is difficult to be applied for high-order systems, such as the QTW [28]. In this study, the size of the DOA is indirectly estimated by calculating the radius of a sphere that is encompassed by the DOA by combining the Monte Carlo method and golden-section search method. To reduce the number of numerical integration, the numerical integration of all the non-linear equations is performed for initial state vectors that are on the surface of a given sphere. The radius of the sphere r cur is updated according to the golden-section search rule until the value converges to its maximum value. The calculation flow of the proposed method is shown in Fig. 6. For an initial sphere radius r cur , N MC of the initial state are randomly generated under limitation of its norm, ∆|x| = r cur . Then, the numerical simulations for a set of initial state, r min , r max , and r cur are updated according to the golden-section rule. This process is iterated until the iteration number is equal to its upper limit of N GS . The proposed method obtains the radius in the DOA with a given probability that depends on the number of the Monte Carlo sample N MC , and thus, reduces the required computational resources when compared with those of the direct investigation of the DOA shape. In this study, N MC = 1000 is used to obtain reliable estimation of the DOA radius. The reliability of the DOA estimation in Fig. 7 depends on the number of the Monte-Carlo sample N MC . The probability that the N MC of sample are encompassed by the DOA is expressed using the binomial distribution. We assumed that a single sample exists inside the DOA with a probability of p. If all the samples converge at the equilibrium point, they are encompassed by the DOA with a probability of p N NC . Thus, if all the samples that exist on the r DOA of sphere converge to the equilibrium point, more than N MC p of samples are inside the DOA with 1 − p N MC of confidence. As we used N MC = 1000, more than 99% of samples exist inside the DOA with 99% of confidence. Although the proposed method cannot investigate the accurate shape of the DOA, it can ensure the existence of the set of stable states using smaller computation resources.

C. DOA ESTIMATION SCHEME FOR QTW-UAV
To estimate the DOA size of the QTW-UAV, the dynamical system is linearized at a given equilibrium point and input that are represented as x eq and u eq . The state vector of the system is x = [ϕ θ ψ U V W P Q R] T and the input vector is u = [ 1 2 3 4 ζ 1 ζ 2 ζ 3 ζ 4 ] T . Then, the linearized QTW-UAV system is expressed asẋ = Ax + Bu, where A = ∂ẋ/∂x| x=x eq ,u=u eq and B = ∂ẋ/∂u| x=x eq ,u=u eq , respectively. The propeller thrust that achieves the equilibrium point is analytically calculated such that the QTW-UAV horizontally flies with a constant speed (ẋ = U = v e m/s, z = W = 0 m/s). The linear quadratic regulator (LQR) is employed to design the controller. The quadratic cost function of the LQR is given by In Eq. (31), Q is a 9 × 9 and R is a 8 × 8 matrix. For these matrices, diagonal matrices diag {Q} = 1 and diag {R} = 1 are used. The radius of the sphere is determined as the norm of the disturbance in the angular velocity, ∆P, ∆Q, ∆R as r cur = ∆P 2 + ∆Q 2 + ∆R 2 . The disturbance is converged to a disturbance moment M d and substituted into Eq. (2). For the convergence condition, we used (δϕ, δθ, δψ) ≤ 0.01 rad [28], (δU , δV , δW ) ≤ 0.01 m/s, and (δP, δQ, δR) ≤ 0.01 rad/s. We further constrained the change rate of each variables such that they must be smaller than 1.0 × 10 −3 unit/s. The time required to analyze the convergence of a trajectory is set to be t conv = 30 s. For better accuracy of the DOA estimation, we extended the time for convergence used in Ref. [28]. Fig. 7 shows an example of the estimation of the radius of the sphere in the DOA of the QTW-UAV proposed in Ref. [26]. We assumed the common tilt angles to be ξ 1 = ξ 2 = 30 • and original parameters summarized in Table 1. In this paper, the lift of the main body is neglected, because it is small compared with the lift generated by the front and rear wings. As shown in Fig. 7, the DOA radius converges to its maximum value at approximately the 12-th iteration. The estimated value of the radius at 20-th iteration is r DOA = 1.087 rad/s. In Fig. 8, the DOA is visualized in the ∆P∆Q∆R-space. The samples contained in the ∆P∆Q∆R-space were numerically calculated for the visualizing the distribution of the DOA in the space, which is independent of the numerical scheme used by the proposed method. The color of the samples represent the time required VOLUME 10, 2022  to satisfy the convergence condition. The samples depicted with red squares represent the states that are outside the DOA. For clarity, the samples in the ∆P > 0 rad/s, ∆Q > 0 rad/s, and ∆R > 0 rad/s region are not plotted. As shown in Fig. 8, the sphere effectively estimates the smallest distance between the equilibrium point and edge of the DOA.

A. OPTIMIZATION PLANNING
Optimization settings used in this study are described as follows. As initial structural parameters of the QTW-UAV, we used the values that are reported in Ref [26]. In this study, for simplicity, S w1 = S w2 , l tx1 = l tx2 , and b 1 = b 2 are assumed. The values of the moment of inertia and product of inertia provided in Ref. [26] are not used but are redefined using Eq. (25). Further, the front and rear wings are titled by the same angle, ξ 1 = ξ 2 . We assume that four rotors can generate a thrust that is equal to 1.2 times the gravity force acting on the centroid of the original QTW-UAV (20% margin). As the mass of the original QTW-UAV is 1.2 kg, each rotor can generate 3.53 N. If the propeller thrust and lift force cannot cancel the gravitational force acting on the QTW-UAV, the system becomes unstable and the sphere radius becomes r DOA = 0 rad/s.

B. PARAMETRIC STUDY FOR l tx AND l ty
As the first step of the optimization, we conduct parametric studies for varying the position of the wings and propellers (l tx and l ty ). The survey ranges are 0.15 m ≤ l tx ≤ 0.6 m and 0.1 m ≤ l ty ≤ 0.85 m. The change in the r DOA that is caused by the change in the position of the front and rear wings and four rotors is investigated. Fig. 9 shows the results for ξ 1,2 = 30 • , 45 • , 60 • , and 75 • . The samples that are emphasized with black squares represent the results of the original l tx and l ty values. The size of the sphere r DOA changes depending on the position of the wings and rotors. The DOA size of the original QTW-UAV can be expanded for every tilt angle. Further, the size of the DOA is affected by the tilt angle. Generally, when the tilt angle is small, the size of the DOA becomes small. When the tilt angle is small, the force points in the z-axis depend on the lift generated by the wings. If the QTW-UAV is changing its flight mode from the rotary-wing mode to the fixed-wing mode, the flight speed U may not be sufficient to generate the lift and flaperon forces, and thus, r DOA can become small. The results show that the DOA size strongly depends on the change in the structural parameters.
To address the simultaneous optimization for several parameters, numerical methods should be employed.

C. STRUCTURAL PARAMETER OPTIMIZATION EMPLOYING SEQUENTIAL QUADRATIC PROGRAMMING METHOD
Similar to the previous parametric study, expanding the DOA size for smaller tilt angles is an efficient way to design a more stable QTW-UAV. Here, we show an example of the structural parameter optimization for ξ 1,2 = 30 • case by employing the numerical optimization scheme. To optimize the structural parameters of the QTW-UAV, the DOA estimation scheme is combined with the sequential quadratic programming (SQP) method [29]. Structural parameters that maximize the DOA radius for a given tilt angle ξ are surveyed using the SQP algorithm. The optimization variables are l tx1,2 , l ty , b 1,2 , and S w1,2 . The parameters l tx1,2 and l ty are chosen so that the propellers and wings are placed within the original body length. The problem is formulated as follows: Subject to:   The body length along x b -axis l bx = 1.3 m and propeller diameter d p = 0.2032 m are used [26]. Also, the body width w b is set to be 0.1 m.
The optimization results are summarized in Table 2.
In addition, change in the QTW-UAV shape is shown in Fig. 10. Fig. 11 shows is the top view of the original and optimized structures of the QTW-UAV. Fig. 10a shows the original QTW-UAV shape and Fig. 10b shows the optimized shape that employed the proposed method, respectively. As the optimized QTW-UAV yields 22.3% larger r DOA when compared with the original QTW-UAV, the stability of the QTW-UAV for ξ 1,2 = 30 • is enhanced using the proposed optimization scheme. It means that the proposed optimization scheme successfully added the safety margin of the QTW-UAV for the disturbance moment. As shown in Table 2, the span is 21.2% longer than that of the original QTW-UAV. Further, the distance between the centroid and rotary axis of the wing l tx and distance between the centroid and rotor l ty are expanded. Such changes in the parameters can enable the QTW-UAV to generate larger control moments. The wing areas become 31.8% smaller than those of the original QTW-UAV. As a  result, the mass of the optimized QTW-UAV is 22.5% smaller than that of the original QTW-UAV. Larger control moments and smaller weight m enable the QTW-UAV to perform an agile attitude control. Fig. 12 compares the time histories of the state variables (a-c) and transitional velocity in the inertial coordinate system (d) when the QTW-UAV faces a disturbance. The disturbance moments are added to change the angular velocities in the body-fixed frame as ∆P = ∆Q = ∆R = 0.5 rad/s. The continuous lines represent the values of the optimized QTW-UAV while the dashed lines represent the values of the original QTW-UAV. As shown in Fig. 12, the values of the optimized QTW-UAV converge faster than those of the original UAV.
We also investigated the change in the r DOA caused by the design of the LQR controller. The diagonal elements of Q and R of the QTW-UAV whose structural parameters were optimized for ξ 1/2 = 30 • were searched so that the sphere radius r DOA was maximized. In this scheme, the proposed DOA estimation technique and SQP were used in combination. The results were summarized in Fig. 13. When the tilt angle was 30 • , the sphere radius was expanded to  Fig. 13 shows, even the QTW-UAV whose shape was optimized for ξ 1/2 = 30 • yields larger r DOA for other tilt angles compared with the original QTW-UAV. As we stated, the size of the DOA directly reflects the safety margin of the system, and thus, the QTW-UAV can deal with larger disturbance moments. The proposed DOA estimation technique provides an efficient objective function in the QTW-UAV optimization problem. By searching the structural parameters and diagonal elements of Q and R in an optimal framework by maximizing the sphere radius r DOA , the stability of the QTW-UAV is more enhanced.

IV. CONCLUSION AND SCOPE FOR FUTURE WORK
In this study, we proposed and investigated a method to optimize the structural parameters for QTW-UAVs by indirectly estimating the size of the domain of attraction. The size of the domain of attraction is modeled by the radius of the sphere that is encompassed by the domain of attraction, and the structural parameters of the QTW-UAV that expand the radius are surveyed. The results showed that, using the proposed optimization system, the size of the domain of attraction can be actively expanded by changing the structural parameters. Further, by optimizing the matrices in the quadratic cost function of the controller, the QTW-UAVs whose structural parameters were optimized for a given tilt angle achieve a large size of the domain of attraction for other tilt angles. The proposed optimization system can deal with the new flexibility in the aircraft design. The method can suggest a possible advantage of optimizing structural parameters of an electric aircraft based on the expansion of the domain of attraction size. In terms of practicality, demonstration experiments should be conducted in the future. In 2008, he joined the Faculty of Nagoya University, Nagoya, where he is currently a Professor with the Department of Aerospace Engineering, Graduate School of Engineering. His current research interests include motion and vibration control of mechanical structures, nonstationary control methods, control problems of man-machine systems, and aerospace equipment.