Sizing of a multirotor drone (ISAE)#

Written by Marc Budinger (INSA Toulouse), Scott Delbecq (ISAE-SUPAERO) and Félix Pollet (ISAE-SUPAERO), Toulouse, France.

The objective of this notebook is to select the best compromise of components (propeller, motor, ESC, battery) of a multi-rotor drone for given specifiations.

# Import libraries
import numpy as np
import openmdao.api as om

Introduction#

We have seen that even at the component level the set of equations of a sizing code can generate typical issues such :

  • Underconstrained set of equations: the lacking equations can come from additional scenarios, estimation models or additional sizing variable

  • Overconstrained equations often due to the selection of a component on multiple critera: the adding of over-sizing coefficients and constraints in the optimization problem can generally fix this issue

  • Algebraic loops often due to selection criteria requiring informations generally available after the selection

The presence of singularities is emphasized when assembling component sizing codes to build a sizing code for the whole system.

Concerning overconstraints components, we have here:

  • Brushless motors with multiple torque and voltage constraints (hover and transient vertical displacement)

Multiple algebraic loops appears in the sizing problem:

  • The thrust depends of the total mass which depend of components required for generating this thrust

More details in the setting up of sizing code can be found in the following paper:

Reysset, A., Budinger, M., & Maré, J. C. (2015). Computer-aided definition of sizing procedures and optimization problems of mechatronic systems. Concurrent Engineering, 23(4), 320-332.

Design variables#

Exercise 15

Summarize the design variables that have been identified previously.

Name

Unit

Description

\(\beta_{pro}\)

[-]

\(\beta_{pro}=pitch/diameter\) ratio to define the propeller

Constraints#

Exercise 16

Summarize the constraints that have been identified previously (satisfied constraint are positive).

Name

Unit

Description

cons_1

[-]

\(cons_1 = M_{total}-M_{total,real}\)

Objectives#

Workflow#

Exercise 17

Propose a workflow using an N2 or XDSM diagram (hand drawn) that represents the ordering and interconnections of the disciplines of the sizing code. You can use the visualization of the design graphs.

Sizing code#

Lets now implement the code to size the multirotor. For this purpose we use the MDAO framework OpenMDAO. To stay simple first, we will embed the sizing in a simple OpenMDAO ExplicitComponent which is not the most efficient way of solving the problem. The best practice would be to decompose the problem into several ExplicitComponent to make the code more modular and the problem solving more efficient. For this purpose you can refer to the following paper:

Delbecq, S., Budinger, M., & Reysset, A. (2020). Benchmarking of monolithic MDO formulations and derivative computation techniques using OpenMDAO. Structural and Multidisciplinary Optimization, 62(2), 645-666.

Specifications#

The first step is to provide the specifications (top-level requirements) for the drone.

Main specifications :

  • a load (video, control card) of mass \(M_{load}\).

  • an autonomy \(t_{hf}\) for the hover flight.

  • an acceleration to take off \(a_{to}\).

# SPECIFICATIONS

# Load
M_pay = 50.0  # [kg] load mass

# Acceleration during take off
a_to = 0.25 * 9.81  # [m/s**2] acceleration

# Autonomy
t_hov_spec = 25.0  # [min] time of hover flight

# MTOW
MTOW = 360.0  # [kg] maximal mass allowed

Architecture definition and design assumptions#

Then, we must provide the architecture definition and design assumptions for the models.

# ARCHITECTURE of the multi-rotor drone (4,6, 8 arms, ...)
N_arm = 4  # [-] number of arms
N_pro_arm = 1  # [-] number of propeller per arm (1 or 2)

# BATTERY AND ESC : reference parameters for scaling laws
# Ref : MK-quadro
M_bat_ref = 0.329  # [kg] mass
E_bat_ref = 220.0 * 3600.0 * 0.329  # [J]
C_bat_ref = 5  # [Ah] Capacity
I_bat_max_ref = 50 * C_bat_ref  # [A] max discharge current

# Ref : Turnigy K_Force 70HV
P_esc_ref = 3108.0  # [W] Power
M_esc_ref = 0.115  # [kg] Mass

# MOTOR : reference parameters for scaling laws
# Ref : AXI 5325/16 GOLD LINE
T_nom_mot_ref = 2.32  # [N*m] rated torque
T_max_mot_ref = 85.0 / 70.0 * T_nom_mot_ref  # [N*m] max torque
R_mot_ref = 0.03  # [ohm] resistance
M_mot_ref = 0.575  # [kg] mass
K_T_ref = 0.03  # [N*m/A] torque coefficient
T_mot_fr_ref = 0.03  # [N*m] friction torque (zero load, nominal speed)

