Multicopter Design Optimization and Validation

This paper presents a method for optimizing the design of a multicopter unmanned aerial vehicle (UAV, also called multirotor or drone). In practice a set of datasheets is available to the designer for the various components such as battery pack, motor and propellers. The designer can not normally design the parameters of the actuator system freely, but is constrained to pick components based on available datasheets. The mixed-integer programming approach is well suited to design optimization in such cases when only a discrete set of components is available. The paper also includes an experimental section where the simulated dynamic responses of optimized designs are compared against the experimental results. The paper demonstrates that mixed-integer programming is well suited to design optimization of multicopter UAVs and that the modeling assumptions match well with the experimental validation

The interest in multicopters (also called multirotors or drones) has increased significantly in recent years, both in academic research and in commercial applications.
One example is the Prime Air multicopter by Amazon, see Fig. 1 and Amazon (2015).Other examples are the multicopters designed and developed to demonstrate acrobatic abilities, see for example Hehn and D'Andrea (2014) and the references therein.Despite the large number of published works, very little is published on design optimization of multicopters given an intended application and desired performance specifications.
In Magnussen et al. (2014) a concept for design optimization of a multicopter based on mixed-integer programming was presented.In the current paper the same concept is described in more detail giving all the necessary information to allow a design engineer to repeat the design procedure.In normal practice, a set of datasheets of components such as battery pack, motors and propellers is available to the designer.Variables such as thrust per propeller or motor torque can not be designed freely, but must be chosen from a fixed set of available datasheets.The mixed-integer programming framework solves optimization problems that fall within this category.Certain variables are constrained to be discrete (for example a Boolean variable indicating selection of a particular component represented by a datasheet) while others are continuous variables (for example total weight and flight time).As demonstrated in the paper, the mixed-integer linear programming framework also allows for approximation of nonlinear functions.A nonlinear function is approximated by a set of discrete regions of the y-axis, each represented by a linear function.The accuracy of the approximation can be adjusted by changing the total number of linearized regions.
Numerical methods for nonlinear optimization normally suffer from drawbacks such as long computation times and the fact that the methods can end up in undesirable sub-optimal local minima.Examples are Newton-Raphson gradient search methods and evolutionary-based search methods, such as the Complex method, see for example Whitney (1969) and Tyapin and Hovland (2009) where the optimization algorithm took several hours to converge towards an acceptable but sub-optimal solution.A linear program (LP) on the other hand, using the interior point method, can be solved in polynomial time, see Karmarkar (1984).The solution of an LP with a large number of variables can be found quickly using commercial solvers such as IBM CPLEX, see IBM (2015).The solution of a mixed-integer linear program (MILP) is built on an LP solver and techniques such as branchand-bound.The optimization solver formulated in this way can handle nonlinear function approximations and discrete design variables and returns a repeatable solution in relatively short time compared to an iterative nonlinear search algorithm where the solution can vary from run to run, depending on the initial conditions.It should be noted, however, that mixed-integer optimization problems are NP-hard, see for example Clausen (1999).
The paper also contains a section presenting experimental validation of the performance of the multicopter actuators.Both the actuator rise-time constants and the total flight time used in the design optimization procedure are compared with the same parameters estimated from real measurements using the various components available to the designer.The validation compares the performance of four propellers and three motors against the simulation model in three different test cases: I) Longest possible flight time, no payload, II) Longest possible flying time, 1.5kg payload and III) Fastest possible motor response, no payload.For the three test cases, optimized solutions were calculated for multicopters with 4, 6 and 8 actuators, giving a total of 9 potentially different designs.The results demonstrate that the simulation model matches well with the experiments.The experimental validation increases the reliability of the optimization results.
The paper is organized as follows: in Section 2 a dynamic model of a multicopter actuator is presented.Section 3 gives a detailed description of the design optimization procedure.Section 4 presents results from the experimental validation based on three different test cases.Discussion and conclusions are presented in Section 5.

