Piecewise Polynomial Modeling for Control and Analysis of Aircraft Dynamics Beyond Stall

TLDR: HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not, for teaching and research institutions in France or abroad, or from public or private research centers.

Abstract: HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Piecewise Polynomial Modeling for Control and Analysis of Aircraft Dynamics beyond Stall Torbjørn Cunis, Laurent Burlion, Jean-Philippe Condomines

Figures

Fig. 2 Axis systems with angles and vectors (projections into the plane are marked by ′).

Fig. 5 Continuation of the longitudinal dynamics of the GTM for polynomial and piecewise defined models of the aerodynamic coefficients.

Fig. 4 (Part 2) Piecewise model of the aerodynamic coefficients in angle of attack and side-slip angle for neutral surface deflections.

Fig. 4 (Part 2) Piecewise model of the aerodynamic coefficients in angle of attack and side-slip angle for neutral surface deflections.

Fig. 3 Comparison of 3rd-order polynomial [10, 11] and piecewise identifications. [12]

Fig. 1 Flow chart for the piecewise polynomial fitting approach: a) without continuity constraint; b) with continuity constraint in x0.

Full text

Piecewise Polynomial Modeling for Control and
Analysis of Aircraft Dynamics beyond Stall
Torbjørn Cunis
*
and Laurent Burlion
†
ONERA – The French Aerospace Lab, Centre Midi-Pyrénées, Toulouse, 31055, France
Jean-Philippe Condomines
‡
French Civil Aviation School, Toulouse, 31055, France
Nomenclature
α
0
= Low-angle of attack boundary (°);
ρ = Pseudo-radius (ρ ∈ R);
ϕ
(
·
)
= Boundary condition function (ϕ : R
m
→ R);
Σ = Positive-denite shape factor (Σ ∈ R
n×n
);
Ω
ϕ
= Boundary curve set (Ω
φ
⊂ R
m
);
C
l,m,n
= Aerodynamic coecients of moments in body axes (·);
C
X,Y,Z
= Aerodynamic coecients of forces in body axes (·);
C = Objective matrix (C ∈ R
k×r
);
d = Vector of measurements (d ∈ R
k
);
f
(
·
)
= Non-linear, open-loop system dynamics (f :
(
X, ·
)
7→
˙
X);
g
(
·
)
, h
(
·
)
= Positive-semi-denite Lagrange multiplier (g, h : R
n
→ R
≥0
);
k = Number of measurements;
m = Number of variables;
n = Number of states; system degree; polynomial degree;
q = Vector of coecients (q ∈ R
r
);
r = Number of coecients;
*
Doctoral Researcher, Department of Information Processing and Systems, e
-
mail:
torbjoern.cunis@onera.fr
; associated
researcher with the French Civil Aviation School, Drones Research Group; AIAA Student Member.
†
Research Scientist, Department of Information Processing and Systems, e-mail: laurent.burlion@onera.fr.
‡
Assistant Professor, Drones Research Group, e-mail: jean-philippe.condomines@enac.fr.

X
∗
, µ
∗
= States and parameters at trim condition;
(
·
)
post
= Domain of high angle of attack;
(
·
)
pre
= Domain of low angle of attack;
X = State space (X ⊆ R
n
);
∂X = Set boundary of X.
I. Introduction
F
ull-envelope
aircraft models require extensive eort to represent the aerodynamic coecients well in
the entire region of the envelope as ight dynamics beyond stall are highly non-linear and often unstable
[
1
,
2
]. With upset recovery approaches found in the literature being model-based ([
3
–
5
], and references
herein) there is a clear need for reliable full-envelope models of ight dynamics. NASA’s Generic Transport
Model (GTM) has contributed signicantly to analysis and control approaches of civil and unmanned aircraft
over the entire ight envelope (see, e.g., [
6
,
7
]). Representing a
5.5 %
down-scaled typical aerial transport
vehicle, the GTM provides exhaustive, full-envelope aerodynamic data from wind-tunnel studies [
8
] and its
open-source aerodynamic model for MATLAB/Simulink [
9
] has given access for development of modeling,
analysis, and control methods to the aerospace community. An overview of research studies on longitudinal
trim conditions, regions of attraction, and upset situations can be found in [
2
,
5
,
10
,
11
]. However, analytical
representations proposed for the full-envelope aerodynamics are still insucient for non-linear analysis and
control design [
12
]. Therefore, improved methods for accurate modeling are imperative, in particular when
developing robust and powerful advanced control strategies for upset recovery. Subsequently, model feedback
designs based on full-envelope aerodynamic models will grant full authority and control eciency for stability
and performances in unmanned aircraft (UA) [13].
Polynomial models of the aerodynamic coecients have provided a constructive method to dene and
evaluate models based on analytical computation due to their continuous and dierentiable nature. Despite
the fact that polynomial models have been published recently [
10
,
11
], none of the results represent the
aerodynamic coecients well in the entire region of the envelope [
12
]. Indeed, at the stall angle of attack, the
laminar ow around the wings in the pre-stall region changes to turbulent ow and remains so in post-stall.
This signicant change of the ow dynamics motivates a piecewise model of the pre-stall and post-stall
dynamics instead.
Piecewise regression theory can be dated to the
1970
s; rst research focused on regression of a few
polynomial functions piecewise over the observations. However, the estimation of suitable switching surfaces
2