# FRAME
sigma_max = (
    280e6 / 4.0
)  # [Pa] Composite max stress (2 reduction for dynamic, 2 reduction for stress concentration)
rho_s = 1700.0  # [kg/m**3] Volumic mass of aluminum

# PROPELLER
# Specifications
rho_air = 1.18  # [kg/m**3] Air density
ND_max = 105000.0 / 60.0 * 0.0254  # [Hz.m] Max speed limit (N.D max) for APC MR propellers

# Reference parameters for scaling laws
D_pro_ref = 11.0 * 0.0254  # [m] Reference propeller diameter
M_pro_ref = 0.53 * 0.0283  # [kg] Reference propeller mass

Optimization variables#

The next step is to give an initial value for the optimisation variables:

# Optimisation variables : initial values
beta_pro = 0.33  # pitch/diameter ratio of the propeller
k_os = 3.2  # over sizing coefficient on the load mass
k_ND = 1.2  # slow down propeller coef : ND = NDmax / k_ND
k_mot = 1.0  # over sizing coefficient on the motor torque
k_speed_mot = 1.2  # adaption winding coef on the motor speed
k_mb = 1.0  # ratio battery/load mass
k_vb = 1.0  # oversizing coefficient for voltage evaluation
k_D = 0.99  # aspect ratio D_in/D_out for the beam of the frame

Sizing code#

Now, complete the following sizing code with the missing design variables (inputs), constraints/objective (outputs) and equations:

class SizingCode(om.ExplicitComponent):
    """
    Sizing code of the multirotor UAV.
    """

    def setup(self):
        self.add_input("beta_pro", val=0.0)

        self.add_output("t_hov", val=0.0)
        self.add_output("M_total_real", val=0.0)
        self.add_output("cons_1", val=0.0)

    def setup_partials(self):
        # Finite difference all partials.
        self.declare_partials("*", "*", method="cs")

    def compute(self, inputs, outputs):
        beta_pro = inputs["beta_pro"]

        #% OBJECTIVES
        # ---
        t_hov = C_bat / I_bat_hov / 60.0  # [min] Hover time
        M_total_real = (M_esc + M_pro + M_mot) * N_pro + M_pay + M_bat + M_frame  # [kg] Total mass

        #% CONSTRAINTS
        cons_1 = M_total - M_total_real

        outputs["t_hov"] = t_hov
        outputs["M_total_real"] = M_total_real
        outputs["cons_1"] = cons_1

Now that the ExplicitComponent is defined we have to add it to a Group, itself added to a Problem.

prob = om.Problem()

group = om.Group()
group.add_subsystem("sizing_code", SizingCode(), promotes=["*"])

prob.model = group

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options["optimizer"] = "SLSQP"
prob.driver.options["maxiter"] = 100
prob.driver.options["tol"] = 1e-8

prob.model.add_design_var("beta_pro", lower=0.3, upper=0.6)

prob.model.add_constraint("cons_1", lower=0)

# prob.model.add_objective('t_hov', scaler=-1)
prob.model.add_objective("M_total_real", scaler=0.1)

# Ask OpenMDAO to finite-difference across the model to compute the gradients for the optimizer
prob.model.approx_totals()

prob.setup()

# Setup initial values
prob.set_val("beta_pro", beta_pro)


prob.run_driver()

print("Design variables")
print("beta_pro :", prob.get_val("beta_pro"))

print("Constraints")
print("cons_1 :", prob.get_val("cons_1"))

print("Objective")
print("t_hov: ", prob.get_val("t_hov"))
print("M_total_real: ", prob.get_val("M_total_real"))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/system.py:2674, in System._call_user_function(self, fname, protect_inputs, protect_outputs, protect_residuals)
   2673 try:
-> 2674     yield
   2675 except Exception:

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/explicitcomponent.py:292, in ExplicitComponent._compute_wrapper(self)
    291 else:
--> 292     self.compute(self._inputs, self._outputs)

Cell In[5], line 22, in SizingCode.compute(self, inputs, outputs)
     20 #% OBJECTIVES
     21 # ---
---> 22 t_hov = C_bat / I_bat_hov / 60.0  # [min] Hover time
     23 M_total_real = (M_esc + M_pro + M_mot) * N_pro + M_pay + M_bat + M_frame  # [kg] Total mass

NameError: name 'C_bat' is not defined

During handling of the above exception, another exception occurred:

NameError                                 Traceback (most recent call last)
Cell In[6], line 29
     25 # Setup initial values
     26 prob.set_val("beta_pro", beta_pro)
---> 29 prob.run_driver()
     31 print("Design variables")
     32 print("beta_pro :", prob.get_val("beta_pro"))

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/utils/hooks.py:131, in _hook_decorator.<locals>.execute_hooks(*args, **kwargs)
    129 def execute_hooks(*args, **kwargs):
    130     _run_hooks(pre_hooks, inst)
