Formal Method Library for Propulsion System of FCS

Zhengpu Shi
 
Draft V1.0

Oct 25, 2020

Contents

1 Introduction
2 Propulsion System (PS)
 2.1 Common Module
  2.1.1 Definition of constants
  2.1.2 Definition of variables
  2.1.3 Variable Relation Graph (VRG)
  2.1.4 Definition of basic functions
  2.1.5 Definitions of basic inverse functions
  2.1.6 Definitions of extended functions
  2.1.7 Definitions of advanced extended functions
  2.1.8 Simplest Form of Functions (SFFs)
 2.2 Formal verification
  2.2.1 Customized automatic tactic
  2.2.2 Verification of equality of function transformation
 2.3 Example 1: Estimate the longest hovering time
  2.3.1 Analysis of the question
  2.3.2 General calculation step
  2.3.3 Direct calculate formulas
  2.3.4 Vefiry direct formulas
 2.4 Example 2: Estimate system effencicy at max throttle
  2.4.1 Analysis of the question
  2.4.2 General calculation step
 2.5 Example 3: Estimate maximum load and maximum pitch Angle
  2.5.1 Analysis of the question
  2.5.2 General calculation step
 2.6 Example 4: Estimate maximum flat flight speed and distance
  2.6.1 Analysis of the question
3 Summary

1 Introduction

This is one part of the FML4FCS project. For detail information, please refer the introduction of the whole project in this website.

2 Propulsion System (PS)

The modeling of Propulsion System (we called PS) is a basic condition for designing the flight control system.

We set up a module named Common Module which related to the general issues. Then we set up some modules named like Example* Module for several specific problems.

2.1 Common Module

2.1.1 Definition of constants

Constants refer to values that do not change after the aircraft is designed, such as propeller size. We have slightly simplified the aerodynamics, and set the flight altitude and ambient temperature to constants, so the air pressure and air density also become constants. The reason for this is that their changes during low-altitude flight are very small, almost negligible, and in order to highlight other more important variables, simplifying them to constants will make the entire model more concise.

Although the Coq system supports to using Unicode characters as identifier, we still use ASCII text represent an identify. Because other languages may not support this feature and we need to maintain naming consistency.

E N V I R O N M E N T
P R O P E L L E R
M O T O R
E S C
B A T T E R Y
Symbol Identifier Meaning Unit Remark
T 0 T_0 standard temperature K 273.15K
p0 p_0 standard atmospheric pressur Pa 101.325 × 105
ρ0 rho_0 standard atmospheric density kg∕m3 1.293
h h Altitude m
T t T_t Environment temperature ∘C
nr n_r Number of rotors piece
G G Total weight N
Iother I_other Other current A
p p atmospheric pressure Pa
ρ rho density kg∕m3
D p D_p Propeller diameter m
Hp H_p Propeller geometric pitch m
Bp B_p Number of blades piece
Gp G_p Propeller weight N
PPA PP_A Propeller parameter* - [1][p67]
PPε PP_epsilon Propeller parameter* - [1][p67]
PPλ PP_lambda Propeller parameter* - [1][p67]
PPζ PP_zeta Propeller parameter* - [1][p67]
PPe PP_e Propeller parameter* - [1][p67]
PPK0 PP_K0 Propeller parameter* - [1][p67]
PPα0 PP_alpha0 Propeller parameter* - [1][p67]
Cfd C_fd Zero-lift drag coefficient -
CT C_T Thrust coefficient -
CM C_M Torque coefficient -
Cd C_d Drag coefficient -
KV 0 K_V0 Nominal no-load motor KV r∕min∕V
Um0 U_m0 Nominal no-load motor voltage V
Im0 I_m0 Nominal no-load motor current A
ImMax I_mMax Maximum motor input current A
Rm R_m Motor resistance Ω
Gm G_m Motor weight N
KE K_E Back-electromotive force constant -
KT K_T Torque constant -
I eMax I_eMax Maximum ESC output current A
R e R_e ESC resistance Ω
G e G_e ESC weight N
R b R_b Battery resistance Ω
U b U_b Battery voltage V
C b C_b Battery capacity mAh
C min C_min Battery minimum capacity mAh
Kb K_b Maximum discharge rate C
Gb G_b Battery weight N

Notice that, some constants are composed by more basic constants. We give all the relationship of them.

