</a>
</a>

# Sizing of a multi-rotor drone

*Written by Marc Budinger, Aitor Ochotorena (INSA Toulouse) and Scott Delbecq (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.

**Scipy** and **math** packages will be used for this notebook for the optimization algorithms of python.

**widgets** provides an interactive layout to display the results.


```python
import scipy
import scipy.optimize

from math import pi
from math import sqrt
from math import sin

import numpy as np
import pandas as pd

import ipywidgets as widgets
from ipywidgets import interactive
from IPython.display import display

pd.options.display.float_format = '{:,.2f}'.format
```

## Contents
1. [Sizing the propeller for vertical climb speed ](#section_1)
2. [Architecture defintion and design assumptions](#section_2)
3. [Sizing code](#section_3)
4. [Optimization problem](#section_4)

<a id='section_1'></a>

## Sizing propeller to vertical climb speed scenario

This exercise targets to minimize the propeller torque for vertical climb speed scenario. This Notebook shows the role of the design variables as well as the constraints in an optimization problem. 


Let us say that $Q_{pro,cl}$ is the climb torque to minimize. 

For the vertical climb, advance ratio ($J=\frac{V_{cl}}{n_{cl}   D_{pro}}$) is a key parameter to estimate the vertical speed. 

Given that $D_{pro}$ and $V_{cl}$ are known from the specifications of the problem, we face a problem with one equation and two undetermined parameters.

**STRATEGY:** One of these parameters is going to be defined within a design range (e.g. J) and thus, the rest of calculations are done. This particular case is called **underconstraint singularity**.

Finally, estimated parameter is compared to the real value through the constraint. 

The sizing code is defined here in a function which can give:
- an evaluation of the objective: here minimize the torque in climb scenario
- an evaluation of the constraints: 



<math>\begin{align}
&\underset{\mathbf{x}}{\operatorname{minimize}}& & Q_{pro,cl}(\mathbf{x}) \\
&\operatorname{subject\ to}
& & 5\% \leq  V_{cl}-J n_{cl}   D_{pro} \leq 0 \\
\end{align}</math>

where it is expected that the product $J n_{cl}   D_{pro}$ calculated does not differ more than a 5\%.

Here is an non-exhaustive XDSM diagram of the whole multirotor sizing code:



More details in the setting up of sizing code can be found in the  [following paper](https://www.researchgate.net/profile/Marc_Budinger/publication/277933677_Computer-aided_definition_of_sizing_procedures_and_optimization_problems_of_mechatronic_systems/links/55969de508ae793d137c7ea5/Computer-aided-definition-of-sizing-procedures-and-optimization-problems-of-mechatronic-systems.pdf):  

> 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.


<a id='section_2'></a>

## Architecture defintion and design assumptions


```python
Mtotal = 1.35 # [kg] mass of the drone 

rho_air=1.18# [kg/m^3] Air density

C_D=1.3 #[-]drag coefficient

pitch=4.5*0.0254#[m] propeller pitch

Dpro=10*0.0254#[m] propeller diameter

k_maxthrust= 3 #Ratio max thrust-hover [-]

A_top=0.045 #[m^2] top surface

N_arm = 4 # [-] number of arm

N_pro_arm = 2 # [-] number of propeller per arm (1 or 2)

N_pro = N_pro_arm * N_arm # [-] Propellers number

V_cl = 5 # [m/s] max vertical speed

bounds=[(0.01,0.5)]#J
```

<a id='section_3'></a>

## Sizing code


```python
# -----------------------
# sizing code
# -----------------------
# inputs: 
# - param: optimisation variables vector (reduction ratio, oversizing coefficient)
# - arg: selection of output  
# output: 
# - objective if arg='Obj', problem characteristics if arg='Prt', constraints other else

def SizingCode(param, arg):

# Design variables
# ---
    J = param[0] # advance ratio

# Hover, Climbing & Take-Off thrust 
# ---
    
    F_pro_hov=Mtotal*(9.81)/N_pro # [N] Thrust per propeller for hover
    
    F_pro_to=F_pro_hov*(k_maxthrust) # [N] Max Thrust per propeller   
    
    F_pro_cl=(Mtotal*9.81+0.5*rho_air*(C_D)*(A_top)*(V_cl)**2)/N_pro # [N] Thrust per propeller for climbing
    
# Propeller models
# ---

    # Statics APC MR
    
    beta=pitch/Dpro # [-] Ratio pitch-to-diameter
    
    C_t_sta=4.27e-02 + 1.44e-01 * beta # Thrust coef with T=C_T.rho.n^2.D^4
    C_p_sta=-1.48e-03 + 9.72e-02 * beta  # Power coef with P=C_p.rho.n^3.D^5
    
    # Dynamics APC MR

    C_t_dyn=0.02791-0.06543*J+0.11867*beta+0.27334*beta**2-0.28852*beta**3+0.02104*J**3-0.23504*J**2+0.18677*beta*J**2 # thrust coef for APC props in dynamics
    C_p_dyn=0.01813-0.06218*beta+0.00343*J+0.35712*beta**2-0.23774*beta**3+0.07549*beta*J-0.1235*J**2 # power coef for APC props in dynamics

# Propeller mass
# ---

    Dpro_ref=11*.0254 # [m] diameter
    Mpro_ref=0.53*0.0283 # [kg] mass

    Mpro=Mpro_ref*(Dpro/Dpro_ref)**3 # [kg] Propeller mass

#Rotational speed 
# ---

    n_pro_hover=sqrt(F_pro_hov/(C_t_sta*(rho_air)*Dpro**4)) # [Hz] hover speed
    Wpro_hover=n_pro_hover*2*3.14 # [rad/s] Propeller speed    

    n_pro_to=sqrt(F_pro_to/(C_t_sta*(rho_air)*Dpro**4)) # [Hz] climbing speed
    Wpro_to=n_pro_to*2*3.14 # [rad/s] Propeller speed

    n_pro_cl=sqrt(F_pro_cl/(C_t_dyn*(rho_air)*Dpro**4)) # [Hz] climbing speed
    Wpro_cl=n_pro_cl*2*3.14 # [rad/s] Propeller speed
        

# Propeller selection with take-off scenario
        
    Ppro_to=C_p_sta*(rho_air)*n_pro_to**3*Dpro**5# [W] Power per propeller
    Qpro_to=Ppro_to/Wpro_to # [N.m] Propeller torque

# Propeller torque& speed for hover
    
    Ppro_hover=C_p_sta*(rho_air)*n_pro_hover**3*Dpro**5# [W] Power per propeller
    Qpro_hover=Ppro_hover/Wpro_hover # [N.m] Propeller torque    

# Propeller torque& speed for climbing
   
    Ppro_cl=C_p_sta*(rho_air)*n_pro_to**3*Dpro**5# [W] Power per propeller
    Qpro_cl=Ppro_cl/Wpro_cl # [N.m] Propeller torque


# Objective and Constraints sum up
# ---
    Qpro_cl # minimize torque in climbing 

    constraints = [(-J*n_pro_cl*Dpro+V_cl),
                    (0.05+J*n_pro_cl*Dpro-V_cl)]

    # Run algorithm slsqp
    if arg == 'Obj':
        return Qpro_cl # for torque optimisation

     # Run algorithm differential evolution
    elif arg == 'ObjP':
        P = 0. # Penalisation nulle
        for C in constraints: 
            if (C < 0.): 
                P = P-1e9*C
        return Qpro_cl + P # for torque optimisation       

    elif arg=='Prt':

        col_names = ['Type', 'Name', 'Value', 'Unit', 'Comment']

        df = pd.DataFrame()
        df = df.append([{'Type': 'Optimization', 'Name': 'J', 'Value': J, 'Unit': '[-]', 'Comment': 'advance ratio '}])[col_names]
        df = df.append([{'Type': 'Optimization', 'Name': 'Const 0', 'Value': constraints[0], 'Unit': '[-]', 'Comment': '(V_cl-J*n_pro_cl*Dpro)'}])[col_names]
        df = df.append([{'Type': 'Optimization', 'Name': 'Const 1',  'Value': constraints[1], 'Unit': '[-]', 'Comment': '-V_cl+J*n_pro_cl*Dpro+0.05'}])[col_names]
        df = df.append([{'Type': 'Optimization', 'Name': 'Objective',  'Value': Qpro_cl, 'Unit': '[kg]', 'Comment': 'Torque vertical climb'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'beta', 'Value': beta, 'Unit': '[-]', 'Comment': 'Ratio pitch to diameter'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'pitch', 'Value': pitch, 'Unit': '[m]', 'Comment': 'Pitch'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'F_pro_to', 'Value': F_pro_to, 'Unit': '[N]', 'Comment': 'Thrust for 1 propeller during Take Off'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'F_pro_cl', 'Value': F_pro_cl, 'Unit': '[N]', 'Comment': 'Thrust for 1 propeller during Climbing'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'F_pro_hov', 'Value': F_pro_hov, 'Unit': '[N]', 'Comment': 'Thrust for 1 propeller during Hover'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'rho_air', 'Value': rho_air, 'Unit': '[kg/m^3]', 'Comment': 'Air density'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Dpro_ref', 'Value': Dpro_ref, 'Unit': '[m]', 'Comment': 'Reference propeller diameter'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'M_pro_ref', 'Value': Mpro_ref, 'Unit': '[kg]', 'Comment': 'Reference propeller mass'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'C_t_sta', 'Value': C_t_sta, 'Unit': '[-]', 'Comment': 'Static thrust coefficient of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'C_p_sta', 'Value': C_p_sta, 'Unit': '[-]', 'Comment': 'Static power coefficient of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'C_t_dyn', 'Value': C_t_dyn, 'Unit': '[-]', 'Comment': 'Dynamic thrust coefficient of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'C_p_dyn', 'Value': C_p_dyn, 'Unit': '[-]', 'Comment': 'Dynamic power coefficient of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'D_pro', 'Value': Dpro, 'Unit': '[m]', 'Comment': 'Diameter of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'n_pro_to', 'Value': n_pro_to, 'Unit': '[Hz]', 'Comment': 'Rev speed of the propeller during takeoff'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'n_pro_hov', 'Value': n_pro_hover, 'Unit': '[Hz]', 'Comment': 'Rev speed of the propeller during hover'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'n_pro_cl', 'Value': n_pro_cl, 'Unit': '[Hz]', 'Comment': 'Rev speed of the propeller during climb'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'P_pro_to', 'Value': Ppro_to, 'Unit': '[W]', 'Comment': 'Power on the mechanical shaft of the propeller during takeoff'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'P_pro_hov', 'Value': Ppro_hover, 'Unit': '[W]', 'Comment': 'Power on the mechanical shaft of the propeller during hover'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'P_pro_cl', 'Value': Ppro_cl, 'Unit': '[W]', 'Comment': 'Power on the mechanical shaft of the propeller during climb'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'M_pro', 'Value': Mpro, 'Unit': '[kg]', 'Comment': 'Mass of the propeller'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Omega_pro_to', 'Value': Wpro_to, 'Unit': '[rad/s]', 'Comment': 'Rev speed of the propeller during takeoff'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Omega_pro_hov', 'Value': Wpro_hover, 'Unit': '[rad/s]', 'Comment': 'Rev speed of the propeller during hover'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Omega_pro_cl', 'Value': Wpro_cl, 'Unit': '[rad/s]', 'Comment': 'Rev speed of the propeller during climb'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Qpro_hov', 'Value': Qpro_hover, 'Unit': '[N.m]', 'Comment': 'Torque on the mechanical shaft of the propeller during hover'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Qpro_to', 'Value': Qpro_to, 'Unit': '[N.m]', 'Comment': 'Torque on the mechanical shaft of the propeller during takeoff'}])[col_names]
        df = df.append([{'Type': 'Propeller', 'Name': 'Qpro_cl', 'Value': Qpro_cl, 'Unit': '[N.m]', 'Comment': 'Torque on the mechanical shaft of the propeller during climb'}])[col_names]
        df = df.append([{'Type': 'Specifications', 'Name': 'M_total', 'Value': Mtotal, 'Unit': '[kg]', 'Comment': 'Total mass'}])[col_names]
        df = df.append([{'Type': 'Specifications', 'Name': 'N_pro', 'Value': N_pro, 'Unit': '[-]', 'Comment': 'Number of propellers '}])[col_names]
        df = df.append([{'Type': 'Specifications', 'Name': 'N_pro_arm', 'Value': N_pro_arm, 'Unit': '[-]', 'Comment': 'Number of propellers per arm '}])[col_names]

        items = sorted(df['Type'].unique().tolist())
        def f(Type):
            return df[df['Type']==Type] 
        widgets.interact(f, Type=items)
        return f
             
    else:
        return constraints

```

<a id='section_4'></a>

## Optimization problem


We will now use the [optimization algorithms](https://docs.scipy.org/doc/scipy/reference/optimize.html) of the Scipy package to solve and optimize the configuration. We can use the SLSQP algorithm without explicit expression of the gradient (Jacobian) or a differential evolution optimization to have a more exhaustive calculation. A course on Multidisplinary Gradient optimization algorithms and gradient optimization algorithm is given [here](http://mdolab.engin.umich.edu/sites/default/files/Martins-MDO-course-notes.pdf):
> Joaquim R. R. A. Martins (2012). A Short Course on Multidisciplinary Design Optimization. University of Michigan


The first step is to give an initial value of optimisation variables:


```python
# Optimisation variables
J = .01 # advance ratio

# Vector of parameters
parameters = np.array([J])
```

We can print of the characterisitcs of the problem before optimization with the initial vector of optimization variables. Through the droplist menu, an interactive display of the parameters grouped by their type of variable is given:


```python
# Initial characteristics before optimization 
print("-----------------------------------------------")
print("Initial characteristics before optimization :")
SizingCode(parameters,'Prt')
print("-----------------------------------------------")
```

    -----------------------------------------------
    Initial characteristics before optimization :



    interactive(children=(Dropdown(description='Type', options=('Optimization', 'Propeller', 'Specifications'), va…


    -----------------------------------------------


Then we can solve the problem and print of the optimized solution:


```python
# Optimization with SLSQP algorithm
contrainte = lambda x: SizingCode(x, 'Const')
objectif = lambda x: SizingCode(x, 'Obj')
objectifP = lambda x: SizingCode(x, 'ObjP')

SLSQP = False # Optimization algorithm choice

if SLSQP == True:
    # SLSQP omptimisation
    result = scipy.optimize.fmin_slsqp(func=objectif, x0=parameters, 
                                   bounds=[(0.01,0.5)],
                                   f_ieqcons=contrainte, 
                                   iter=1500, acc=1e-12)
else:
    # Differential evolution omptimisation
    result = scipy.optimize.differential_evolution(func=objectifP,
                                   bounds=[(0.01,0.5)],
                                   tol=1e-12)

# Final characteristics after optimization 
print("-----------------------------------------------")
print("Final characteristics after optimization :")

if SLSQP == True:
    SizingCode(result,'Obj')
    SizingCode(result, 'Prt')
else:
    SizingCode(result.x,'Obj')
    SizingCode(result.x, 'Prt')
print("-----------------------------------------------")


```

    -----------------------------------------------
    Final characteristics after optimization :



    interactive(children=(Dropdown(description='Type', options=('Optimization', 'Propeller', 'Specifications'), va…


    -----------------------------------------------



```python
parameters.shape

```




    (1,)




```python

```

f797e4366cf3dc2628a2b90c3c10aff5b72b776f
ipynb
Jupyter Notebook
notebooks/.ipynb_checkpoints/08a_SizingCodeVerticalFlight-checkpoint.ipynb
aitorochotorena/multirotor-all
e0820a0fd32f852549ebde54425177156fdf40db
Apache-2.0
2021-01-28T09:57:42.000Z
2021-01-28T09:57:42.000Z
notebooks/.ipynb_checkpoints/08a_SizingCodeVerticalFlight-checkpoint.ipynb
SizingLab/droneapp-legacy
b7844810107051bc4bb1861ecedc71188d09e881
Apache-2.0
notebooks/.ipynb_checkpoints/08a_SizingCodeVerticalFlight-checkpoint.ipynb
SizingLab/droneapp-legacy
b7844810107051bc4bb1861ecedc71188d09e881
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES