Online Estimation and Adaptive Control for a Class of History Dependent Functional Differential Equations

TL;DR

This paper develops conditions for convergence and stability of online estimation and adaptive control methods for history dependent functional differential equations, relevant to robotic systems with complex nonlinearities.

Contribution

It introduces new theoretical results on convergence and stability for distributed parameter systems with history dependence, including finite-dimensional approximation guarantees.

Findings

Convergence of online estimation methods is established.

Stability of adaptive control strategies is proven using Barbalat's lemma.

Finite-dimensional approximations converge to infinite-dimensional solutions.

Abstract

This paper presents sufficient conditions for the convergence of online estimation methods and the stability of adaptive control strategies for a class of history dependent, functional differential equations. The study is motivated by the increasing interest in estimation and control techniques for robotic systems whose governing equations include history dependent nonlinearities. The functional differential equations in this paper are constructed using integral operators that depend on distributed parameters. As a consequence the resulting estimation and control equations are examples of distributed parameter systems whose states and distributed parameters evolve in finite and infinite dimensional spaces, respectively. suWell-posedness, existence, and uniqueness are discussed for the class of fully actuated robotic systems with history dependent forces in their governing equation of motion. By deriving rates of approximation for the class of history dependent operators in this paper, sufficient conditions are derived that guarantee that finite dimensional approximations of the online estimation equations converge to the solution of the infinite dimensional, distributed parameter system. The convergence and stability of a sliding mode adaptive control strategy for the history dependent, functional differential equations is established using Barbalat's lemma.

Equations

$$M(q(t))\ddot{q}(t)+C(q(t),\dot{q}(t))\dot{q}(t)+\frac{\partial V}{\partial q}=Q_a(t,\mu )+\tau (t)$$

$$\frac{d}{dt}\left[\begin{array}{c}q(t)\\ \dot{q}(t)\end{array}\right]=\left[\begin{array}{cc}0& I\\ -G_0& -G_1\end{array}\right]\left[\begin{array}{c}q(t)\\ \dot{q}(t)\end{array}\right]+\left[\begin{array}{c}0\\ I\end{array}\right](M^{-1}(q)Q_a(t,\mu )+u(t)).$$

$$\frac{d}{dt}\left[\begin{array}{c}e(t)\\ \dot{e}(t)\end{array}\right]=\left[\begin{array}{cc}0& I\\ -G_0& -G_1\end{array}\right]\left[\begin{array}{c}e(t)\\ \dot{e}(t)\end{array}\right]+\left[\begin{array}{c}0\\ I\end{array}\right](M^{-1}(e+q_d)Q_a(t,\mu )+u(t))$$

$$\dot{X}(t)=AX(t)+B((\mathcal{H}X)(t)\circ \mu +u(t)).$$

$$\mathcal{H}:C([0,T),{\mathbb{R}}^m)\to C([0,T),P^{\ast })$$

$$y(t)=(\mathcal{H}X)(t)\circ \mu$$

$$(h_{i}f)(t)\circ {\mu }_i:={\iint }_{\Delta }{\kappa }_i(s,t,f){\mu }_i(s)ds$$

$$(HX)(t):=\left\{\begin{array}{ccc}{h}_1(a(X))(t)& & 0\\ & \ddots & \\ 0& & h_{\ell }(a(X))(t)\end{array}\right\}$$

$$(\mathcal{H}X)(t)=b(X(t))(HX)(t),$$

$$y_i(t):=\sum _{j=1}^l{b}_{ij}(X(t))[h_j(a(X))](t)\circ {\mu }_j$$