Definition of p, see [1][p.66, 4.4]

( h )5.2561 p = p0 1- 0.0065 T-+-T-- 0 t
(2.1)

Definition of ρ, see [1][p.66, 4.3]

ρ = ---T0p---ρ p0(T0 + Tt) 0
(2.2)

Definition of CT, see [1][p.67, 4.6]

3 2 εarctanπHDpP-- α0 CT = 0.25π λ ζBpK0 ----πA-+-K0----
(2.3)

Definition of Cd, see [1][p.67, 4.7]

( Hp- )2 πAK20--εarctan-πDp --α0-- Cd = Cfd + e (πA + K0 )2
(2.4)

Definition of CM, see [1][p.67, 4.6]

C = -1-π2C ζ2λB2 M 8A d p
(2.5)

Definition of KE, see [1][p.80, 4.64]

KE = Um0---Im0Rm- KV 0Um0
(2.6)

Definition of KT, see [1][p.80, 4.65]

KT = 9.55KE
(2.7)

2.1.2 Definition of variables

Variable means that its value may changed at different or in different problems.

E N V I R O N M E N T
P R O P E L L E R
M O T O R
E S C
B A T T E R Y
Symbol Identifier Meaning Unit Remark
p p atmospheric pressure Pa
ρ rho density kg∕m3
Tb T_b Time of endurance min
V V Flight speed m∕s
V max V_max Maximum flight speed m∕s
Z Z Flight distance m
Zmax Z_max Maximum flight distance m
η eta System efficiency m 0˜1
Gmaxload G_maxload Maximum load N
θmax theta_max Maximum pitch angle rad
T T Propeller thrust N
M M Propeller torque N.m
N N Motor speed r∕min
E a E_a Back-electromotive force V
U m U_m Equivalent motor input voltage V
Im I_m Equivalent motor input current A
I e I_e ESC input current A
Ue U_e ESC input voltage V
Ueo U_eo Equivalent ESC output voltage V
σe sigma_e ESC throttle - 0˜1
I b I_b Battery current A

2.1.3 Variable Relation Graph (VRG)

We study the engineering relationship between the various components of the PS, and give the VRG of the most important 14 variables in Fig 1, which can help us visually design function expressions between the variables.

Here is the dependency relationship between all variables. From this picture, you can intuitively see how to calculate between two variables or multiple variables. The double-headed arrow means that two variables can be freely converted. The dashed line represents the derived simple relationship. A one-way arrow refers to a situation where there are multiple inputs and a single output. Although the reverse solution cannot be obtained directly, because there are also relationships between the input variables, the solution process can still be constructed.

This picture can help us deal with the forward and reverse problems in the propulsion system. Forward solution refers to solving the aircraft’s performance indicators (longest hovering time, maximum load, maximum pitch angle, etc.) based on known hardware parameters (propeller, motor, ESC, battery, etc.). The reverse solution is to solve the hardware parameters according to the given aircraft performance index.

PIC

Figure 1: Variable Relationship Graph of Propulsion Subsystem

The graph is divided into four parts, namely the propeller, motor, ESC and battery. Each circle in the graph represents a variable, and the connection represents a function. It can be seen from the graph that there is a one-to-one function obtained by one or more steps between any nodes, although some function forms can be more complex.

2.1.4 Definition of basic functions

Basic functions are given by engineering principles that describe the most basic conversion relationship between variables.

Calculate T by N, see [1][p.66, 4.1]

( N )2 get T-by N (N) = ρCTD4P 60
(2.8)

Calculate M by N, see[1][p.66, 4.2]

( N )2 get M-by N (N) = CM ρ -- D5P 60
(2.9)

Calculate Ea by N, see[1][p.80, 4.61]

get E a-by N (N) = KEN
(2.10)

Calculate M by Im, see[1][p.80, 4.59]

get M-by-I m (I ) = K (I - I ) m T m m0
(2.11)

Calculate Um by Ea and Im, see[1][p.72]

get U-m by-E-a and-I-m(Ea,Im) = Ea + RmIm
(2.12)

Calculate Ueo by Um and Im, see[1][p.69, 4.12]

get U eo-by U m-and I-m(Um, Im) = Um +ImRe
(2.13)

Calculate Ueo by sigmae, see[1][p.69, 4.13]

get U eo-by sigma-e(sigma ) = sigma U e e b
(2.14)