or joints for the piecewise functions usually adds computational diculty and load [
14
–
17
]. Later, multivariate
splines were introduced [
18
,
19
]; using simplices and baryocentric coordinates for the base, the so-called
B-splines bear the advantage of generalized continuity, dimensional exibility, and ecient evaluation as
well as a stable local basis [
20
]. Recently, a further approach combining splines with fuzzy logic has been
presented in [21].
While splines today present a powerful yet complex tool for accurate and smooth interpolation, they lack
an underlying physical model justifying the partition. Moreover, for functional analysis of trim conditions
and stable sets, as in [
10
,
11
], splines weren’t used but polynomials. Motivated by the practical problems
encountered with mini-UAs ight control and guidance, civil aircraft fault detection and isolation, and upset
recovery, we aim to derive a simple yet powerful aerodynamic model still suitable for functional analysis. A
novel approach for piecewise polynomial modeling aerodynamic coecients, the
pwpfit
toolbox for MATLAB,
was recently proposed in [
22
]. Here, we have proven feasibility of tting both a piecewise polynomial model
and its joint surface using linear least-square (LSQ) optimization techniques. While this approach is limited
to a single joint without dierential continuity, the switch in the dynamics is motivated by the change from
laminar to turbulent ow at stall and the resulting model is found to t the full-envelope aerodynamics well.
This article focuses on the recent research detailing the theoretical aspects in the sequel and their
application to functional analysis. The main contributions of this paper are therefore: to address (in
§II) a concise bibliographical review of the polynomial based-methods used for full-envelope identication;
to introduce (in §III) a novel and generic formulation of the piecewise polynomial tting method which
approximates a piecewise polynomial function and its joint; to provide (in §IV) a six-degrees-of-freedom model
of an aircraft and its aerodynamic coecients, accounting for both pre-stall and post-stall characteristics by
piecewise identication; and nally to demonstrate and assess the extension of functional analysis tools for
the piecewise polynomial model (in §V).
II. State of the Art
A. Polynomial regression
Polynomial regression is a general approach similar to linear curve regression, where a polynomial function
f
is to be found in order to approximate best a set of measured data points. Here, the coecients of
f
are
subject to the optimality problem of minimal sum of squared residuals of
f
with respect to the measurements
(goodness of t, GoF). The formulation of optimal coecients as a linear least-square problem dates back
to Legendre (1805) and Gauss (1809); a rst application can be found by Gergonne in 1815 [
23
]. It has
been shown that on average, the residuals of such an optimal polynomial vanish and their deviation is
3

minimized [
24
]. Furthermore, polynomial functions are dened by basic mathematical operations of addition
and multiplication and thus provide by their nature smoothness to innite dierentiation. Recent polynomial
models of the full-envelope aerodynamic coecients of the GTM have been presented in [10, 11].
B. Multi-variate splines
Splines are piecewise sequences of polynomial functions, where each polynomial is active only in the
respective partition. These partitions are chosen before tting instead of being subject of the t. The
polynomials sub-functions are computed such that at the boundaries of the selected partitions the overall
spline function is smooth to a certain degree of continuity. Thus, spline functions show characteristics of both
lookup tables and polynomials, as noted by de Visser et al. [18, p. 3]:
Eectively, spline functions [...] combine the global nonlinear modeling capability of lookup tables
with the analytic, continuous nature of polynomials.
While for a single-variable spline function, the boundaries equal point-wise joints, the partitions of
multi-variate splines can be more complex. In addition to simple rectangles (or rather rectangular polytopes),
triangular partitions have recently proposed by [
20
]. However, the high accuracy of splines in terms of their
residuals is opposed by their computational costs for further analytical investigation. Multi-variate splines
have been used in, among others, [18, 19] in order to model full-envelope aerodynamics.
III. Methodology
In vector notation, optimal coecients for piecewise polynomial ts are expressed as linear least-squares.
We introduce a polynomial notation by the vectors of monomials and coecients and thus reduce the goodness
of t to a function of the latter. The joint will be given by the scalar eld
ϕ
(
·
)
in the variables of the model
and the scalar bound
x
0
; here, we assume
ϕ
to be linear matrix inequality and
x
0
will be determined by the
tting. The resulting model then has a single joint with value continuity, i.e., the model is not dierentiable
at the joint. For models in several variables and outputs, such as the aerodynamic coecients, it is desirable
to have further constraints to the t enforced. We will add those desired properties of the piecewise t as
constraint matrices.
Denition 1 A linear least-square (LSQ) problem is given as the optimization problem
lsq
(
C, d, A, 0
)
= arg min
q∈Ω
A
k
Cq − d
k
2
2
. (1)
with q ∈ R
r
, C ∈ R
k×r
, d ∈ R
k
, and Ω
A
=
{
q
|
Aq = 0
}
for a constraint matrix A.
4