--> 131     ret = f(*args, **kwargs)
    132     _run_hooks(post_hooks, inst)
    133     return ret

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/problem.py:740, in Problem.run_driver(self, case_prefix, reset_iter_counts)
    738     else:
    739         with SaveOptResult(driver):
--> 740             return driver.run()
    742 finally:
    743     self._recording_iter.prefix = old_prefix

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/drivers/scipy_optimizer.py:281, in ScipyOptimizeDriver.run(self)
    279 # Initial Run
    280 with RecordingDebugging(self._get_name(), self.iter_count, self) as rec:
--> 281     model.run_solve_nonlinear()
    282     self.iter_count += 1
    284 self._con_cache = self.get_constraint_values()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/system.py:4558, in System.run_solve_nonlinear(self)
   4552 """
   4553 Compute outputs.
   4554 
   4555 This calls _solve_nonlinear, but with the model assumed to be in an unscaled state.
   4556 """
   4557 with self._scaled_context_all():
-> 4558     self._solve_nonlinear()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/group.py:3669, in Group._solve_nonlinear(self)
   3666 name = self.pathname if self.pathname else 'root'
   3668 with Recording(name + '._solve_nonlinear', self.iter_count, self):
-> 3669     self._nonlinear_solver._solve_with_cache_check()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/solvers/nonlinear/nonlinear_runonce.py:26, in NonlinearRunOnce._solve_with_cache_check(self)
     25 def _solve_with_cache_check(self):
---> 26     self.solve()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/solvers/nonlinear/nonlinear_runonce.py:45, in NonlinearRunOnce.solve(self)
     41             subsys._solve_nonlinear()
     43 # If this is not a parallel group, transfer for each subsystem just prior to running it.
     44 else:
---> 45     self._gs_iter()
     47 rec.abs = 0.0
     48 rec.rel = 0.0

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/solvers/solver.py:803, in NonlinearSolver._gs_iter(self)
    801 if subsys._is_local:
    802     try:
--> 803         subsys._solve_nonlinear()
    804     except AnalysisError as err:
    805         if 'reraise_child_analysiserror' not in self.options or \
    806                 self.options['reraise_child_analysiserror']:

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/explicitcomponent.py:318, in ExplicitComponent._solve_nonlinear(self)
    316 with self._unscaled_context(outputs=[self._outputs], residuals=[self._residuals]):
    317     self._residuals.set_val(0.0)
--> 318     self._compute_wrapper()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/explicitcomponent.py:292, in ExplicitComponent._compute_wrapper(self)
    289     self.compute(self._inputs, self._outputs,
    290                  self._discrete_inputs, self._discrete_outputs)
    291 else:
--> 292     self.compute(self._inputs, self._outputs)

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/contextlib.py:137, in _GeneratorContextManager.__exit__(self, typ, value, traceback)
    135     value = typ()
    136 try:
--> 137     self.gen.throw(typ, value, traceback)
    138 except StopIteration as exc:
    139     # Suppress StopIteration *unless* it's the same exception that
    140     # was passed to throw().  This prevents a StopIteration
    141     # raised inside the "with" statement from being suppressed.
    142     return exc is not value

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/system.py:2680, in System._call_user_function(self, fname, protect_inputs, protect_outputs, protect_residuals)
   2678         raise
   2679     else:
-> 2680         raise err_type(
   2681             f"{self.msginfo}: Error calling {fname}(), {err}").with_traceback(trace)
   2682 finally:
   2683     self._inputs.read_only = False

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/system.py:2674, in System._call_user_function(self, fname, protect_inputs, protect_outputs, protect_residuals)
   2671 self._residuals.read_only = protect_residuals
   2673 try:
-> 2674     yield
   2675 except Exception:
   2676     err_type, err, trace = sys.exc_info()

File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/openmdao/core/explicitcomponent.py:292, in ExplicitComponent._compute_wrapper(self)
    289     self.compute(self._inputs, self._outputs,
    290                  self._discrete_inputs, self._discrete_outputs)
    291 else:
--> 292     self.compute(self._inputs, self._outputs)

Cell In[5], line 22, in SizingCode.compute(self, inputs, outputs)
     18 beta_pro = inputs["beta_pro"]
     20 #% OBJECTIVES
     21 # ---
---> 22 t_hov = C_bat / I_bat_hov / 60.0  # [min] Hover time
     23 M_total_real = (M_esc + M_pro + M_mot) * N_pro + M_pay + M_bat + M_frame  # [kg] Total mass
     25 #% CONSTRAINTS

NameError: 'sizing_code' <class SizingCode>: Error calling compute(), name 'C_bat' is not defined