Calculate Ie by sigmae and Im, see[1][p.69, 4.14]

get-I e-by sigma-e and-I-m(sigmae,Im) = σeIm
(2.15)

Calculate Ib by Ie, see[1][p.70, 4.17]

get I-b by I-e(Ie) = nrIe + Iother
(2.16)

Calculate Ue by Ib, see[1][p.70, 4.20]

get U e-by-I-b(I ) = U - I R b b b b
(2.17)

Calculate Tb by Ib, see[1][p.70, 4.22]. Notice that, we simplify the battery model. First, we assume the battery voltage is a fixed value when it is discharging. Second, we assume the remain capacity of battery is linear changed when it is discharging.

get-T-b by I-b(I ) = Cb --Cmin-60 b Ib 1000
(2.18)

Calculate Gmaxload by T, see[1][p.74, 4.35]

get G maxload-by T (T ) = nrT - G
(2.19)

Calculate θmax by T, see[1][p.74, 4.36]

-G-- get theta max-by-T(T) = arccosnrT
(2.20)

Calculate η by M and N and Ib, see[1][p.73, 4.32]

26π0nrM N get eta-by-M-and N-and I-b(M, N, Ib) =--UbIb--
(2.21)

2.1.5 Definitions of basic inverse functions

The basic inverse functions are some inverse functions of unary functions which belongs to basic function.

Calculate N by T, it is an inverse function of (2.8).

∘ ------- T get N-by-T(T) = 60 ρC-D4-- T p
(2.22)

Calculate N by M, it is an inverse function of (2.9).

∘ ------- get N-by-M (M ) = 60 --M---- ρD5pCM
(2.23)

Calculate N by Ea, it is an inverse function of (2.10).

get N-by E a(E ) = Ea a KE
(2.24)

Calculate Im by M, it is an inverse function of (2.11).

get-I m by-M (M ) =-M + I KT m0
(2.25)

Calculate sigmae by Ueo, it is an inverse function of (2.14).

get sigma-e by-U-eo(U ) = Ueo Ub
(2.26)

Calculate Ie by Ib, it is an inverse function of (2.16).

get I e-by I-b(I ) = Ib---Iother e nr
(2.27)

Calculate Ib by Ue, it is an inverse function of (2.17).

get I-b by-U-e(U ) = Ub---Ue e Rb
(2.28)

Calculate Ib by Tb, it is an inverse function of (2.18).

get I-b by-T-b(Tb) = Cb --Cmin-60 Tb 1000
(2.29)

Calculate T by Gmaxload, it is an inverse function of (2.19).

get-T-by-G maxload(G ) = Gmaxload +-G- maxload nr
(2.30)

Calculate T by θmax, it is an inverse function of (2.20).

get T-by theta max (θ ) =----G----- max nrcosθmax
(2.31)

2.1.6 Definitions of extended functions

Extended function are some composited functions by basic functions and basic inverse functions, and the compositions are one or more steps.

Calculate T by M, it is composed by (2.23) and (2.8).

get-T by-M (M ) = CT-*-M- DpCM
(2.32)

Calculate M by T, it is an inversion function of (2.32).

C D T get-M by-T(T) = -M--p-- CT
(2.33)

Calculate Im by T, it is composed by (2.33) and (2.25).

CM DpT get I-m by-T(T) = K-C---+ Im0 T T
(2.34)

Calculate T by M, it is composed by (2.22) and (2.10).

∘ ------- get E a-by T (T ) = 60KE--T---- ρCTD4p
(2.35)

Calculate Um by T, it is composed by (2.35) and (2.34) and (2.12).

∘ ------- ( ) get U m-by-T(T ) = 60KE --T---+ Rm CM-DpT- + Im0 ρCTD4p KT CT
(2.36)

Calculate Um by N and M, it is composed by (2.10) and (2.25) and (2.12).

( ) get-U-m-by N-and M (N, M ) = K N + R M--+ I E m KT m0
(2.37)

Calculate Um by N, it is composed by (2.9) and (2.37).

( ( ) ) CM ρ N60 2D5P get-U-m-by N (N) = KEN + Rm ----KT------+ Im0
(2.38)

Calculate Um by M, it is composed by (2.23) and (2.37).

∘ ------- ( ) get U m-by M (M ) = 60KE --M----+ Rm -M- +Im0 ρD5pCM KT
(2.39)