A. Polynomials
A monomial of degree
n
is a single product of powers where the exponents add up to the total degree
n
,
without any scalar coecient. We notate a monomial of x =
(
x
1
, . . . , x
m
)
in degrees n =
(
n
1
, . . . , n
m
)
as
x
n
= x
n
1
1
. . . x
n
m
m
, (2)
where x
n
has the total degree n =
k
n
k
1
= n
1
+ · · · + n
m
.
Denition 2 P
n
(
x
)
is the vector of monomials
x
ν
in variables
x =
(
x
1
, . . . , x
m
)
with degrees
ν ∈ N
m
and
total degrees
k
ν
k
1
≤ n; and the number of elements in P
n
(
x
)
is denoted by r
[
n
]
, i.e., P
n
∈ R
[
x
]
r
[
n
]
.
By this notation, a polynomial f is expressed as scalar product of its monomials and coecients,
f
(
x
)
= hP
n
(
x
)
, qi (3)
with the vector of coecients q ∈ R
r
[
n
]
.
B. Piecewise polynomial tting
Consider the k observations
(
x
i
, z
i
)
i
given as sequences over i ∈
[
1, k
]
:
z
i
= γ
(
x
i
)
+ 
i
, (4)
where
(
x
i
, z
i
, 
i
)
1≤i ≤k
⊂ R
m
× R × R
and
γ
(
·
)
and
(

i
)
i
are an unknown function and measurement error,
respectively; we will nd coecients q
1
, q
2
as well as a scalar x
0
∈ R such that
f : x 7−→











hP
n
(
x
)
, q
1
i if ϕ
(
x
)
≤ x
0
;
hP
n
(
x
)
, q
2
i else;
with ϕ: R
m
→ R minimizes the sum of squared residuals
GoF
(
f
)
=
def
k
X
i=1


f
(
x
i
)
− z
i


2
(5)
for an
n > 0
. We note the sub-polynomials of
f
by
f
1,2
: X
1,2
→ R, x 7→ hP
n
(
x
)
, q
1,2
i
with
X
1
∩ X
2
= ∅
and
call
X
1
∪ X
2
the entire domain of
f
. The joint of
f
is given by
Ω
ϕ
=
def
∂X
1
∩ ∂X
2
=

x


ϕ
(
x
)
= x
0

; if
ϕ
(
·
)
is
a linear matrix inequality, the boundary is convex. Re-writing the goodness of t using matrix calculus, we
reduce the cost functional to a cost function and polynomial data tting to a linear least-square problem.
5

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "Piecewise polynomial modeling for control and analysis of aircraft dynamics beyond stall" ?

In this paper, a piecewise polynomial fitting of the aerodynamic coefficients with a single non-smooth joint representing the change of dynamics at high angles of attack is proposed. 

Q2. What is the way to reduce the goodness of fit?

Re-writing the goodness of fit using matrix calculus, the authors reduce the cost functional to a cost function and polynomial data fitting to a linear least-square problem. 

Q3. What is the importance of knowing the flight envelope?

When looking for the full-envelope, non-linear behaviour of an aircraft, knowledge about the flight envelope is vitally important for the vehicle’s safety. 

Q4. How many iterations of the polynomial short-period model?

After 94 iterations, the authors obtain the stable setXρ0 = { X P(X) ≤ ρ0 } withρ0 = 1.4404 (29)and the dissipative region { X V0(X) ≤ λ0 } ⊆ {X V̇0(X) < 0 } of the quartic Lyapunov function V0(·) with λ0 = 0.3265. 

Q5. What is the simplest way to prove that a set of aerial vehicles is stable?

In order to compute and prove minimal stable sets of aerial vehicles based on Lyapunov function theory, researchers have successfully applied sum-of-squares programming for smooth polynomial models and ellipsoid-shaped sets [11, 29]. 

Q6. what is the simplest equation for the LSQ problem?