$$\displaystyle\kappa(s,t,f):=\left\{\begin{array}[]{ccc}\max\left\{\kappa^{n-1},\gamma_{s_{2}}(f(t))\right\}&&f\text{ increasing on }[t_{n-1},t_{n}],\\ \min\left\{\kappa^{n-1},\gamma_{s_{1}}(f(t))\right\}&&f\text{ decreasing on }[t_{n-1},t_{n}].\end{array}\right.$$

$${\gamma }_{s_2}(\cdot )$$

$${\gamma }_{s_1}(\cdot )$$

$$\displaystyle\phi(x)=\left\{\begin{array}[]{lll}1&\text{if }&\text{}x\in[0,1)\\ 0&\text{otherwise.}\end{array}\right.$$

$${\varphi }_{j,k}(x)=2^{j/2}\varphi (2^{j}x-k)=2^{j/2}{1}_{{\Delta }_{j,k}}(x)$$

$$⟨{\varphi }_{j,k},{\varphi }_{j,l}⟩={\int }_{\mathbb{R}}{\varphi }_{j,k}(x){\varphi }_{j,l}(x)\mathrm{d}x={\delta }_{kl}.$$

$$V_j=\underset{k=0,\cdots ,2^{j-1}}{\mathrm{span}}\{{\varphi }_{j,k}\}.$$

$$\displaystyle\psi(x)=\left\{\begin{array}[]{lll}1&x\in[0,\frac{1}{2}),\\ -1&x\in[\frac{1}{2},1).\end{array}\right.$$

$${\psi }_{j,k}(x)=2^{j/2}\psi (2^{j}x-k),$$

$$V_j=V_{j-1}⨁W_{j-1}.$$

$$\displaystyle 1_{\Delta_{s}}(x)=\left\{\begin{array}[]{lll}1&x\in\Delta_{s}\\ 0&\text{otherwise.}\end{array}\right.$$

$${\varphi }_{j,k}(x)=\frac{1_{{\Delta }_{i_1{i}_2\dots i_j}}(x)}{\sqrt[]{m({\Delta }_{i_1{i}_2\dots i_j})}}$$

$$\Delta:=\left\{(s_{1},s_{2})\in\mathbb{R}^{2}\biggl{|}\underline{s}\leq s_{1}\leq s_{2}\leq\overline{s}\right\}.$$

$$\underset{4^{j}functions}{\underbrace{V_j}}=\underset{4^{j-1}functions}{\underbrace{V_{j-1}}}⨁\underset{3\times 4^{j-1}functions}{\underbrace{W_{j-1}}}.$$

$$\displaystyle\begin{array}[]{lll}(\psi_{j,k},\psi_{m,\ell})=\delta_{j,m}\delta_{k,\ell}&&\text{for all }j,k,m,\ell,\\ (\psi_{j,k},\phi_{m,\ell})=0&&\text{for $j\geq m$ and all $k,\ell$.}\end{array}$$

$${\Pi }_j:P\to V_j.$$

$$\mathcal{A}^{\alpha}_{2}:=\left\{f\in P\biggl{|}\|f\|_{\mathcal{A}^{\alpha}_{2}}:=\left(\sum_{j=0}^{\infty}2^{2\alpha j}\|(\Pi_{j}-\Pi_{j-1})f\|^{2}_{P}\right)^{1/2}\right\}.$$

$${\left\{{\Delta }_{j,k}\right\}}_{k\in {\Lambda }_j}:={\left\{{\Delta }_{i_1\dots i_j}\right\}}_{i_1,\dots ,i_j\in 1,2,3,4}$$

$$(h_{j}f)(t)\circ \mu :={\iint }_{\Delta }\left(\sum _{\ell \in {\Gamma }_j}{1}_{{\Delta }_{j,\ell }}(s)\kappa ({\xi }_{j,\ell },t,f)\right)\mu (s)ds.$$

$$\left|(h_{j}f)(t)\circ \mu -(hf)(t)\circ \mu \right|\le C{2}^{-(\alpha +1)j}$$

$$\left|(h_{j}f)(t)\circ {\Pi }_j\mu -(hf)(t)\circ \mu \right|\le \tilde{C}{2}^{-(\alpha +1)j}$$

Full text

1 Introduction

It is typical in texts that introduce the fundamentals of modeling, stability, and control of robotic systems to assume that the underlying governing equations consist of a set of coupled nonlinear ordinary differential equations. This is a natural assumption when methods of analytical mechanics are used to derive the governing equations for systems composed of rigid bodies connected by ideal joints. A quick perusal of the textbooks [30], [29], or [22], for example, and the references therein gives a good account of the diverse collection of approaches that have been derived for this class of robotic system over the past few decades. Theses methods have been subsequently refined by numerous authors. Over roughly the same period, the technical community has shown a continued interest in systems that are governed by nonlinear, functional differential equations. These methods that helped to define the direction of initial efforts in the study of well-posedness and stability include [23], [20],[21], and their subsequent development is expanded in [12], [26], [25]. More recently, specific control strategies for classes of functional differential equations have appeared in [27], [16], and [14]. The research described in some cases above deals with quite general plant models. These can include classes of delay equations and general history dependent nonlinearities. One rich collection of history dependent models includes hysteretically nonlinear systems. General discussions of nonlinear hysteresis models can be found in [31] or [4], and some authors have studied the convergence and stability of systems with nonlinear hysteresis. For example, a synthesis of controllers for single-input / single-output functional differential equations is presented in [27] and [14], and these efforts include a wide class of scalar hysteresis operators.

The success of adaptive control strategies in classical manipulator robotics, as exemplified by [30], [22], [29], can be attributed to a large degree to the highly structured form of the governing system of nonlinear ordinary differential equations. As is well-known, much of the body of work in adaptive control for robotic systems relies on traditional linear-in-parameters assumptions.

The purpose of this paper is to explore the degree to which the approaches that have been so fruitful in adaptive control of robotic manipulators can be extended to robotic systems governed by certain history dependent, functional differential equations. Emulating the strategy used for robotic systems modeled by ordinary differential equations, we restrict attention to a class of hysteresis operators that satisfy a linear in distributed parameters condition. That is, the contribution to the functional differential equations takes the form of a nonlinear, history dependent operator that acts linearly on an infinite dimensional and unknown distributed parameter.

We illustrate the class of models that are considered in this paper by outlining a variation on two familiar problems encountered in robotic manipulator dynamics, estimation, and control. Consider the task of developing a model and synthesizing a controller for a flapping wing, test robot that will be used to study aerodynamics in a wind tunnel. See [1] for such a system that has been developed by researchers at Brown University over the past few years. Dynamics for a ground based flapping wing robot can be derived using analytical mechanics in a formulation that is tailored to the structure of a serial kinematic chain [22], [30], [29]. The equations of motion take the form

[TABLE]

where $M(q)(t){\in}{\mathbb{R}}^{N{\times}N}$ is the generalized inertia or mass matrix, $C(q(t),\dot{q}(t)){\in}{\mathbb{R}}^{N{\times}N}$ is a nonlinear matrix that represents Coriolis and centripetal contributions, $V$ is the potential energy, $Q_{a}(t,\mu){\in}{\mathbb{R}}^{N}$ is a vector of generalized aerodynamic forces, and ${\tau(t)}{\in}{\mathbb{R}}^{N}$ is the actuation force or torque vector. The generalized forces $Q(t,\mu)$ due to aerodynamic loads are assumed to be expressed in terms of history dependent operators that are carefully discussed below in Section 2, and $\mu$ is the distributed parameter that defines the specific history dependent operator. For the current discussion, it suffices to note that the aerodynamic contributions are unknown, nonlinear, unsteady, and notoriously difficult to characterize.

We consider two specific sets of equations in this paper that are derived from the robotic Equations 1.1, both of which have similar form. We are interested in online identification problems in which we seek to find the final state and distributed parameters from observations of the states of the evolution equation. We are also interested in control synthesis where we choose the input to drive the system to some desired configuration, or to track a given input trajectory. To simplify our discussion, and following the standard practice for many control synthesis problems for robotics, we choose the original control input to be a partial feedback linearizing control that that reformats the control problem in a standard form. In the case of online identification, we choose the input $\tau=M(q)(u-G_{1}\dot{q}-G_{0}q)-(C(q,\dot{q})\dot{q}+\frac{\partial V}{\partial q}(q))$ so that the governing equations take the form

[TABLE]

in terms of a new input $u$. The goal in the online identification problem is to learn the parameters $\mu$ and limiting values $q_{\infty},\dot{q}_{\infty}$ from knowledge of the inputs and states $(u,q,\dot{q})$. We are also interested in tracking control problems. When the desired trajectory is given by $q_{d}$, we choose the input $\tau=M(q)(u+\ddot{q}_{d}-G_{1}\dot{e}-G_{0}e)-(C(q,\dot{q})\dot{q}+\frac{\partial V}{\partial q}(q))$,and the equations governing the tracking error $e:=q-q_{d}$ take the form

[TABLE]

In either of the above two cases, we will show in the next section that the equations can be written in the general form

[TABLE]

where ${A}\in\mathbb{R}^{m\times m}$ is the system matrix, ${B}\in\mathbb{R}^{m\times q}$ is the control input matrix, ${u}(t)\in\mathbb{R}^{q}$ is the corresponding input, and ${(\mathcal{H}X)}(t)$ is a history dependent operator that acts on the distributed parameter ${\mu}$.

2 History Dependent Operators

There is a significant body of research to model and study the unsteady aerodynamic phenomena in flapping flight. Many different models have been presented in the last twenty years to study the aerodynamics and control of flapping flight. Numerically intensive computational fluid dynamics (CFD) presents a precise method to simulate and study the unsteady lift and drag aerodynamic forces. Generally CFD methods exploit high dimensional models that incorporate computationally expensive moving boundary techniques for the Navier-Stokes equations. They are powerful tools to explain some of the characteristics of the aerodynamic forces. One of the characteristics that has inspired the approach here is the history dependence of the aerodynamic lift and drag functions. We refer the interested reader to [35] to study this phenomena in detail. Although CFD methods are advantages in several aspects, they suffer from curse of dimensionality which makes them a very unfavorable choice for online control applications. In this section, we model the unsteady aerodynamics using history dependent operators. Moreover, we present a method that provides an alternative to a high dimensional aerodynamic model that typically evolves in a much lower dimensional space. We also study the accuracy of the presented method with respect to the resolution level of the lower dimensional model.

2.1 A Class of History Dependent Operators

Methods for modeling history dependent nonlinearities can be formulated using a wide array of approaches. Analytical methods for the study of such systems can be based on ordinary or partial differential equations, differential inclusions, functional differential equations, delay differential equations, or operator theoretic approaches. See references [15],[31],[30]. This paper treats evolution equations that are constructed using a specific class of history dependent operators $\mathcal{H}$ that are defined in terms of integral operators constructed from history dependent kernels. These operators are studied in general in [15] and [31]. In this paper the history dependent operators are mappings

[TABLE]

where the $T$ is the final time of an interval under consideration, $m$ is the number of input functions, $q$ is the number of output functions, $P$ is a Hilbert space of distributed parameters and its topological dual space $P^{*}$. We limit our consideration to input$-$output relationships that take the form

[TABLE]

for each $t\in[0,T)$ where $y(t)\in\mathbb{R}^{q}$,$(\mathcal{H}X)(t)\in P^{*}$, and $\mu\in P$.

The definition of $\mathcal{H}$ in this paper is carried out in several steps. All of our history dependent operators $\mathcal{H}$ are defined by a superposition or weighting of elementary hysteresis kernels $\kappa_{i}$ that are continuous as mappings $\kappa_{i}:\Delta\times[0,T)\times C[0,T)\rightarrow C[0,T)$ for $i=1,\ldots,\ell$. We first define the operator $h_{i}:C[0,T)\rightarrow C([0,T),P^{*})$

[TABLE]

for $\mu_{i}\in P_{i}$ and $P=P_{1}\times\cdots\times P_{\ell}$. When we consider problems such as in our motivating examples and numerical case studies, we must construct vectors $H$ of history dependent operators where we define the diagonal matrix

[TABLE]

for each $t\in[0,T)$ where $a:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is some nonlinear smooth map. Finally, our applications to robotics require that we consider

[TABLE]

where $b:\mathbb{R}^{m}\rightarrow\mathbb{R}^{q\times\ell}$ is some nonlinear, smooth map. In terms of our entrywise definitions of the input–output mappings, we have

[TABLE]

for $i=1,\cdots,q$.

In the following discussion, let $\kappa$ be a generic representation of any of the kernels $\kappa_{i}$ for $i=1,\ldots,\ell$. We choose a typical kernel $\kappa(s,t,f)$ to be a special case of a generalized play operator [31]. We suppose that $f$ is a piecewise linear function on $[0,t]$ with breakpoints $0=t_{0}<t_{1}<\cdots<t_{N}=t$. The output function $t\mapsto\kappa(s,t,f)$, for a fixed $s=(s_{1},s_{2})\in\Delta\subset\mathbb{R}^{2}$ and piecewise linear $f:[0,t]\rightarrow\mathbb{R}$, is defined by the recursion where $\kappa^{n-1}:=\kappa(s,t_{n-1},f)$ and for $t\in[t_{n-1},t_{n}]$ we have

[TABLE]

The recursion above depends on the choice of the left and right bounding functions $\gamma_{s_{1}},\gamma_{s_{2}}$ that are depicted in Figure 1. There are given in terms of a single ridge function $\gamma:\mathbb{R}\rightarrow\mathbb{R}$ with

[TABLE]

As noted in [31], the definition of $\kappa$ is extended for any $f\in C[0,T)$ by a continuity and density argument.

2.2 Approximation of History Dependent Operators

The integral operator introduced in Equation 2.2 allows for the representation of complex hysteretic response via the superposition or weighting of fundamental kernels $\kappa_{i}$. These fundamental kernels, each of which has simple input-output relationships, play the role of building blocks for modeling much more complex response characteristics. See [32] for studies of history dependent active materials, [35] for applications that represent nonlinear aerodynamic loading, or Section 6 of this paper to see an example of richness of this class of models. In this section we emphasize another important feature of this particular class of history dependent operators. We show that relatively simple approximation methods yields bounds on the error in approximation of the history dependent operator that are uniform in time and over the class of functions $\mu\in P$.

2.3 Approximation Spaces $\mathcal{A}^{\alpha}_{2}$

The approximation framework we follow in this paper is based on a straightforward implementation of approximation spaces discussed in detail in [11] or [10], and further developed by Dahmen in [6]. We will see that approximation of the class of history dependent operators under consideration exploit a well-known connection between the class of Lipschitz functions and certain approximation spaces as described in [10].

2.4 Wavelets and Approximation Spaces

Multiresolution Analysis ( MRA ) techniques use results from wavelet theory to model multiscale phenomena. To motivate our discussion, we begin by constructing Haar wavelets in one spatial dimension and subsequently discuss how the process can be easily extended to piecewise constant functions over triangulations in two dimensions. The Haar scaling function is defined as follows:

[TABLE]

The dilates and translates $\phi_{j,k}$ of $\phi$ are defined over $\mathbb{R}$ as

[TABLE]

for $j=0,\cdots,\infty$ and $k\in\mathbb{Z}$. It is important to note that with this normalization the functions $\{\phi_{j,k}\}_{k=0}^{2^{j}-1}$ are $L^{2}[0,1]$ orthonormal so that

[TABLE]

In this equation $\Delta_{j,k}=\{x|x\in[2^{-j}k,2^{-j}(k+1)]\}$ and $1_{\Delta_{j,k}}$ is the characteristic function of $\Delta_{j,k}$. For any fixed integer $j$, the $\phi_{j,k}$ form a orthonormal basis that spans the space of piecewise constants over $\{\Delta_{j,k}\}_{k=0}^{2^{j}-1}$. We let $V_{j}$ denote space of piecewise constant functions

[TABLE]

Corresponding to Haar scaling function, the Haar wavelet $\psi$ is defined as

[TABLE]

Again, the translates and dilates of $\psi_{j,k}$ of $\psi$ are given by

[TABLE]

and the complement spaces $W_{j}$ are defined by $W_{j}=\underset{k}{\mathrm{span}}\{\psi_{j,k}\}$. It is straightforward to verify that the spaces $V_{j-1}$ and $W_{j-1}$ form an orthogonal direct sum of $V_{j}$. That is, we have

[TABLE]

This process is well-known and standard in the literature as a means of constructing multiscale bases for $L^{2}[0,1]$. We will follow an analogous strategy to construct multiscale bases over the triangular domain depicted in Figure 2. We first denote the characteristic functions over the triangular domain as shown in the Figure 2

[TABLE]

We next consider the regular refinement shown in Figure 2 where $\Delta_{i_{1}i_{2}}$ is the $i_{2}$ child of $\Delta_{i_{1}}$. In general $\Delta_{i_{1}i_{2}\ldots i_{m}i_{m+1}}$ is the $(m+1)^{st}$ child of $\Delta_{i_{1}i_{2}\ldots i_{m}}$. The multiscaling function $\phi_{j,k}$ is defined as

[TABLE]

where $j$ refers to the level of refinement in the grid.

Since the history dependent operators $(\mathcal{H}X)(t)$ act on the infinite dimensional space $P=P_{1}\times\cdots\times P_{\ell}$ of functions $\mu=(\mu_{1},\ldots,\mu_{\ell})$, we need approximations of these operators for computations and applications. In the discussion that follows we choose each function $\mu_{i}\in P_{i}:=L^{2}(\Delta)$ where the domain $\Delta\subset\mathbb{R}^{2}$ is defined as

[TABLE]

The modification of the construction that follows for different domains $\Delta_{i}$ for the functions $\mu_{i}\in L^{2}(\Delta_{i})$ is trivial, but notationally tedious, and we leave the more general case to the reader. Given the domain $\Delta$ we introduce a regular refinement depicted in Figure 2 and disscused in more detail in Appendix A. The set $\Delta$ is subdivided into $\Delta_{1},\Delta_{2},\Delta_{3},\Delta_{4}$ as shown, and each $\Delta_{i}$ is subdivided into $\Delta_{i1},\Delta_{i2},\Delta_{i3},\Delta_{i4}$. Further subdivision recursively introduces the sets $\Delta_{{i_{1}}\ldots{i_{j}}}$ for $i_{j}=1,\ldots,4$ that are the children of $\Delta_{{i_{1}}\ldots{i_{j-1}}}$.

The characteristic functions $1_{\Delta_{1}},\ldots,1_{\Delta_{4}}$ define a collection of multiscaling functions $\phi^{1},\ldots,\phi^{4}$ as defined in [19]. We define the space of piecewise constant functions $V_{j}$ on grid refinement level $j$ to be the span of the characteristic functions of the sets $\Delta_{i_{1},\ldots,i_{j}}$, so that the dimension of $V_{j}$ is $4^{j}$. We denote by $\left\{\phi_{j,k}\right\}_{k=1,\ldots,4^{j}}$ the orthonormal basis obtained from these characteristic functions on a particular grid level, each normalized so that $(\phi_{j,k},\phi_{j,\ell})_{L^{2}(\Delta)}=\delta_{k,l}$. Each of the basis functions $\phi_{j,k}$ will be proportional to the characteristic function $\phi^{\ell}(2^{j}x+d)$ for some $\ell\in,{1,2,3,4}$ and displacement vector $d$. It is straightforward in this case [19] to define $3$ piecewise constant multiwavelets $\psi^{1},\psi^{2},\psi^{3}$ that that are used to define functions $\psi_{j,m}$ for $m=1,\ldots,3\times 4^{j}$ that span the complement spaces $W_{j}=\text{span}\left\{\psi_{j,m}\right\}_{m=1,\ldots,3\cdot 4^{j}}$ that satisfy

[TABLE]

It is a straightforward exercise to define $L^{2}(\Delta)-$orthonormal wavelets that span $W_{j}$ for each $j\in\mathbb{N}_{0}$, but the nomenclature is lengthly. Since we do not use the wavelets specifically in this paper, the details are omitted. Each function $\psi_{j,m}$ is proportional to one of the three scaled and translated multiwavelet functions and satisfies the orthonormality conditions

[TABLE]

In the next step, we denote the orthogonal projection onto the span of the piecewise constants defined on a grid of resolution level $j$ by $\Pi_{j}$ so that

[TABLE]

Finally, we define the approximation space $\mathcal{A}^{\alpha}_{2}$ in terms of the projectors $\Pi_{j}$ as

[TABLE]

Note that this is a special case of the more general analysis in [6]. We define our approximation method in terms of one point quadratures defined over the triangles $\Delta_{i_{1}\ldots i_{j}}$ that constitute the grid of level $j$ that defines $V_{j}$. For notational convenience, we collect all triangles at a fixed level $j$ in the singly indexed set

[TABLE]

where $\Lambda_{j}:=\left\{k\in\mathbb{N}\quad|1\leq k\leq 4^{j}\right\}$, and the quadrature points are chosen such that $\xi_{j,k}\in\Delta_{j,k}$ for $k=1,\cdots,\Lambda_{j}$. We now can state our principle approximation result for the class of history dependent operators in this paper.

Theorem 1**.**

Suppose that the function $\gamma$ that defines the history dependent kernel in Equation 2.5 is a bounded function in $C^{\alpha}(\mathbb{R})$, and define the approximation $h_{j}$ associated with the grid level $j$ of the history dependent operator $h$ to be

[TABLE]

Then there is a constant $C>0$ such that

[TABLE]

for all $f\in C[0,T]$, $t\in[0,T]$, and $\mu\in P$. If in addition $\mu\in\mathcal{A}^{\alpha+1}_{2}$, there is a constant $\tilde{C}>0$ such that

[TABLE]

for all $f\in C[0,T]$ and $t\in[0,T]$.

Proof.

We first prove the inequality in Equation 2.7. By definition of the operator $h$, we can write

[TABLE]

Since the ridge function $\gamma$ is a bounded function in $C^{\alpha}(\mathbb{R})$, the output mapping $s\mapsto\kappa(s,t,f)$ is also a bounded function in $C^{\alpha}(\Delta)$ where the Lipschitz constant is independent of $t\in[0,T]$ and $f\in C[0,T]$. Using Proposition 2.5 of [31], we have

[TABLE]

Since we have

[TABLE]

the second inequality in Equation 2.8 follows from the first Equation 2.7 provided we can show that

[TABLE]

for some constant $C$. But it is a standard feature of the approximation spaces that if $\mu\in\mathcal{A}^{\alpha+1}_{2}$, then $\|\mu-\Pi_{j}\mu\|_{P}\leq 2^{-(\alpha+1)j}\|\mu\|_{\mathcal{A}^{\alpha+1}_{2}}$. To see why this is so, suppose that $\mu\in\mathcal{A}^{\alpha+1}_{2}$. We have

[TABLE]

When we apply this to our problem, the upper bound follows immediately

[TABLE]

since the boundedness of the ridge function $\gamma$ implies the uniform boundedness of the history dependent operators $(h_{j}f)(t)$ over $[0,T]$. ∎

Theorem 1 can now be used to establish error bounds for input-output maps that have the form in Equation 2.4.

Theorem 2**.**

Suppose that the hypotheses of Theorem 1 hold. Then we have

[TABLE]

Where $\mathcal{H}$ is defined in Equations 2.1,2.4, and $\mathcal{H}_{j}$ is defined in Equations 2.9, 2.10 and 2.11 below.

Proof.

Recall that for $i=1\dots q$ we had

[TABLE]

In matrix form this equation can be expressed as

[TABLE]

It follows that,

[TABLE]

By assumption $X\in C([0,T],\mathbb{R}^{m})$. The construction of $H$ and $\mathcal{H}$ guarantees that $(HX)(t)\in\mathcal{L}(P,\mathbb{R}^{l}),$

$H:C\left([0,T],\mathbb{R}^{m}\right)\to C\left([0,T],\mathcal{L}(P,\mathbb{R}^{l})\right),$ and $\mathcal{H}:C([0,T],\mathbb{R}^{m})\to C([0,T],\mathbb{R}^{q}).$ In this proof we denote by $(\mathbb{R}^{l},\|.\|_{u})$ the norm vector space that endows $\mathbb{R}^{l}$ with the $l^{m}$ norm $\|v\|_{u}:=\left(\sum^{l}_{i=1}|v_{i}|^{u}\right)^{\frac{1}{u}}$ for $1\leq u\leq\infty$. The normed vector space $(\mathbb{R}^{q\times l},\|.\|_{s,u})$ denotes the induced operator norm on matrices that map $(\mathbb{R}^{l},\|.\|_{u})$ into $(\mathbb{R}^{q},\|.\|_{s})$. Now we define an approximation on the mesh level $j$ of $\mathcal{H}$ to be

[TABLE]

[TABLE]

and

[TABLE]

for $i=1\dots\ell$ and $\nu_{i}\in P_{i}$. To simplify the derivation or an error bound for approximation of $\mathcal{H}X(t)\circ\mu$, let $(HX)(t)\circ\mu$ be denoted by $g(t)$. We have assumed that $X\to b(X)$ and $t\to X(t)$ are continuous mappings. There fore $t\mapsto b(X(t))$ is continuous and on a compact set $[0,T]$, and $b(X(\cdot))$ $\in$ $C([0,T],\mathbb{R}^{q\times l})$. We therefore by definition have

[TABLE]

with the norms explicitly denoted in the subscript. For $t\in[0,T]$, and applying these definitions,

[TABLE]

with

[TABLE]

Therefore we can now write

[TABLE]

Hence, recalling Theorem 1 we can now derive the convergence rate

[TABLE]

Therefore we obtain the final bound

[TABLE]

for all $t\in[0,T]$. ∎

3 Well-Posedness: Existence and Uniqueness

The history dependent governing equations studied in this paper are a special case of the more general class of abstract Volterra equations or functional differential equations. A general treatise on abstract Volterra equations can be found in [5], while various generalizations of theory for the existence and uniqueness of functional differential equations have been given in [12], [25], [16]. We have noted in Section 1 that the general form of the governing equations we consider in this paper have the form

[TABLE]

where the state vector $X(t)\in\mathbb{R}^{m}$, the control inputs $u(t)\in\mathbb{R}^{q}$, $A\in\mathbb{R}^{m{\times}m}=\mathbb{R}^{2n{\times}2n}$ is a Hurwitz matrix, and $B\in\mathbb{R}^{m{\times}q}$ is the control input matrix. We make the following assumptions about the history dependent operators $\mathcal{H}$:

1. H1)

$\mathcal{H}:C([0,\infty),\mathbb{R}^{m})\mapsto C([0,\infty),P^{*})$ 2. H2)

$\mathcal{H}$ is causal in the sense that for all $x,y\in C([0,\infty);\mathbb{R}^{m})$,

[TABLE] 3. H3)

Define the closed set consisting of all continuous functions $f$ that remain within radius $r$ of the initial condition $X_{0}$ over the closed interval $[t,t+h]$,

[TABLE]

for a fixed $X_{0}\in\mathbb{R}^{m}$. For each $t\geq 0$, we assume that there exist $h,r,L>0$ such that

[TABLE]

for all $X,Y\in\overline{\mathcal{B}}_{[t,t+h],r}(X_{0})$.

Our first result guarantees the existence and uniqueness of a local solution to Equation 1.4, and also describes an important case when such local solutions can be extended to $[0,\infty)$. This theorem can be proven via the existence and uniqueness Theorem 2.3 in [16] for functional delay-differential equations. However, since we are not interested in delay differential equations in this paper, but rather on a highly structured class of integral hysteresis operators, the proof can be much simplified.

Theorem 3**.**

Suppose that the history dependent operator $\mathcal{H}$ satisfies the hypotheses (H1),(H2),(H3). Then there is a $\delta>0$ such that Equation 3.1 has a solution $X\in C([0,\delta),\mathbb{R}^{m})$. Suppose the interval $[0,\delta)$ is extended to the maximal interval $[0,\omega)\subset[0,\delta)$ over which such a solution exists. If the solution is bounded, then $[0,\omega)=[0,\infty)$.

Corollary 1**.**

Suppose that the history dependent operator $\mathcal{H}$ in Equation 3.1 is defined as in Equation 2.3 and 2.4 in terms of a globally Lipschitz, bounded continuous ridge function $\gamma:\mathbb{R}\rightarrow\mathbb{R}$ in Equation 2.5. Then Equation 3.1 has a unique solution $X\in C([0,\infty),\mathbb{R}^{m})$ for each $\mu\in P$.

Proof.

For completeness, we outline a simplified version the proof of Theorem 3 for our class of history and parameter dependent equations. As a point of comparison, the reader is urged to compare the proof below to the conventional proof for systems of nonlinear ordinary differential equations, such as in [17]. If we integrate the equations of motion in time, we can define an operator $T:C([0,h),\mathbb{R}^{m})\rightarrow C([0,h),\mathbb{R}^{m})$ from

[TABLE]

for all $t\in[0,h]$. As introduced in hypothesis (H3), we select $h,r>0$ and define

[TABLE]

such that the local Lipschitz condition in Equation 3.2 holds. Now we consider restricting the equation to a subinterval $[0,\delta]\subseteq[0,h]$, and investigate conditions on $T$ that enable the application of the contraction mapping theorem. We first study what conditions on $\delta>0$ are sufficient to guarantee that $T:\overline{\mathcal{B}}_{[0,\delta],r}(X_{0})\rightarrow\overline{\mathcal{B}}_{[0,\delta],r}(X_{0})$. We have

[TABLE]

where $M_{T}=M_{H}+M_{u}$. Now we restrict $\delta$ so that

[TABLE]

which implies

[TABLE]

We thereby conclude that

[TABLE]

and it follows that $T:\overline{\mathcal{B}}_{[0,\delta],r}\rightarrow\overline{\mathcal{B}}_{[0,\delta],r}$. Next we study conditions on $\delta$ that guarantee that $T:\overline{\mathcal{B}}_{[0,\delta],r}\rightarrow\overline{\mathcal{B}}_{[0,\delta],r}$ is a contraction. We compute directly a bound on the difference of the output as

[TABLE]

If we choose

[TABLE]

it is apparent that $T$ is a contraction that maps the closed set $\overline{\mathcal{B}}_{[0,\delta],r}$ into itself. There is a unique solution in $\overline{\mathcal{B}}_{[0,\delta],r}$ on $[0,\delta]$. ∎

4 Online Identification

A substantial literature has emerged that treats online estimation problems for linear or nonlinear plants governed by systems of ordinary differential equations. Approaches for these finite dimensional systems that are based on variants of Lyapunov’s direct method can be found in any of a number of good texts including, for instance, [24], [28], or [18]. The general strategies that have proven fruitful for such finite dimensional systems have often been extended to classes of systems whose dynamics evolve in an infinite dimensional space: distributed parameter systems. A discussion of the general considerations for identification of distributed parameter systems can be found in [2], for example, while studies that are specifically relevant to this paper include [7], [8], [9], and [3].

In this section we adapt the framework introduced in [3] to our class of history dependent, functional differential equations. The approach in [3] assumes that the state equations for the distributed parameter system have first order form, and they are cast in terms of a nonlinear, parametrically dependent bilinear form that is coercive. The resulting equations that govern the error in state and in distributed parameter estimates is a nonlinear function of the state trajectory of the plant. In contrast, a similar strategy in this paper yields error equations that depend nonlinearly on the history of the state trajectory.

The general online estimation problem discussed in this section assumes that we observe the value of the state $X(t)\in\mathbb{R}^{m}$ at each time $t\geq 0$ that depends on some unknown distributed parameter $\mu\in P$, and subsequently use the observed state to construct estimates $\hat{X}$ of the states and $\hat{\mu}$ of the distributed parameters. We construct online estimates that evolve on the state space $\mathbb{R}^{m}\times P$ according to the time varying, distributed parameter system equations

[TABLE]

for $t\geq 0$ where the initial conditions are $\hat{X}_{0}:=X_{0}$, $\hat{\mu}(0):=\mu_{0}$. In these equations, we denote the adjoint operator $L^{*}$ for any bounded linear operator $L$. These equations can be understood as incorporating a natural choice of a parameter update law. The learning law above can be interpreted as generalization of the conventional gradient update law that features prominently in approaches for finite dimensional systems [18] and that has been extended to distributed parameter systems in [3]. It is immediate that the error in estimation of the states $\tilde{X}:=X-\hat{X}$ and in the distributed parameters $\tilde{\mu}:=\mu-\hat{\mu}$ satisfy the homogeneous system of equations

[TABLE]

4.1 Approximation of the Estimation Equations

The governing system in Equations 4.1 constitute a distributed parameter system since the functions $\hat{\mu}(t)$ evolve in the infinite dimensional space $P$. In practice these equations must be approximated by some finite dimensional system. We define $\tilde{X}_{j}=\hat{X}-\hat{X}_{j}$ and $\tilde{\mu}_{j}=\hat{\mu}-\hat{\mu}_{j}$ where $\tilde{X}_{j}$ and $\tilde{\mu}_{j}$ express approximation errors due to projection of solutions in $\mathbb{R}^{m}\times P$ to a finite dimensional approximation space. We construct a finite dimensional approximation of the the online estimation equations using the results of Section 2.2 and obtain

[TABLE]

Theorem 4**.**

Suppose that the history dependent operator $\mathcal{H}$ in Equation 3.1 is defined as in equation 2.3 and 2.4 in terms of a globally Lipschitz, bounded continuous ridge function $\gamma:\mathbb{R}\rightarrow\mathbb{R}$ in Equation 2.5. Then for any $T>0$, we have

[TABLE]

as $j\rightarrow\infty$.

Proof.

Define the operators $G(t):P\rightarrow\mathbb{R}^{m}$ and $G_{j}(t):P\rightarrow\mathbb{R}^{m}$ for each $t\geq 0$ as

[TABLE]

The time derivative of the error in approximation can be expanded as follows:

[TABLE]

We will next use a common inequality that can be derived from two applications of the triangle inequality. We have

[TABLE]

We conclude from this pair of inequalities that

[TABLE]

The specific form that we apply this theorem is written as

[TABLE]

We apply the inequality in Equation 4.4 to each term in which $\tilde{\mu}_{j}$ and $\tilde{X}_{j}$ appear in a product.

[TABLE]

Then

[TABLE]

We integrate this inequality in time from [math] to $t$ to obtain

[TABLE]

Choose $a,b>0$ large enough so that $(2c+a+b)>0$ and set $\gamma>1$. If we define

[TABLE]

then the inequality can be written as

[TABLE]

Gronwall’s Inequality now completes the proof of the theorem (see Appendix C). ∎

We also further investigate $\lambda_{j}(t)$ to derive the convergence rate for the approximate states and parameters evolving associated with level $j$ resolution. According to the convergence results obtained in Theorem 1 we have $\|G(s)-G_{j}(s)\|_{\mathcal{L}(P,\mathbb{R}^{m})}=\|B(\mathcal{H}X)(t)-B(\mathcal{H}_{j}X)(t)\Pi_{j}\|\leq C_{2}2^{-(\alpha+1)j}$. Therefore, $\|G(s)-G_{j}(s)\|^{2}\leq C_{2}^{2}2^{-(\alpha+1)2j}$. It then follows that

[TABLE]

If $t$ $\simeq$ $C_{3}2^{(\alpha+1)j}$, then $\lambda_{j}(t)<\mathcal{O}(2^{-(\alpha+1)j})\quad\text{for}\quad t\in[0,C_{3}2^{(\alpha+1)j}].$

5 Adaptive Control Synthesis

In order to estimate the function $\mu$ that weighs the contribution of history dependent kernels to the equations of motion, we first map it to an n-dimensional subspace of square integrable functions using a projection operator $\Pi^{n}:P\mapsto P^{n}$. Let

[TABLE]

be the governing equation of a robotic system after applying a feedback linearization control signal as mentioned in Equation 1.3 with $u=v-(\mathcal{H}X)\circ\hat{\mu}$. We substitute $\mu=\Pi^{n}\mu+(I-\Pi^{n})\mu$ and write

[TABLE]

Finally, by replacing $d=\{(\mathcal{H}X)(I-\Pi^{n})\circ\mu\}$ we obtain

[TABLE]

where

[TABLE]

Theorem 5**.**

Suppose the state equations have the form of Equation 1.3 and the matrix $\mathcal{P}$ is a symmetric positive definite solution of the Lyapunov equation $A^{T}\mathcal{P}+\mathcal{P}A=-Q$ where $Q>0$. Then by employing the update law $\dot{\tilde{\mu}}=-((\mathcal{H}X)\Pi^{n})^{*}B^{T}\mathcal{P}X$, the control signal

[TABLE]

with $k>\|d\|$ drives the tracking error dynamics of the closed loop system is uniformly ultimately bounded and its norm is eventually $O(\epsilon)$.

Proof.

We choose the Lyapunov function

[TABLE]

where $\mathcal{P}$ is the solution of the Lyapunov equation $A^{T}\mathcal{P}+\mathcal{P}A=-Q$. The derivative of the Lyapunov function $V$ along the closed loop system trajectory is

[TABLE]

Therefore we have

[TABLE]

By Theorem 4.18 in [17] we conclude that there is a $\bar{T}>0$ and $\tau>0$ such that $\|X(t)\|\leq\bar{C}\epsilon$ for all $t\geq\bar{T}$.

∎

6 Numerical Simulations

Our principle approximation result, the proposed online identification, and adaptive control of systems with history dependent forces are verified in this section. In the first experiment, we validate the operator approximation error bound presented in Theorem 1. In the second experiment, we model a wind tunnel single wing section with a leading and trailing edge flaps and apply the proposed sliding mode adaptive controller presented in Theorem 4. We illustrate the stability of the closed loop system and convergence of the closed-loop system trajectories to the equilibrium point.

6.1 Operator Approximation Error

In this section we consider a collection of numerical experiments to validate the operator approximation rates derived in Theorem 1. In order to show that Equation 2.8 holds, we choose a function $\mu(s)$ over $\Delta$ and then calculate $(h_{j}f)(t)\circ\mu_{j}$ for different levels of refinement. Since the computation of $(hf)(t)\circ\mu$ exactly is numerically infeasible, we choose $J\gg j$ as the finest level of refinement in our simulation. According to Theorem 1, we have

[TABLE]

and for $j\ll J$ we see that

[TABLE]

Assuming $C=\max\{C_{j},C_{J}\}$ and using the triangle inequality, we obtain

[TABLE]

Therefore, given the weights $\mu_{J}$ for the finest level of refinement $J$, we can evaluate $\mu_{j}=\Pi_{j}\mu_{J}$ and numerically verify Equation 6.1.

Figure 3 shows the simulation results for $J=7$ and $j=2,3,4,5$. The error term attenuates with increasing j. In order to investigate the rate of attenuation, we evaluate constant $C$ for different levels of refinements. As shown in figure 4, $C$ is approximately constant with respect to $j$ which agrees with the result from Equation 6.1.

6.2 Online Identification of History Dependent Aerodynamics and Adaptive Control for a Simple Wing Model

The reformatted governing equations of the system take the form of Equation 1.2 where $Q_{a}(t,\mu)$ is the vector of generalized history dependent aerodynamic loads. The dynamic equation of the system can be written in the form of Equation 1.4, where the history dependent term $M^{-1}(q)Q_{a}(t,\mu)$ is rewritten in terms of a history dependent operator $(\mathcal{H}X)(t)$ acting on the distributed parameter function $\mu$. The history dependent operator includes a family of fixed history dependent kernels and the distributed parameters $\mu$ act as a weighting vector that determines the contribution of a specific history dependent kernel to the overall history dependent operator.

We perform an offline identification based on a set of experimental data collected from a wind tunnel experiments or CFDsimulations. These define a nominal model for the history dependent aerodynamic loads that appear in the governing equations of the system. We can exploit the model in the numerical simulations to perform an online estimation of the history dependent aerodynamics and adaptive control of a simple wing model. The details of offline identification of history dependent aerodynamics follow the steps explained in [35].

The model developed in Figure 5 is chosen to validate our proposed adaptive sliding mode controller where $w$ is the velocity of wind, $k_{h}$ is spring constant in plunge, $k_{\theta}$ is a spring constant in pitch, $\theta$ is the pitch angle, $h$ is the plunge displacement, $c_{\theta}$ and $c_{h}$ are viscous damping coefficients, $m$ and $I_{\theta}$ are the mass and moment of inertia and, $x_{\theta}$ is the non-dimensionalized distance between center of mass and the elastic axis. Finally, $L$ and $M$ are lift and moment generated by the leading and trailing edge flaps. The angles $\beta_{1}$ and $\beta_{2}$ define the rotation of the trailing edge and leading edge flaps respectively. The dynamic equations of the wing model is derived in the appendix C as

[TABLE]

We have assumed the aerodynamic moment $M$ to be zero and the distance $x_{a}$ between the aerodynamic center $A$ and hinge point to be negligible to simplify the simulation. The unsteady aerodynamic lift is $L=Q_{a}(t,\mu)$ where $Q_{a}(t,\mu)=(\mathcal{H}X)\circ\mu$ reflects the history dependent nature of aerodynamic loads. We rewrite Equation 6.2 to achieve the standard form presented in Equation 1.2.

The adaptive controller presented in Theorem 4 is composed of two parts. The first part compensates for the flutter generated by the history dependent aerodynamic forces through online identification of the aerodynamics. The second part employs an sliding mode controller to compensate for modeling errors.

It is noteworthy that numerical time integration of the evolution equations must accommodate history dependent terms. Since the dynamics of such systems are given via functional differential equations, the ordinary integration rules are not directly applicable. We exploit the predictor-corrector integration rule that has been introduced first in [37]. We also refer the interested reader to our previous paper [39] for details of such integration rules.

Figure 6 Shows the simulation results for the case where $\epsilon=0.01$ and $t_{h}=0.001$. The system response eventually enters in a $\epsilon$ neighborhood of the sliding manifold. However, as depicted in the figure a chattering behavior occurs in the control signal and system trajectories. We trace this behavior back to the integration error induced by the size of time step. When we increase $\epsilon$ or reduce the integration time step, the control signal and system trajectories become smooth. The simulation results for $\epsilon=0.01$ and $t_{h}=0.0005$ are depicted in Figure 7. The system trajectories converge to a neighborhood of zero or the set $\mathcal{M}$ in Equation LABEL:eq:M with time and the control signals are relatively smooth. Also, Figure 8 shows the case when $\epsilon=0.1$ and $t_{h}=0.001$. The convergence rate of the signals to zero is slower but the results do not show any chattering. Therefore, the proposed smooth sliding mode adaptive controller proves to be effective to identify and compensate for the unknown history dependent aerodynamic forces.

7 Results and Conclusion

In this paper, we have derived an explicit bound for the error of approximation for certain history dependent operators that are used in construction of robotic FDE’s in  [35] and this paper. The numerical simulations presented validate our results. We establish uniform upper bounds on their accuracy of the approximations. The uniform $\mathcal{O}(2^{-(\alpha+1)j})$ rates of approximation for grid resolution $j$ depend on the Holder coefficient $\alpha$ that describes the smoothness of the ridge functions that define the history dependent kernels. In Section 3 we prove the existence and uniqueness of a local solution for the special case of functional differential equations with history dependent terms shown in Equation 3.1. Since the functional differential equation of interest evolves in an infinite dimensional space, we construct finite dimensional approximations with grid resolution $j$. We further show that the solution of the finite dimensional distributed parameter system converges to the solution of the infinite dimensional FDE as the resolution is refined. Finally, we propose an adaptive control strategy to identify and compensate the unknown history dependent dynamics.

Appendix A: Wavelets and Approximation Spaces over the Triangular Domain

We define the multiscaling functions

[TABLE]

in which

[TABLE]

and $m(\Delta_{i_{1},i_{2},\ldots,i_{j}})$ is the area of a triangle in the level $j$ refinement. We have defined $(hf)(t)\circ\mu=\iint_{\Delta}\kappa(s,t,f)\mu(s)ds.$ The approximation $(h_{j}f)(t)\circ\mu$ of this operator is given by

[TABLE]

where $\xi_{j,l}$ is the quadrature point of number $l$ triangle of grid level $j$. We approximate $\mu(s)\approx\sum_{m\in\Gamma_{j}}\mu_{j,m}\phi_{j,m}(s)$. Therefore,

[TABLE]

For an orthonormal basis $\left\{\phi_{k}\right\}_{k=1}^{\infty}$ of the separable Hilbert space $P$, we define the finite dimensional spaces for constructing approximations as $P_{n}:=\text{span}\left\{\phi_{k}\right\}_{k=1}^{n}$. The approximation error $E_{n}$ of $P_{n}$ is given by

[TABLE]

The approximation space $\mathcal{A}^{\alpha}_{2}$ of order $\alpha$ is defined as the collection of functions in $P$ such that

[TABLE]

For our purposes, the approximation spaces are easy to characterize: they consist of all functions $f\in P$ whose generalized Fourier coefficients decay sufficiently fast. That is, $f\in\mathcal{A}^{\alpha}_{2}$ if and only if

[TABLE]

for some constant $C$.

Appendix B: The Projection Operator $\Phi_{J\rightarrow j}$

The orthogonal projection operator $\Phi_{J\rightarrow j}:V_{J}\rightarrow V_{j}$ maps a distributed parameter $\mu_{J}$ to $\mu_{j}$ i.e. $\Phi_{J\rightarrow j}:\mu_{J}\mapsto\mu_{j}$.

By exploiting the orthogonality property of the operator we have

[TABLE]

Therefore, we can write

[TABLE]

Since orthogonality implies $\iint_{\Delta}\phi_{j,m}(s)\phi_{j,n}(s)ds=\delta_{m,n},$ we conclude that

[TABLE]

From Theorem 1 we have

[TABLE]

with

[TABLE]

where $\mu\in P=L^{2}(\Delta)$ and we approximate $\mu(s)\approx\sum_{l\in\Gamma_{J}}\mu_{J,l}\phi_{J,l}(s)\in V_{J}$. To implement this for the finest grid $J$, we compute

[TABLE]

when $\sqrt[]{m(\Delta_{J,l})}$ is the area of the corresponding triangle $\Delta_{J,l}$ in the grid having resolution level $J$.

Appendic C: Gronwall’s Inequality

We employ the integral form of Gronwall’s Inequality to obtain our final convergence result. Many forms of Gronwall’s Inequality exist, and we will use a particularly simple version. See Section 3.3.4 in [18]. If the piecewise continuous function $f$ satisfies the inequality

[TABLE]

with some piecewise continuous functions $\alpha,\beta$ where $\alpha$ is nondecreasing, then

[TABLE]

Appendix D: Modeling of a Prototypical Wing Section

Figure 5 shows a simplified model of the wing. In the figure we denote the center of mass by $c.m.$, $A$ is the aerodynamic center, and $O$ is the elastic axis of the wing. The constants $K_{h}$ and $K_{\theta}$ are the linear and torsional stiffness, and $h$ is the distance from origin to point $O$ in the fixed reference frame. We denote by $x_{\theta}$ the distance between point $O$ and center of mass, whereas $x_{a}$ is the distance between $O$ and $A$. Point $O$ is the origin for the body fixed reference frame.

We employ The Euler-Lagrange technique to derive the equation of motion for the depicted wing model. The function $L(\theta,\dot{\theta})$ is the history dependent lift force acting at the aerodynamic center, and $M(\theta,\dot{\theta})$ is the history dependent aerodynamic moment about point $A$. The variables $L_{\beta_{1}}$ and $L_{\beta_{2}}$ are the actuating forces acting at point $D$, and $\beta_{1}$, $\beta_{2}$ are the angles between the midchord of the wing and the trailing edge and leading edge flaps, respectively.

The position vector of the mass center is given as

[TABLE]

and therefore the corresponding velocity of point $C$ is

[TABLE]

The rotation matrix for transformation between inertial frame of reference to body fixed frame of reference is

[TABLE]

The kinetic energy is computed to be

[TABLE]

[TABLE]

and the corresponding potential energy is

[TABLE]

therefore we can write Lagrangian as $L=T-V$. We apply Euler-Lagrange equations to write the equation of motion as follows

[TABLE]

The above equation is written in the form of a standard robotic equations of motion $M(q(t))\ddot{q}(t)+C(q(t),\dot{q}(t))\dot{q}(t)+K(q(t))=Q_{a}(t)+\tau(t)$, where $q=[h\>\theta]^{T}$. We have discussed control applications for such systems in detail in Section 1. In addition, we employ a simplified version of this equation to validate our online identification and adaptive control strategy in Section 6.2.