Calculate Im by Ea, it is composed by (2.12) and (2.13).

get-U-eo by-E-a and-I-m(Ea,Im) = Ea + (Rm + Re )Im
(2.40)

Calculate Ueo by N, it is composed by (2.13), (2.12), (2.10), (2.25) and (2.9).

5 get U eo-by-N (N ) = (Rm-+-Re)CM-ρDpN 2 + KEN + (Rm + Re)Im0 KT602
(2.41)

Calculate sigmae by Ea and Im, it is composed by (2.26) and (2.40).

E + (R + R )I get sigma-e by-E-a and I-m(Ea,Im ) =-a--m----e--m Ub
(2.42)

Calculate Ie by Ea and Im, it is composed by (2.42) and (2.15).

(E +(R +R )I )I get I-e by-E-a and I-m(Ea,Im ) =-a---m----e-m--m- Ub
(2.43)

Calculate Ie by T, it is composed by (2.43), (2.35) and (2.34).

∘ --T--- CMDpT- CMDpT- get-I e-by T(T ) = ((60KE-ρCTD4p)+-(Rm-+-Re)(KT-CT-+-Im0))(-KTCT--+-Im0-) Ub
(2.44)

Calculate Ue by Ie, it is composed by (2.17) and (2.16).

get U-e by-I e(Ie) = Ub - (nrIe + Iother)Rb
(2.45)

Calculate Tb by Ie, it is composed by (2.18) and (2.16).

-Cb --Cmin--60- get-T-b by I-e(Ie) = nrIe + Iother 1000
(2.46)

2.1.7 Definitions of advanced extended functions

Advanced extended function are some extended functions that would got several solutions. Usually, these functions are complex higher order equations. In actual situations, only one solution will be left based on the context constraint.

In fact, these definitions are too difficult so that we don’t like it. We will deform and simplify these formulas appropriately before using it.

Calculate M by Um, it is an inverse function of (2.38).

( ∘ ----------------------------) 1800KT RmCM ρD5p(RmIm0 - Um ) get-N-by U m (Um ) = R-C-ρD5-( - KE ± K2E - --------900K----------) m M p T
(2.47)

Calculate N by Ueo, it is an inverse function of (2.48).

( ∘ ----------------------------------------) ----1800KT------( 2 (Rm-+-Re-)CM-ρD5p ) get N by-U-eo(Ueo) = (Rm + Re)CM ρD5p - KE ± K E - 900KT ((Rm + Re)Im0 - Ueo)
(2.48)

2.1.8 Simplest Form of Functions (SFFs)

SFF is the simplest function in mathematical form relative to the variables wer care about. This is done by rearranging the order of operations, separating the arguments and making the whole coffecient expression as a constant. Then the calculation cost will be reduced to minimize. It is very useful to improve calculation performance, especially when the function will be called frequently.

Another benefit to us is that we can see the type of the new function more clearly contrast to the old form.

For example, the new SFF for calculating motor input voltage Um by propeller speed N is (2.49), it is an equivelent transformation of (2.38).

RmCM--ρD5p- 2 get U m-by N-sff(N) = 3600KT N + KEN + RmIm0
(2.49)

We can analyze the difference between (2.38) and (2.49) in calculation.

PIC

Figure 2: Comparison of traditional functions and SFF

As can be seen from Fig 2, the new SFF are calculated less and the number of intermediate variables are reduced too. At the same time, the type of the new SFF is more clearly than old form, it is an unary quadratic function form.

This comparison shows that SFF is a better approach and may greatly improve the efficiency of engineering computing.

We have done equivalent transformations to many original functions, and given the SFF table, see Fig 3.

PIC

Figure 3: Simplest Form of Function Table of Propulsion Subsystem

The first column of the table is the input, and the first row is the output, the a, b, c in the table represent a constant expression in SFF, and the a, b, c in different functions is not the same value, just means they are constant expressions. We haven’t filled all blank grids because some functions only have mathematical meaning, haven’t engineering meaning, and the form is more complicated, so it is not the first priority to complete. Because we know from VRG that there are calculation paths between any nodes, there should be a function for each grid, and there must be an SFF. It should be pointed out that any pair of functions that are symmetrical on the main diagonal are mutually inverse.