Modeling
A Simulink model of the multicopter actuator is shown in Fig. 2. The model generates the simulated thrust response for the 9 different test cases used in this paper (test cases I, II and III with 4, 6 and 8 actuators).The propeller thrusts from the simulation model are later compared with experiments (force measurements using a load cell mounted underneath the actuator).The propeller thrust is calculated from eq. ( 1): where C t is the thrust coefficient, ρ is the density of air, n is the rotor speed (rev/sec) and D is the propeller diameter.The actuator power W and torque T are calculated from eqs. ( 2)-( 3): where C p is the power coefficient and w = 2πn is the rotor speed (rad/sec).The motor controller and actuator dynamics are given by the equations below: where J is the total inertia, V m , i m , τ m are the motor voltage, current and torque, K e , K t and R a are motor parameters and V i is the supply voltage from the motor controller.Mechanical friction is not included in the dynamic model in eq. ( 7), since this information is difficult to obtain from the manufacturer's datasheets.
Eqs. ( 1)-( 7) are implemented in the Simulink model in Fig. 2. By combining eqs.( 2)-( 7) and using the fact that for DC motors the motor torque and back emf constants are equal, that is, K t = K e , the following transfer function can be defined: where . Assuming constant parameters in eq. ( 8), the time-constant would be t r = JRa K 2 t +DωRa .Since D ω is a function of the propeller speed, this time-constant is an estimation.

Design Optimization
A quadratic program is defined as follows: where x ∈ R n is the state vector to be solved, Q ∈ R n×n is a positive semi-definite penalty matrix, g ∈ R n is a penalty vector, A ∈ R m×n is the constraint matrix while b ∈ R m is the constraint vector.When the penalty matrix Q = 0, eqs.( 9)-( 10) reduce to a linear program.When some elements are constrained to be integer variables, eqs.( 9)-( 10) are called mixed integer quadratic program (MIQP) or, when Q = 0, mixed integer linear program (MILP).By constraining the integer variables further to accept only Boolean values (0 or 1), logical constraints can be incorporated and linked with the continuous variables in the optimization problem.The following rules 1 and 4 are taken from Bemporad and Morari (1999) and Mignone (2002).Rules 19a and 19b are modified compared to Mignone (2002) (δ = 0 instead of δ = 1).In the rule definitions below the following notation is used: denotes a small, real, positive constant, typically the machine precision.
It should be noted that equality constraints of the type ax = b must be converted to two inequality constraints ax ≤ b and −ax ≤ −b to satisfy eq. ( 10).However, the solver used (CPLEX) allows specification of both equality and inequality constraints.
Product of Boolean variable and function of continuous variables.Product of several Boolean variables and function of continuous variables.
Rule 4: Product of several δ's and f (x) The scalar y ≥ 0 equals the nonlinear function g(•) of the scalar x ∈ [0, x max ].The set z approximates g(x) by using a set of straight line Only one of the Boolean variables in the set ).The nonlinear function approximation builds on rules 1, 19a and 19b in the equations below: Implication. [ with three regions (δ 12 , δ 13 and δ 14 ).
Fig. 3 illustrates how a nonlinear function can be approximated using a combination of discrete and continuous variables.In this example three regions δ 12 , δ 13 and δ 14 are chosen to approximate the nonlinear function g(x).The δ's are chosen to represent an equal spacing on the y-axis in the figure.In each region represented by δ i a linear approximation The objective function g T x in eq. ( 9) can be defined as shown in Table 1 depending on the application's requirement specifications.For payload transfer and package delivery, for example, the longest possi-

Criterion
Objective function Min.Time-constant g T x = x 16 Min.Power Table 1: Different candidate objective functions to be minimized/maximized depending on the requirement specifications.
Variables Description For the multirotor design optimization, the Boolean decision variables δ i are defined as shown in Table 2 and the continuous variables x i are defined as shown in Table 3.The complete mixed-integer state vector x introduced in eq. ( 9) is defined as follows: The indices i and j of discrete and continuous variables are summarized in Table 4. Three different MILPs are created.One with the number of actuators constrained to N a = 4, one with N a = 6 and one with N a = 8.The total number of discrete and continuous variables varies between the three programs and is equal to 1795, 2444 and 3091, respectively.To avoid multiplication of several free variables in the constraints, it was decided to create three separate MILPs rather than creating one optimization program with N a as a free variable.To find the optimal solution, the three MILPs are run separately and the solution with the lowest objective function value is chosen.
In the following, the complete set of constraints used in the multicopter design optimization is listed.The constraints make use of the rules 1, 4, 19a, 19b  Set of time-constants Distance from cg perp.to roll axis Distance from cg perp.to pitch axis Distance 2 from cg perp.to roll axis Distance 2 from cg perp.to pitch ax.
Set of actuator roll inertias Set of actuator pitch inertias Set of frame roll inertia Set of frame pitch inertia Set of frame lengths nonlinear function approximation defined earlier in this section.
Boolean variable constraints: Total weight and price: (m m,i δ i ) ( 16) Individual propeller thrust and n 2 : x . . .
x 16 = 96 i=49 x i (40) It should be noted here that the time-constants calculated in eqs.( 35) and ( 39) are estimates based on the linear transfer function in eq. ( 8).If real measurements of the time-constants were available for each combination of motor and propeller, these timeconstants would be used in x 49 • • • x 96 instead of the estimates based on the datasheet information.
Calculate square of roll & pitch lengths: These constraints use the nonlinear function approximation rules.The corresponding Boolean variables are x j2+1 = x 2 97 (42) . . .