aTx ≤ x0 be a linear matrix inequality (LMI) with aT = [ a1 · · · am ] and a1 , 0; a piecewise polynomial function f with continuity in Ωϕ is subject to the constrained LSQ problem given by the continuity constraint matrix C , i.e.,C q1 q2 = 0⇐⇒ ∀x ∈ Ωϕ . 〈Pn(x) , q1〉 = 〈Pn(x) , q2〉. (10)The constraint matrix C is constructed by separation of the constrained variables into Λ0 such that〈Pn(x0) , q1,2〉 = 〈Pn(x̃) ,Λ0q1,2〉 (11)for all x0 ∈ { x a Tx = x0 } , where x̃ are the remaining free variables; the authors then have that〈Pn(x̃) ,Λ0q1〉 = 〈Pn(x̃) ,Λ0q2〉for all x̃ ∈ Rm−1 if and only if Λ0q1 = Λ0q2 and hence, Equation (10) holds for C = [ Λ0 −Λ0 ] . 

Q7. What are the tools that can be used to compute the continuation and bifurcation of continuous?

Toolboxes like the Continuation Core and Toolboxes (COCO) [26] offer computation of continuation and bifurcation of continuous functions. 

Q8. What is the way to model a polynomial fitting?

Due to measurement errors or modelling flaws, a polynomial fitting may have relations that either do not exist or shall not be modeled; e.g., for a symmetric aircraft aligned to the flow, there is no side-force— regardless its angle of attack. 

Q9. What is the implementation of the LSQ problem?

〈Pn ( x′, 0m−j ) , q〉 = 0 (12)and 0m−j ∈ {0}m−j is subject to the constrained LSQ problem given by the zero constraint matrix Z.D. Implementation 

Q10. What is the aerodynamic coefficients of the GTM?

the changes of air speed V , side-slip β, inclination γ, and azimuth χ are subject to lift, drag, thrust, and side-force (given by the aerodynamic force coefficients CX,CY,CZ and the thrust input F); the changes of angular body rates ṗ, q̇, ṙ are given by the aerodynamic moment coefficients Cl,Cm,Cn; and the changes of attitude Φ,Θ,Ψ are functions of the angular body rates. 

Q11. What is the simplest way to compute a stable set of a piecewise defined system?

the authors compute a stable set of the piecewise model for the short-period motion, α̇ q̇ = fpresp (α, q, η = η∗) if α ≤ α0; fpostsp (α, q, η = η∗) else; (28)in the neighbourhood of the trim condition V ∗ = 45.7 m/s, γ∗ = 0, α∗ = 3.75°, η∗ = 1.49°, and F∗ = 21.44 N with the shape factorΣ = diag(20°, 50 °/s)−2accounting for the physical operation range of the Generic Transport Model at the selected trim condition. 

Q12. What is the corresponding figure for the piecewise model?

The piecewise model shows additionally a section of stable trim conditions along the branch of increasing angles of attack, corresponding to thrust inputs larger than 135 N—i.e., 100 % throttle—and elevator deflections of −19° to −23°. 

Q13. What is the way to find the coefficients of f?

the coefficients of f are subject to the optimality problem of minimal sum of squared residuals of f with respect to the measurements (goodness of fit, GoF). 

Q14. What is the linear matrix inequality of aTx x0?

If the linear matrix inequality aTx ≤ x0 of Proposition 1 is given with a single non-zero element of a – suppose, a1 , 0 –, separation of the assigned variable x1 ≡ x0 yieldspN (x0, x̃) = λN (x0) T pN (x̃)with x̃ ∈ Rm−1 and λN = diag pN (x0, 1m−1) for 1m−1 ∈ {1}m−1. 

Q15. What is the way to reduce the goodness of fit to a function of the latter?

The authors introduce a polynomial notation by the vectors of monomials and coefficients and thus reduce the goodness of fit to a function of the latter. 

Q16. What is the way to evaluate the cost functional for f?

The cost functional for f can be evaluated piecewise toGoF( f ) = ∑xi ∈X1f1(xi) − zi 2+ ∑xi ∈X2f2(xi) − zi 2 (6)with X1 = {x1, . . . , xi′ }, X2 = { xi′+1, . . . , xk } chosen a priori. 

Q17. What is the definition of the continuity constraint matrix C?

With the definition Pn = [ p0(x) · · · pN (x) ] , the continuity constraint matrix C is derived from (11) by linear separation of each block. 

Q18. What is the simplest way to determine the flight envelope?

Within such a strict envelope, one would have the largest control-invariant set (or safe set) [27, 28], i.e., the largest set of initial conditions such that, given suitable control input, the aircraft is kept within the flight envelope.