There are now 14 important variables in total. After a rough estimate, the number of formulas that can be obtained is 13 + 12 + ... + 1 = 91, but in fact some of them are meaningless and there is no need to list them.

2.2 Formal verification

2.2.1 Customized automatic tactic

In the field of engineering control, the most common is the lengthy mathematical formula of the real number type. The real number type in the Coq system is constructed with axioms, and the related proofs require a lot of lemmas to rewrite and simplify. This is a manual proof example of correctness of a function which composited by two functions about thrust T and propeller speed N.

Lemma verify_trans_N_and_T N : 0 < N -> get_N_by_T (get_T_by_N N) = N.  
Proof.  
  intros H1. unfold get_N_by_T, get_T_by_N.  
  remember (rho * C_T * D_p ^ 4) as a.  
  unfold Rdiv. rewrite Rinv_r_simpl_m.  
  - rewrite sqrt_pow2.  
    + rewrite -> Rmult_assoc. rewrite Rinv_r_simpl_m; trivial.  
    + apply Rle_mult_inv_pos. trivial. lra.  
  - rewrite Heqa.  
    apply Rmult_integral_contrapositive. split.  
    + apply Rmult_integral_contrapositive. split.  
      * apply Rgt_not_eq. apply Rlt_gt; apply gt0_rho.  
      * apply Rgt_not_eq. apply Rlt_gt; apply gt0_C_T.  
    + apply pow_nonzero; apply Rgt_not_eq;  
      apply Rlt_gt;apply gt0_D_p.  
Qed.

It can be seen that this script is cumbersome and poor readable. We found that many proofs are same pattern, first expand the definition, then use some lemmas to apply, try to use the Ring or Field strategy to auto prove the equality, and try to use the Fourier strategy to auto prove the general inequality.

We packaged these automatic proof technologies and provide a real number proof module named R_prove. This module consists of three parts: a table of commonly used lemmas, see Table 2.2.1; custom lemmas; custom strategies.

Lemma name Lemma type
Rlt_0_1 0 < 1
Rmult_le_pos 0 ≤ r1 → 0 ≤ r2 → 0 ≤ r1
Rmult_lt_0_compat 0 < r1 → 0 < r2 → 0 < r1 * r2
Rgt_not_eq r1 > r2 → r1≠r2
Rlt_le r1 < r2 → r1 ≤ r2
pow_lt 0 < x → 0 < xn
pow_nonzero x≠0 → xn≠0
Rsqr_sqrt 0 ≤ x → (sqrtx)2 = x
sqrt_pow2 0 ≤ x → sqrt(x2) = x

We use the proof strategy library management command Hint to add strategies to a specified strategy library. For example, “Hint Resolve Rlt_le : flyctrl.” adds the strategy “apply Rlt_le” to the flyctrl strategy library. Then you can enter “auto with flyctrl” to use this automatic proof strategy library.

Among the custom strategies, one of the most commonly used strategies is “simple_equation”, which can solve most of the proof problems in PS.

Ltac simple_equation := intros;  
  repeat autounfold with flyctrl; auto with flyctrl; unfold Rdiv;  
  autorewrite with flyctrl;try field;try fourier; auto with flyctrl;  
  try apply Rmult_le_pos; auto with flyctrl;  
  try apply Rmult_integral_contrapositive; auto with flyctrl;  
  repeat try split; auto with flyctrl.

The example of manual proof just demonstrated can be done in one step with this strategy.

Lemma verify_trans_N_and_T’ N :  
  0 < N -> get_N_by_T (get_T_by_N N) = N.  
  Proof. simple_equation. Qed.

2.2.2 Verification of equality of function transformation

It is mainly to verify the definitions of basic inverse functions and derivative functions.

The first is the transformation of the paired basic function and the inverse function, named verify_trans_XX_and_YY, where XX and YY represent two variables respectively.

Lemma verify_trans_N_and_T N : get_N_by_T (get_T_by_N N) = N.  
Lemma verify_trans_N_and_M N : get_N_by_M (get_M_by_N N) = N.  
Lemma verify_trans_N_and_E_a N : get_N_by_E_a (get_E_a_by_N N) = N.  
Lemma verify_trans_I_m_and_M I_m : get_I_m_by_M (get_M_by_I_m I_m) = I_m.  
Lemma verify_trans_sigma_e_and_U_eo sigma_e:  
    get_sigma_e_by_U_eo (get_U_eo_by_sigma_e sigma_e) = sigma_e.  