Experimental Results
The purpose of the experiments presented in this section is to validate the models and assumptions used in the mixed-integer optimization presented in section 3. The experimental setup is shown in Fig. 4. The validation presented in this paper has chosen both the 63% rise time of the actuator thrust as well as the total flight time as the parameters to compare against the simulation model in Fig. 2. The rise times are measured directly from step responses in thrust, while the total flight time is estimated via the measured current when the step responses have reached steady-state.

Summary of Datasheets
The available datasheets are summarized in Table 5, 6 and 7.In total three motors (see Fig. 5), four propellers (see Fig. 6) and four batteries were available (the same as in the design optimization procedure in section 3, Table 2).Since all the batteries have a supply voltage of 11.1V and the experiments focused on validation of the 63% rise-time constant and flight time estimated via the measured current, the same battery could be used in all the tests, see Fig. 7.The selection of battery in the experiments did not have an impact on the estimated time-constants and total flight time.

Test Case I
In this test case the design was optimized for the longest possible flying time and no payload.As seen in Table 8, propeller 4 was chosen when the design was constrained to use 4 and 8 actuators, while propeller 3 was chosen with 6 actuators.Motor 3 was chosen with 4 and 6 actuators, while motor 2 was chosen with 8 actuators.Battery 3 was chosen regardless of the number of actuators being 4, 6 or 8. Battery 3 has the highest capacity (16000 mAh), but also the highest weight.
The selection of the battery with the highest capacity is not obvious.The solution with 4 actuators has both the longest flight time and the lowest price, so it is the prefered design in Test Case I.
Fig. 8 illustrates the experimental results from Test Case I.The black curves show the simulated thrust response when using the Simulink model in Fig. 2 with either 4, 6 or 8 actuators.The red curves show the measured thrust response under acceleration, while the green curves show the measured curves under deceleration.Note that the green curves are inverted about the steady-state value for easier comparison with the acceleration and the simulated response.The circles in the figure represent the values at the 63% rise time.The experimental results confirm that a linear model assumption is appropriate for the multicopter actuator, consisting of battery pack, controller, motor and propeller as shown in Fig. 4. Tables 11 and 12 summarize

Test Case II
In Test Case II the design was optimized for the longest possible flying time with a 1.5kg payload.As seen in Table 9, propeller 4 and battery 3 were chosen regardless of the number of actuators being 4,6 or 8. Motor 3 was chosen when the design was constrained to 4 and 6 actuators, while motor 1 was chosen with 8 actuators.The solution with 8 actuators has the longest flight time.However, the price with 8 actuators is significantly higher than the price for the solutions with 4 and 6 actuators.11 and 12 summarize the three different rise times (simulation, acceleration, deceleration).For Test Case II the simulated rise times are slightly faster than the measured values, ranging from 10 to 15ms faster.Overall, the match between the simulation model and the experiments is good.

Test Case III
In Test Case III the design was optimized for the fastest possible motor response and no payload.As seen in Table 10 propeller 1, motor 1 and battery 1 were chosen regardless of the number of actuators being 4,6 or 8. Since the time-constant (39ms) is the same for 4, 6 and 8 actuators, the natural choice is to use the solution with 4 actuators since the flight time is the highest Table 13 shows the differences between the timeconstants found from the experiments and the estimates based on the transfer function in eq. ( 8) and also used in the optimization, eqs.( 35), (39).The results are satisfactory with differences in the range -41 to +25ms.One error source is the time-constant assumption made in eq. ( 8) for a linear system.
The actual flight times can be estimated by dividing the battery capacity (Ah) by the measured current (A) when the step-responses have reached steady-state and by the number of actuators (N a ).Table 14 shows the estimated flight times vs. the (inverted) flight times calculated in eqs.( 72)-( 74).The standard deviation between measured and estimated flight times in Table 14 is 12.8%.Note that the battery voltage as stated in the datasheets is 11.1V, while the voltage when the battery is fully charged is more than 12V.In addition, the battery is not capable of keeping the voltage higher than 11.1V when discharged.Hence, both the estimated flight times and the ones found from experiments in Table 14 probably overstate the actual flight times slightly.

Conclusions
This paper has presented an optimization framework for multirotors based on mixed-integer programming.A simulation model of a multirotor actuator is presented and this model has been validated against experiments in three different test cases (longest possible flying time with or without payload, as well as fastest motor response).The experiments have used a step input in thrust and compared the 63% rise-time against the simulation model.The results are good with steptime differences between simulation and experiments in the range -21 to +42ms.The relatively small differences may be caused by unmodelled effects, such as motor friction and measurement delay.The difference between the estimated time-constants used in the design optimization and the time-constants estimated from the experiments have also been evaluated.The results are good with differences in the range -42 to +25ms.The total flight times have also been validated and have a standard deviation of 12.8% compared to the modeled flight times.
Overall, the results presented in this paper demonstrate that mixed-integer programming provides both a relatively accurate and efficient approach to multirotor design optimization from available datasheets.On an Intel i7-4770S 3.1GHz processor the different optimization problems were solved in typically 5-25 seconds using the IBM CPLEX solver.The experimental validation confirms that the modelling assumptions made in the design optimization formulation are reasonable and hence increase the confidence that the optimized designs actually meet the intended application's requirement specifications.

Figure 1 :
Figure 1: Example of multirotor UAV design by Amazon for transport and delivery of parcels.

Figure 4 :
Figure 4: Experimental setup for validation.A: Digital/analog IO card, B: Load cell amplifier, C: Current sensor, D: Battery pack, E: Motor controller, F: Propeller, G: Motor, H: Load cell.

Figure 5 :
Figure 5: Three different motors used in the design optimization and validation tests.

Figure 6 :
Figure 6: Four different propellers used in the design optimization and validation tests.

Figure 7 :
Figure 7: Battery used in the design validation tests.

Fig. 9
Fig.9illustrates the experimental results from Test Case II.The results for Test Case II are better than for Test Case I. Tables11 and 12summarize the three different rise times (simulation, acceleration, deceleration).For Test Case II the simulated rise times are slightly faster than the measured values, ranging from 10 to 15ms faster.Overall, the match between the simulation model and the experiments is good.

Table 2 :
Boolean decision variables.If δ = 1, then the corresponding component/region is selected.
If δ = 0, the component/region is not selected.ble flight time or smallest possible power consumption may be desirable objectives.For acrobatics, on the other hand, the designer may instead want to select the solution giving the smallest possible time-constant or smallest possible roll and pitch inertias.
and the

Table 3 :
Continuous variables used in the design optimization."cg perp.to" = "center of gravity perpendicular to".

Table 5 :
Datasheets: Motors.m m,i is the motor weight, D m,i is the motor diameter, S l,i ,S w,i are the shaft length and width, ρ m,i is the material density of the motor and P m,i is the motor price.

Table 6 :
Datasheets: Batteries.m b is the weight of the battery, while P b is the price.

Table 7 :
Datasheets: Propellers.m p,i is the weight of the propeller, p i is the pitch angle of the blade and P p,i is the propeller price.

Table 8 :
Optimization Results: Test Case I.

Table 9 :
Optimization Results: Test Case II.

Table 10 :
Optimization Results: Test Case III.

Table 11 :
Test T a,4 T d,4 T a,6 T d,6 T a,8 T d,8 Measured rise times (63%, in ms) for the different test cases and number of propellers.T a is for acceleration, T d is for deceleration.

Table 12 :
Simulated rise times (63%, in ms) for the different test cases and number of propellers.T i is the rise time with i propellers.∆T i is the time difference between T i and the average of T a,i and T d,i .

Table 13 :
Time-constants found from experiments (mean value of T a,i and T d,i in Table11) for the different test cases and number of propellers.∆Ti is the time difference (in ms) between T i and the time-constant estimates based on eq.(8).Test T e,4 T o,4 T e,6 T o,6 T e,8 T o,8 I 66.6 62.5 57.5 60.0 62.0 55.5 II 29.4 25.4 28.5 27.1 33.2 27.3 III 33.9 37.2 25.3 29.4 20.5 23.8

Table 14 :
Estimated flight times in minutes from the experiments T e,i vs. flight times calculated in the optimization T o,i for test cases I, IIand III and 4, 6 and 8 actuators.combinations of components from a set of available datasheets.Nonlinear functions can be approximated by using a combination of discrete and continuous variables.The designs can be optimized towards a set of different criteria, such as flight time, power consumption or dynamic performance, depending on the designer's preferences.Optimized designs with 4, 6 and 8 propellers are presented using real components available at, for example, HobbyKing (2015).