Drone Conduct: A New Paradigm in Flight Controls and Redundancy for Multi-Rotor Aerial Vehicles
Arinze N.N. Eze CAMC, MRAeS, MA, MSc
Ex-Rolls Royce | Ex-Airbus | Ex-Bayer | Intensifying mechanisation on farm fields with the utilisation of aerial (and ground) robotics
KEYWORDS
eVTOL, Drone, Aerospace, Aerial Robotics, Situational, Signal-Based Modelling, Control Allocation, System Modeller, Fuel Cell, Hydrogen, Battery, Energy Systems, Computer-Aided Engineering (CAE), Computer-Aided Analysis (CAA), Model-Based Design, Dynamic Modelling, Guidance, Navigation
ABSTRACT
High-count rotor-propelled eVTOL is designed to overcome significant challenges at low and very low altitudes. These challenges stem from the underactuated nature of current drone designs, which leads to complex control dynamics and instability, particularly during near-ground flights where most accidents occur. To address these issues, the eVTOL features an over-actuated design with 76 rotors and employs innovative DaLI control algorithms. These algorithms distribute control signals across redundant actuators to balance loads without causing instability or over-stressing components during flight.
Fundamentally, our four apriori art claims are as follows:
Our design takes advantage of a larger spread of discloading area and load distribution from a structural design point of view. Unlike other drones and eVTOLs with either high speed & low torque or high torque & low speed, the flight agility of our system is attributed to the design for the smallest possible moment of inertia the rotors can achieve with variational torque assignment based on dynamic inversion incremental compute of the DaLI inner loop. We are also patenting the realigning algorithm (DaLI) that mimics ‘flex’, ‘twist’ and ‘stiff’ morphological states of bird in flight embedded into the control commands. Inspired by avian somatosensory system that actively sums up the perpendicular vectors of the body frame’s CG to the relative position of its gravitational force for continuous stability via thrust vectoring, we are confident this will provide a new view point of aerial vehicle stability and flight control.
INTRODUCTION
Across the globe, crippling traffic in many cities’ circuitous roads has created a need for a path to decongestion. Admittedly, this problem has multiple implementable solutions in the developed world from speed trains, mass transit programs to tunneled/ underground passes.? However, emerging countries have no path to leapfrog 30 yrs of conventional infrastructure of any kind of global standard to build on without hefty price tags of US $4-10? trillion.? Much like the ideas around underground tunnels but without the high infrastructure costs, utilising the airspace with short-haul eVTOLs provides a relief valve for the 2D over-land traffic by operating in the less intrusive third dimension.
The availability of batteries with high energy density, the distributed propulsion technologies developed from drone-like transport systems opens up possibilities air vehicles/ taxis reducing time spent in traffic and opportunities to declutter major city roads. The 76-rotor propelled air vehicle presented in this article takes advantage of its technical redundancies, the minimal infrastructure requirements to the starting and landing spots of an eVTOL. In this article, we will show that the high volume of parts does not pose a high level of complexity in the control system. Neither does it drive high energy cost, memory/ storage nor lag of any sort.? Its real disadvantage of?low energy efficiency specific to stationary cruise flight is compensated by the transformation from electric takeoff to hydrogen.
Powered by a fuel cell and a battery stack, the hybrid powertrain combines the main advantages of the two power sources. The battery can deliver high currents, quick response to load changes and has a high volumetric power density. The fuel cell has higher specific energy densities (stored kWh/kg) which leads to reasonable safety margins like flight deviations, ‘return-to-home' and other critical flight requirements. Due to the latent behavior of fuel cell as a power source, the architecture in this article is for the PEM to deliver energy primarily for stationary cruise-flight. However, it doubles as a recharging energy source for the battery after takeoff guaranteeing energy surplus for landing.? It is without saying, that this high-pressure hydrogen energy delivered by the fuel cell significantly reduces size and weight of the battery required for the short haul flight envelope of the eVTOL.
CONCEPT DESIGN
In Simulink, the AriO Drone library (ADL) is developed for the polyrotor-styled nesting of the rotors in each quadrant for efficient command assignments. The Hytern library is to model the fuel cell and synergise its role with the battery management system. Models from the Wolfram System Modeller and the Simulink complement the eVTOL overall model. With the help of Mathematica and MATLAB practical computations are validated for the choice of airframe, loads, /noise and perturbations, battery pack size and fuel cell stack.
Partitioning the auxiliary battery with a switch function to accommodate recharging the main unit with the risk of disrupting flight power from the rotor system.? The assumption is the main (partitioned) battery recharging occurs during the cruise phase within the flight envelope. Other calculated assumptions are as follows:
Vertical Take-off?
The hovering force for constant hovering equals the gravitational force. Specific values like gravity is measured in m/s^2 and mass in kg.
m = Quantity[1.5, "Kilograms"];
g = Quantity[9.81, "Meters"/"Seconds"^2];
For a mass of 1.5 kg, as an example, calculated at hover state with a force:
F_h = m * g;
result = F_h == m * g // N
We assume the power required for hovering the eVTOL at altitude to be:
F_h = Quantity[1000, "Newtons"]; (* Example hover force *)
v_c = Quantity[5, "Meters"/"Seconds"]; (* Example climb velocity *)
η_p = 0.8; (* Example propeller efficiency *)
For a hover force of 1000 N, climb velocity of 5 m/s, and propeller efficiency of 0.8 produces a potential power of:
P_pot = (F_h * v_c) / η_p;
result = P_pot == (F_h * v_c) / η_p // N
The rotor area of all rotors of the eVTOL is calculated based on the number of rotors and specific values like diameter.
n = 76; (* Number of rotors *)
d = Quantity[1, "Inches"]; (* Diameter of each rotor *)
In our case of 76 rotors with 17-inch diameter each, the total rotor area is calculated by summing up the computed value of the rotor area:
A = n * Pi * (d/2)^2;
result = A // N
Also, the power to add the potential energy during the climbing phase is accounted for by the Hover force, Air density, Rotor area and, a constant k:
kInt = 1.2; (* Example constant k_int *)
Fh = Quantity[1000, "Newtons"]; (* Hover force F_h *)
rho = Quantity[1.225, "Kilograms"/"Meters"^3]; (* Air density rho *)
A = Quantity[1.5, "Meters"^2]; (* Rotor area A *)
Let's take the k_int = 1.2, hover force of 1000 N, air density of 1.225 kg/m^3, and rotor area of 1.5 m^2: The power required to add the potential energy during climbing is calculated as:
Ph = kInt * (Fh^(3/2)) / (2 * Sqrt[rho * A]);
result = Ph // N
Then the total power for the climbing phase (P_h) can be estimated:
kInt = 1.2; (* Example constant k_int *)
Fh = Quantity[1000, "Newtons"]; (* Hover force F_h *)
rho = Quantity[1.225, "Kilograms"/"Meters"^3]; (* Air density rho *)
A = Quantity[1.5, "Meters"^2]; (* Rotor area A *)
Now, calculate the power required for hovering.
Ph = kInt * (Fh^(3/2)) / (2 * Sqrt[rho * A]);
We can define specific values for the potential power (P_pot) as:
vc = Quantity[5, "Meters"/"Seconds"]; (* Climb velocity v_c *)
etaP = 0.8; (* Propeller efficiency η_p *)
Then calculate the potential power
Ppot = (Fh * vc) / etaP;
For the given values, the total power for the climbing phase is calculated as:
Pc = Ph + Ppot;
result = Pc // N
Stationary Hover
The propulsion thrust force to keep the aerial vehicle in hover position depends on its mass and the lift-to-drag ratio calculated as:
(* Define specific values *)
m = Quantity[1.5, "Kilograms"]; (* Mass of the vehicle *)
g = Quantity[9.81, "Meters"/"Seconds"^2]; (* Gravitational acceleration *)
E = 15; (* Lift-to-drag ratio *)
For a mass of 1.5 kg and a lift-to-drag ratio of 15, the propulsion thrust force can be calculated as follows:
Ft = (m * g) / E;
result = Ft // N
Such force should be multiplied by the velocity to get the aerodynamic power. Dividing the propulsion thrust force by the propeller efficiency results in the mechanical power on the drive shaft.
Ft = Quantity[200, "Newtons"]; (* Propulsion thrust force F_t *)
vt = Quantity[20, "Meters"/"Seconds"]; (* Velocity v_t *)
etaP = 0.85; (* Propeller efficiency η_p *)
For a propulsion thrust force of 200 N, velocity of 20 m/s, and propeller efficiency of 0.85, we can compute the mechanical power on the drive shaft as:
Pt = (Ft * vt) / etaP;
result = Pt // N
Ideally, this computational approach can be used to obtain all relevant parameterical values in sizing and designing the eVTOL unit.
Dynamic Inverses with Deafferented aerial Loci Increments (DaLI)
In vertical take-off, for instance, we can independently control the rotational speed of each rotor through engine control units by way of a selected group within 28 nested rotor models of the 76 individual rotors.?
These units ensure the necessary torque is applied to maintain the rotor speed at a specific reference point, a value assigned by our Deafferented dynamic inversion algorithm within the flight controller inner feedback.?Like in this vertical takeoff example, each propulsor contributes to the total thrust and/ torque desired, this active incremental inversion also applies for targeted altitude; pitch and forward flight, bank to given angles and yaw command and multiple combinations. These inputs are the guiding parameters for us to navigate and control our multi-rotor craft.?And here is the initial nonlinear state vector (x) equation for the open loop plant dynamics:
To define the state vector x and the reference rotor speeds Ω_ref:
x[t_] := {x1[t], x2[t], x3[t]} (* Example state vector, x1, x2, x3 are functions of time *)
OmegaRef[t_] := {omega1[t], omega2[t], omega3[t]} (* Example reference rotor speeds *)
Then define the nonlinear function f:
f[x_, OmegaRef_] := {
-x[[1]] + x[[2]] * OmegaRef[[1]],
x[[1]] - x[[3]] * OmegaRef[[2]],
x[[2]] + x[[3]] * OmegaRef[[3]]
} (* Example nonlinear function, customize as needed *)
The initial nonlinear state vector differential equation for the open loop plant dynamics will be:
stateVectorEquation = D[x[t], t] == f[x[t], OmegaRef[t]]
Essentially, by the time the vehicle reaches the desired altitude in our example, the state feedback for high torque and moderately higher speed and inadvertently highest temperature of the drivetrain is reached then tilting for propelling the vehicle in the forward direction commences.?This ramp up, as a handover process, happens whilst the lead rotors simultaneously wind down contributing low to zero power.?This is a welcomed redundancy with our 76 rotors providing over-actuation linearly related inputs and outputs we’ve defined with the matrix within our DaLI algorithm. According to the state equation defined earlier, input vector is "Ωref ∈ ?76"
We'll define the reference rotor speeds Ω_ref as a list of 76 elements. The reference rotor speeds (Ω_ref) for the multi-rotor craft are:
OmegaRef = Array[Subscript[Ω, ref, #] &, 76];
Our poly rotor grouping or nesting approach suggests for thrust coefficient and torque coefficient to constant values obtained from isolated single-rotor analysis but they change with the flight condition for more accurate aerodynamic response. So, by intermittent throttling of the nested rotor groups, the vehicle maintains temperature, control and stability.? However, at stationary phase of cruise flight, altitude is maintained by distributing command inputs based on our control assignment matrices producing a non-square matrix. To define the non-square matrix (R) and the vector (D):
R = Array[Subscript[r, #1, #2] &, {76, 76}]; (* Example 76x76 matrix *)
D = Array[Subscript[d, #] &, 76]; (* Example vector with 76 elements *)
The calculated reference rotor speeds (Ω_ref') for stationary cruise flight will be computed as:
OmegaRefPrime = Inverse[R].D;
result = OmegaRefPrime // MatrixForm
Nesting Inverse with Incremental Deafferentation
The measurement or estimates of states associated with high order modelling as the 76 multi-inflow counts from our eVTOL with lead-lag, torsion and dynamics is impractical.?Hence, the order reduction on the pitch and roll of the nested rotor-thrust demand from the ‘flex’, ‘twist’ and ‘stiff’ morphologies of the propulsion rings is derived from individual perturbations. Instead of determining the total vectors directly, the required control increment at a given moment with respect to the conditions of the aerial vehicle in the instant of time-step just before, is computed.?The linearised outcome is then dynamically wrapped and minimised.
Defining the matrices and vectors:
领英推荐
Q = Array[Subscript[q, #1, #2] &, {76, 76}]; (* Example 76x76 matrix Q *)
R = Array[Subscript[r, #1, #2] &, {76, 76}]; (* Example 76x76 matrix R *)
F = Array[Subscript[f, #1, #2] &, {76, 76}]; (* Example 76x76 matrix F *)
OmegaRefPrimeD = Array[Subscript[ω, ref, d, #] &, 76]; (* Desired speed changes vector with 76 elements *)
D = Array[Subscript[d, #] &, 76]; (* Desired thrust/torque vector with 76 elements *)
To calculate the adjusted reference rotor speeds (Ω'_ref) as follows:
OmegaRefPrime = Inverse[Q + Transpose[R].F.R].(Q.OmegaRefPrimeD + Transpose[R].F.D);
result = OmegaRefPrime // MatrixForm
Compute Overlapping due to Deafferented Signals
The challenge with saturation or overwriting of feedback will be mitigated by a combined effort of modifying the inverse algorithms and constructing weight matrices for the control increments within the specified time horizons
For the matrices and vectors defined:
Q = Array[Subscript[q, #1, #2] &, {76, 76}]; (* Example 76x76 matrix Q *)
Rmod = Array[Subscript[r, mod, #1, #2] &, {76, 76}]; (* Example 76x76 modified matrix R_mod *)
F = Array[Subscript[f, #1, #2] &, {76, 76}]; (* Example 76x76 matrix F *)
OmegaRefPrimeD = Array[Subscript[ω, ref, d, #] &, 76]; (* Desired speed changes vector with 76 elements *)
D = Array[Subscript[d, #] &, 76]; (* Desired thrust/torque vector with 76 elements *)
R = Array[Subscript[r, #1, #2] &, {76, 76}]; (* Example 76x76 matrix R *)
c = Array[Subscript[c, #] &, 76]; (* Vector c with 76 elements *)
Now, the calculations for the adjusted reference rotor speeds (Ω'_ref) are:
OmegaRefPrime = Inverse[Q + Transpose[Rmod].F.Rmod].(Q.OmegaRefPrimeD + Transpose[Rmod].F.(D - R.c));
result = OmegaRefPrime // MatrixForm
De-fluttering Electric Motor Dynamics
High perturbations are experienced with the four quadrants of the flight vehicle. The center of the 4 propulsion rings are connected to the center of the chassis and kevlar tension struts connect the outer ring to the undercarriage. Alleviating the unsteady loads from the CW and CCW rotations requires a harmonics decomposition algorithm.?For this, the lengthToCenter (arm length) and angle (angle of the ring1 with respect to the positive x direction) is provided as a parametrical input for the cancellation mechanism to determine the 0th to the selected higher-harmonics of the rigid state and operational state of the eVTOL body frame
1. De-flutter Matrix R∈R4×(76?f): - To de-flutter matrix, define the number of motor states after excluding f motor states:
f = 2; (* Example value for f *)
motorStates = 76 - f;
Then, define the de-flutter matrix (R) with dimensions 4 x (76 - f) as:
R = Array[Subscript[r, #1, #2] &, {4, motorStates}];
2. Input Weight Matrix Q∈R(76?f)×(76?f)
We can define the input weight matrix (Q) with dimensions (76 - f) x (76 - f) to get:
Q = Array[Subscript[q, #1, #2] &, {motorStates, motorStates}];
3. Desired Virtual Vector Ωref?d′∈R(76?f)
For the Desired Virtual Vector, define Ω'_ref_d with dimensions (76 - f) using:
OmegaRefPrimeD = Array[Subscript[ω, ref, d, #] &, motorStates];
4. Input Offset Vector c∈R(76?f)
And finally, for the input offset vector (c) with dimensions (76 - f) compute with the following:
c = Array[Subscript[c, #] &, motorStates];
Energy Supply
The model takes energy in the DC voltage-Battery input from the battery, convert to AC form and injects controlled amounts at a preset power factor regulated from the Battery Management System. As the fuel cell is not capable of fast answers on current cycles, the converter keeps the fuel cell regulated via the main switching inductor, and output capacitor ??. The required current is set according to the SOC of the battery, keeping it around 85-90%. Additional energy is required for cooling is simplified to an increased hydrogen consumption of the fuel cell.
The? DC bus subsystem also takes input signals from the throttle input off of the controller. The total output current is summed up and sent to the battery management system as? the desired load current. The seven groups of the 28 poly rotor subsets will then have in ring 1 and ring 3 a positive voltage applied whereas for ring 2 and ring 4 the polarity is reversed. This gives clockwise rotation to the assigned rotors in ring 1 and ring 3 and counter clockwise rotation in ring 2 and ring 4 For the model, voltage and current of a single cell are multiplied with the number of cells connected in series and in parallel. This simple model depends on the state of charge.? We are in preliminary designs of a simplified representation of leakage inductance in Δ ? ?? output isolation transformers and/or harmonic filters inferring the angular position per-sample terminal voltage measurements.
Battery and BMS
A battery pack consists of many single cells, connected both serial and parallel. With the calculated powers, a cell type and the number of cells can be defined. The number of cells depends on the chosen voltage level and the maximum power needed during the climbing phase. For this application, a voltage level of U = 29.6 V was assumed. The chosen cell type is a Lithium-Ion Celera 15 Ah cell with NiMnO2/graphite electrodes. Its voltage of about 1.44V per single cell this compute is to get the desired battery pack voltage U and the voltage of a single cell U_cell:
U = Quantity[400, "Volts"]; (* Example desired battery pack voltage *)
Ucell = Quantity[3.7, "Volts"]; (* Example voltage of a single cell *)
The number of cells (n_serial) needed in series to achieve the desired battery pack voltage is easily calculated here as:
nSerial = U / Ucell;
result = nSerial // N
In order to avoid excessive aging of the cells, the maximum current drawn from the cell should be kept below a 5C-rate. This corresponds to a current that would discharge the nominal capacity of a cell within 1 of an hour. To compute the maximum current for the chosen battery cell firstly, we define the cell capacity Q_cell in ampere-hours.
Qcell = Quantity[15, "AmpereHours"]; (* Cell capacity *)
We then define the discharge rate (5h^-1):
dischargeRate = Quantity[5, "1/Hours"]; (* Discharge rate *)
The maximum current a Lithium-Ion Celera 15 Ah cell can safely deliver is calculated as:
imaxCell = Qcell * dischargeRate;
result = imaxCell // UnitConvert
The maximum of the total current imax can be estimated from the lift-off power. The battery has to deliver most of the necessary power, the fuel cell is assumed to deliver its maximum power of about 85 kW.
To compute, first define the total lift-off power P and the maximum power provided by the fuel cell P_FC_max:
P = Quantity[500, "Kilowatts"]; (* Example total lift-off power *)
PFCmax = Quantity[200, "Kilowatts"]; (* Example maximum power provided by the fuel cell *)
Also, define the battery pack voltage U:
U = Quantity[400, "Volts"]; (* Example battery pack voltage *)
The maximum total current (i_max) the battery must supply will be determined as follows:
imax = (P - PFCmax) / U;
result = imax // UnitConvert
Dividing the total current by the current tolerated by one cell row yields about 6, resulting in 100 single cells in the battery pack.(SOC) and the current. Lookup tables parameterize the components of a 1RC circuit and an additional resistor, a simple electric equivalent circuit model for a cell. The lookup tables were created using measurement data from simulation.
Fuel Cell
The PEM fuel cell design is modeled within the same electric circuit. The model from our Hytern Libraries complement with buses and sensors that provide electrical power via voltage and current signals to the aforementioned main signal bus. The fuel cell voltage is described as a difference between the open circuit voltage and spikes in the potentiometer due to activation and ohmic losses. The influence of partial pressures and temperature is considered negligible. So, to define the open-circuit voltage E_OCV:
Eocv = Quantity[1.2, "Volts"]; (* Example open-circuit voltage *)
Followed by defining the activation losses η_AC:
etaAC = Quantity[0.05, "Volts"]; (* Example activation losses *)
We then define the current i and the internal resistance (R_i) as:
i = Quantity[100, "Amperes"]; (* Example current *)
Ri = Quantity[0.01, "Ohms"]; (* Example internal resistance *)
The actual voltage (E_FC) delivered by the fuel cell is calculated as:
Efc = Eocv - etaAC - i * Ri;
result = Efc // UnitConvert
Tafel equations are used to compute the current spikes due to activation ηAC and resistance is used for all low currents the polarization.
We define the constants:
R = Quantity[8.314, "Joules"/("Moles"*"Kelvin")]; (* Universal gas constant *)
T = Quantity[298, "Kelvin"]; (* Temperature *)
n = 1; (* Number of electrons transferred *)
F = Quantity[96485, "Coulombs"/"Moles"]; (* Faraday's constant *)
N = 2; (* Example constant N *)
A = Quantity[0.5, "Volts"]; (* Example constant A *)
i0 = Quantity[0.01, "Amperes"]; (* Exchange current density *)
The we define the actual current (i):
i = Quantity[0.005, "Amperes"]; (* Example current *)
The activation overpotential (η_AC) for the given current is defined as the Tafel equation:
etaAC = Piecewise[
{
{i / i0 * (R * T) / (n * F), i < i0},
{N * A * Log[i / i0], True}
}
];
result = etaAC // UnitConvert
The inner resistance Ri and the parameters for the Tafel equation are extracted from the polarization curve provided by the manufacturer.
If successful, DaLI (Deafferented aerial Locus Imbeds) to encode dimensionally high input and output inference at the ESC and motor edge of UAVs/ eVTOLs will make real-time abrupt flight modulation signals feasible.? As a compact, implantable device integrated into the ESCs of the UAV, the DaLI microprocessor unit, which is programmed at the flight control inner feedback loop level will be a failsafe mechanism – in the event of a master flight controller losing procedural control. The collaborative algorithm of DaLI enables it to seamlessly take over and maintain control of the aerial vehicle.
OBJECTIVE CONCLUSION
The goal of this article was to show the conceptual design of a 76-rotor propelled eVTOL (electric vertical take-off and landing) aircraft using a 1D system model. It was created in the combined analyses of signal-based modeling of MATLAB/ Simulink and the equation-based object-oriented modeling language of Wolfram’s System Modeler. By overlapping models within each of the individual modeling languages a robust and variable analyses of the vehicle’s capabilities was realised.
The all-electric lower speed and enhanced torque direct-drive motors for the 76-rotor propelled eVTOL will explore a combination of Li-ion batteries and PEM (Proton Exchange Membrane) fuel cell as its hybrid drivetrain. The fuel cell delivers the energy for the stationary cruise flight. From basic parameters that are known for the eVTOL, a preliminary choice for battery pack and fuel cell stack is made. Multiple subsystems are modelled using custom made Simulink signal blocks while a coordinated set of Modelica libraries are forged within Wolfram System Modeller. While Simulink unlocks customised libraries for more accurate sensor signal output analyses for the control plant, the Modelica code is part of System Modeler in the Wolfram package that gathers design parameters from mission simulations. Combining the platform serves as the source of truth providing the design parameters to later design stages of the eVTOL model including CFD simulations, airframe construction, load/noise analysis, costing and systems design. Gathered knowledge from downstream design process circles back to the system model to refine the assumptions. The system model can further be used to integrate the flight control algorithms into autopilot hardware and conduct SiL (software in the loop) tests.
Thanks for taking the time to read through another one of what started as a stream of consciousness into an article inspired by my stint and blogpost at the Wolfram Institute on Flight Control Agent Trained in a Reinforcement Learning Environment
Other References
I’m Arinze N. N. Eze, an aerospace engineer who's been around the block more times than a confused pigeon. After bagging over 20 years of experience, I think I’ve mastered making robots fly like caffeinated birds and destructive testing by crashing them. I had employment at big-shot companies like Airbus and Bayer, Exhibit A. I'm sure my colleagues would say I spend my days there trying to convince them that flying cars are just around the corner.
I convinced a bunch of smart engineers from the Singapore University of Technology and Design on the ALKONOST project to create aerial robots mimicking birds—because nothing says "innovation" like robots on tons and tons of energy. The H.A.W.K. project generates energy to power a village or a really big toaster, and IPAMS project aims to send drones to Mars, because Earth was getting boring.
I’ve presented papers at UNESCO, appeared on CNN and Al Jazeera, and collected awards like they’re going out of style. My career is about making the world better, one flying robot at a time.
Cool project, well done!
业务拓展与销售 | 产品管理与创新 | 创业与初创公司 | 沃顿商学院 MBA | 新加坡永久居民 (PR)
4 个月Great work, team!