Lemma verify_trans_I_e_and_I_b I_e : get_I_e_by_I_b (get_I_b_by_I_e I_e) = I_e.  
Lemma verify_trans_I_b_and_T_b I_b : get_I_b_by_T_b (get_T_b_by_I_b I_b) = I_b.  
Lemma verify_trans_I_b_and_U_e I_b : get_I_b_by_U_e (get_U_e_by_I_b I_b) = I_b.  
Lemma verify_trans_T_and_G_maxload T :  
    get_T_by_G_maxload (get_G_maxload_by_T T) = T.  
Lemma verify_trans_T_and_sigma_max T :  
    get_T_by_theta_max (get_theta_max_by_T T) = T.

Then let’s verify the definition of the derived function, and named it verify_XXX, where XXX represents the name of the derived function.

Lemma verify_get_T_by_M M : get_T_by_M M = get_T_by_N (get_N_by_M M).  
Lemma verify_get_M_by_T M : get_M_by_T (get_T_by_M M) = M.  
Lemma verify_get_I_m_by_T T : get_I_m_by_T T = get_I_m_by_M (get_M_by_T T).  
Lemma verify_get_E_a_by_T T : get_E_a_by_T T = get_E_a_by_N (get_N_by_T T).  
Lemma verify_get_U_m_by_T T : get_U_m_by_T T =  
    get_U_m_by_E_a_and_I_m (get_E_a_by_T T) (get_I_m_by_T T).  
 
Lemma verify_get_U_m_by_N_and_M N M : get_U_m_by_N_and_M N M =  
    get_U_m_by_E_a_and_I_m (get_E_a_by_N N) (get_I_m_by_M M).  
Lemma verify_get_U_m_by_N N : get_U_m_by_N N =  
    get_U_m_by_N_and_M N (get_M_by_N N).  
Lemma verify_get_U_m_by_M M : get_U_m_by_M M =  
    get_U_m_by_N_and_M (get_N_by_M M) M.  
 
Lemma verify_get_U_eo_by_E_a_and_I_m E_a I_m :  
    get_U_eo_by_E_a_and_I_m E_a I_m =  
    get_U_eo_by_U_m_and_I_m (get_U_m_by_E_a_and_I_m E_a I_m) I_m.  
Lemma verify_get_U_eo_by_N N : get_U_eo_by_N N =  
    get_U_eo_by_E_a_and_I_m (get_E_a_by_N N)  
    (get_I_m_by_M (get_M_by_N N)).  
 
Lemma verify_get_sigma_e_by_E_a_and_I_m E_a I_m :  
    get_sigma_e_by_E_a_and_I_m E_a I_m =  
    get_sigma_e_by_U_eo (get_U_eo_by_E_a_and_I_m E_a I_m).  
Lemma verify_get_I_e_by_E_a_and_I_m E_a I_m : get_I_e_by_E_a_and_I_m E_a I_m =  
    get_I_e_by_sigma_e_and_I_m (get_sigma_e_by_E_a_and_I_m E_a I_m) I_m.  
Lemma verify_get_I_e_by_T T : get_I_e_by_T T =  
    get_I_e_by_E_a_and_I_m (get_E_a_by_T T) (get_I_m_by_T T).  
Lemma verify_get_U_e_by_I_e I_e : get_U_e_by_I_e I_e =  
    get_U_e_by_I_b (get_I_b_by_I_e I_e).  
Lemma verify_get_T_b_by_I_e I_e : get_T_b_by_I_e I_e =  
    get_T_b_by_I_b (get_I_b_by_I_e I_e).  
Lemma verify_get_N_by_U_m U_m : get_U_m_by_N (get_N_by_U_m U_m) = U_m.  
Lemma verify_get_N_by_U_eo U_eo : get_U_eo_by_N (get_N_by_U_eo U_eo) = U_eo.  

2.3 Example 1: Estimate the longest hovering time

2.3.1 Analysis of the question

Question: In the hover mode, according to the known parameters, estimate the time of endurance Tb, see [1] P75.

Principle: In the hovering state, the thrust generated by the PS is exactly equal to the total weight of the aircraft.

2.3.2 General calculation step

Calculation steps, see [1] p81-p82:

(1).
Calculate the thrust required by a single propeller, simply divide the total mass by the number of propellers.
T = G∕n-r

(2).
Calculate the propeller speed
N = get N-by T (T )

(3).
Calculate the propeller torque
M = get M by-T(T)

(4).
Calculate the back-electromotive force of the motor
Ea = get-E-a by-N (N )

(5).
Calculate the input current of the motor
Im = get I-m by-M (M )

(6).
Calculate the input voltage of the motor
Um = get U m by-E a-and-I m (Ea,Im)

(7).
Calculate the output voltage of the ESC
Ueo = get U-eo by-U-m and-I-m(Um, Im )

(8).
Calculate the ESC throttle
σe = get-sigma-e by-U-eo and U b(Ueo,Ub)

(9).
Calculate the input current of the ESC
Ie = get I e-by I-m and-sigma-e(Im, sigmae )

(10).
Calculate the output current of the battery
Ib = get I-b by-I-e(Ie)

(11).
Calculate the hover time
Tb = get T b-by I-b(Ib)

2.3.3 Direct calculate formulas

The formula of direct calculation is given here, which provides another way of calculation.

A direct formula, calculate N.

∘ --------- direct N = 60 ----G---- ρD4pCTnr
(2.50)

A direct formula, calculate M.

direct M = CM-DpG-- CT nr
(2.51)

A direct formula, calculate Ea.

∘ --------- G direct E-a = 60KE ρD4C--n-- p T r
(2.52)

A direct formula, calculate Im.

CM DpG direct I m = C-n-K--+ Im0 T r T
(2.53)

A direct formula, calculate Um.

∘ --------- ( ) direct U-m = 60KE ---G4-----+ Rm CM-DpG--+Im0 ρD pCTnr CTKT nr
(2.54)

A direct formula, calculate Ueo.

∘ --------- G ( CM DpG ) direct U-eo = 60KE ρD4-C-n--+ (Rm + Re) C-K--n--+ Im0 p T r T T r
(2.55)

A direct formula, calculate σe.

∘ ---G---- (CMDpG-- ) 60KE---ρD4pCTnr +-(Rm-+-Re-)-CTKTnr-+-Im0- direct-sigma e = Ub
(2.56)

A direct formula, calculate Ie.

∘ -------- ( ) 60KE ρD4GC-n-+ (Rm + Re) CCMDKpGn--+Im0 ( ) direct-I e =---------p-T-r--------------T-T-r------- CM-DpG--+ Im0 Ub CTnrKT
(2.57)

A direct formula, calculate Ib.

( ∘ ---G---- (-CMDpG- ) ( ) ) |60KE---ρD4pCTnr-+-(Rm--+-Re)-CTKT-nr +-Im0 -CM-DpG- | direct I-b = nr ( Ub CT nrKT + Im0 ) + Iother
(2.58)

A direct formula, calculate Ue.

( ( ∘ -------- ( ) ) ) | | 60KE ρD4pGCTnr +(Rm + Re ) CCMTDKpTGnr-+ Im0 (CM DpG ) | | direct U e = Ub-(nr (-------------------U-------------------- C--n-K--+ Im0 ) + Iother) Rb b T r T
(2.59)

A direct formula, calculate Tb.

Cb - Cmin 60 direct-T-b =---(60K-∘----G--+-(R--+R-)(-CMDpG+I--)---------------)--------1000 n (---E--ρD4pCTnr---m---e-CTKT-nr--m0- (CMDpG--+ I )) + I r Ub CTnrKT m0 other
(2.60)

Note that these formulas are not so handsome but really correct. If we transform them using the SFF method, then we will get a more elegant form.

2.3.4 Vefiry direct formulas

Verify that the results calculated by general method are equal to the new given direct fomulas, so we can ensure that the new formulas are reliable.

Lemma direct_N_correct : direct_N = N.  
Lemma direct_M_correct : direct_M = M.  
Lemma direct_E_a_correct : direct_E_a = E_a.  
Lemma direct_I_m_correct : direct_I_m = I_m.  
Lemma direct_U_m_correct : direct_U_m = U_m.  
Lemma direct_U_eo_correct : direct_U_eo = U_eo.  
Lemma direct_sigma_e_correct : direct_sigma_e = sigma_e.  
Lemma direct_I_e_correct : direct_I_e = I_e.  
Lemma direct_I_b_correct : direct_I_b = I_b.  
Lemma direct_U_e_correct : direct_U_e = U_e.  
Lemma direct_T_b_correct : direct_T_b = T_b.

2.4 Example 2: Estimate system effencicy at max throttle

2.4.1 Analysis of the question

Question: Solve the motor voltage, motor current, propeller speed and system efficiency under the maximum throttle mode.

Principle: Under the maximum throttle, the electric throttle is equal to 1. At this time, multiple equations are combined to solve the result.

For the solution of the key variable ”propeller speed N”, literature [1][p71, 4.30] gives a system of equations, A numerical iteration method is proposed to solve the problem, but the final expression is not given. The author gives the formula of direct calculation.

2.4.2 General calculation step

(1).
Given the throttle, see [1][p.71, 4.30]
sigmae = 1
(2.61)

(2).
Calculate Ueo
Ueo = get U-eo by-sigma-e(sigmae)
(2.62)

(3).
Calculate N
N = get N by-Ueo(Ueo)
(2.63)

(4).
Calculate M
M = get M-by N (N)
(2.64)

(5).
Calculate Im
Im = get I-m by M (M )
(2.65)

(6).
Calculate Ie
Ie = get-I e-by I-m-and-sigma e(Im,sigmae)
(2.66)

(7).
Calculate Ie
Ib = get I-b by-I e(Ie)
(2.67)

(8).
Calculate η
η = get eta by-M and N and I-b(M, N, Ib)
(2.68)

2.5 Example 3: Estimate maximum load and maximum pitch Angle

2.5.1 Analysis of the question

Question: Calculate the maximum load weight and the maximum pitch Angle of multiple rotors with a certain safety margin (e.g., 20%) reserved for throttle commands.

Principle: Use the maximum throttle in safe mode, here 0.8, then combine multiple equations to solve the result.

For the solution of the key variable ”propeller speed N”, literature [1][p73, 4.33] gives a system of equations, A numerical iteration method is proposed to solve the problem, but the final expression is not given. This problem is similar to the ”maximum throttle mode” problem, and the direct calculation formula is given.

2.5.2 General calculation step

(1).
Given the throttle, see [1][p.73, 4.33]
sigmae = 0.8
(2.69)

(2).
Calculate Ueo
Ueo = get U-eo by-sigma-e(sigmae)
(2.70)

(3).
Calculate N
N = get N by-Ueo(Ueo)
(2.71)

(4).
Calculate T
T = get T-by N (N )
(2.72)

(5).
Calculate Gmaxload
Gmaxload = get G-maxload by-T(T)
(2.73)

(6).
Calculate θmax
θmax = get θ max-by-T(T )
(2.74)

2.6 Example 4: Estimate maximum flat flight speed and distance

2.6.1 Analysis of the question

Question: Designers are concerned about two performance metrics: maximum flat flight speed and distance. The function relationship between pitch Angle and flight speed can be obtained after the model is established by force analysis. The solving method given by the author of the original book is numerical iteration method, but the author wants to try to find a better formula for direct calculation. The knowledge of derivative and extreme point is introduced to try to solve the extreme point. There is no conclusion and the issue is on hold.

Principle: The function relationship between flight speed and pitch Angle is established firstly, and then the pitch Angle when the maximum flight speed is obtained is solved. Then other parameters are easy to be solved. Similarly, the maximum flight distance can also be established as a function of pitch Angle, which can be solved in the same way.

3 Summary

In order ot realize control system as computer software, it is necessary to extract the exisiting engineering knowledge into the reference manual of software programming model required by computer programming. This process is traditionally completed by people manully, and its reliability cannot be guaranteed. We use the theorem prover COQ to gradually establish a reliable knowledge base for the control system from the bottom level in the way of mathematical reasoning. At the same time, we also get a reliable model that has been formally verified. We take the flight control system of small unmanned aerial vehicle as the research object, and analyze the control principle and algorithm of each subsystem. The research object of this chapter is propulsion system modeling.

In future work, we will use program generation technology to automatically generate high reliability control system software from this model.

References

[1]   Quan, Q. . Introduction to Multicopter Design and Control. Springer Singapore, 2017.

[2]   Shi, D. , et al. ”A Practical Performance Evaluation Method for Electric Multicopters.” IEEE/ASME Transactions on Mechatronics 22.3(2017):1337-1348.