Asymmetric Airfoil Morphing via Deep Reinforcement Learning

Abstract

Morphing aircraft are capable of modifying their geometry configurations according to different flight conditions to improve their performance, such as by increasing the lift-to-drag ratio or reducing their fuel consumption. In this article, we focus on the airfoil morphing of wings and propose a novel morphing control method for an asymmetric deformable airfoil based on deep reinforcement learning approaches. Firstly, we develop an asymmetric airfoil shaped using piece-wise Bézier curves and modeled by shape memory alloys. Resistive heating is adopted to actuate the shape memory alloys and realize the airfoil morphing. With regard to the hysteresis characteristics exhibited in the phase transformation of shape memory alloys, we construct a second-order Markov decision process for the morphing procedure to formulate a reinforcement learning environment with hysteresis properties explicitly considered. Subsequently, we learn the morphing policy based on deep reinforcement learning techniques where the accurate information of the system model is unavailable. Lastly, we conduct simulations to demonstrate the benefits brought by our learning implementations and validate the morphing performance of the proposed method. The simulation results show that the proposed method provides an average 29.8% performance improvement over traditional methods.

1. Introduction

While unmanned aerial vehicles (UAVs) have played a crucial role in various civil and military missions, studies demonstrate that birds usually possess higher flight maneuverability and agility than comparatively-sized aircraft in complex and varying environments [1,2]. One of the critical advantages of birds is that they morph their wings and tails intricately to perform efficient behaviors including perching, hovering and maintaining stability under different flight conditions [3]. Such aerodynamic adaptability has aroused flourishing interest in the design and control of avian-inspired morphing UAVs [4,5,6]. In this work, we focus on aircraft capable of morphing their wings by modifying the geometric configuration of the airfoil shape, which refers to more specifically camber morphing [7]. Bird wings are usually cambered to generate sufficient lift force at a low angle of attack. The camber does not remain constant through their flight, and many observations show that birds actively control the camber of their proximal wing via remiges and modify the distal wing airfoil shape via their primary feathers [8,9]. These investigations provide insight into the study of airfoil-morphing aircraft. The benefits brought by camber morphing for aircraft include increasing lift, reducing drag and airframe noise mitigation [10]. Applications of camber morphing mechanisms include a lead-edge morphing combined smart droop nose design, which achieves high-lift performance with significantly reduced complexity and mass [11], and a flexible morphing trailing edge design with deformable ribs, which is used to enhance the Fowler flaps and act as a substitution for ailerons for civil transport aircraft [12].

The ideal airfoil is usually generated using shape optimization techniques including gradient-based methods and gradient-free methods [13], which optimize some aerodynamic performance parameters of aircraft, such as the drag coefficient, the lift coefficient and the lift-to-drag ratio [14,15,16]. Recently, deep reinforcement learning (DRL) approaches such as proximal policy optimization (PPO) [17] have been exploited to learn the airfoil shape directly according to performance metrics computed by computational fluid dynamics (CFD) solvers [18]. In [19], a 3D-printed morphing airfoil model is developed, and the optimal configuration is generated via Q-learning to match the desired pitching moment. In [20], the transfer learning technique is combined with DRL for shape optimization, which formulates a multi-fidelity framework and reduces the computational cost.

Although an optimized airfoil shape can be calculated given a certain flight condition, it is challenging to morph into optimal shapes instantaneously during the flight procedure due to the uncertainty and inconsistency of the environments and aircraft dynamics [2,21]. This galvanizes the utilization of data-driven methods including deep learning and reinforcement learning for morphing control. In [22], the morphing air vehicle was modeled as a smart block in the shape of a rectangular parallelopiped, and the control policy was learned by actor-critic methods. This framework was extended to an ellipsoid-shaped aircraft with Q-learning methods to produce a more efficient policy [23]. The thickness and camber of the airfoil shape were adjusted via Q-learning in [24], where the rewards were related to aerodynamic parameters. The constant-strength source doublet panel method was utilized in [25] to calculate the aerodynamic forces on the morphing air vehicle. In [26], airfoil morphing by vertically moved control points was designed, where the unknown drag and lift coefficients were estimated using neural networks.

Concerning the realization of morphing mechanisms, biologically inspired mechanical joints are often adopted to control sweep-, dihedral- and twist-morphing wings [27,28,29]. However, the development and actuation for camber-morphing aircraft are more challenging, as they are required to permit smooth transitions in the airfoil shape [2]. Conventional actuators for camber morphing include servo motors and hydraulic actuators [30,31], where the wing is constructed by articulated rigid-linked components [32] or a deformable skin with internal compliant mechanisms [33]. Recent advances in material technologies have led to a proliferation of applications of smart materials, especially shape memory alloys (SMA), for morphing aircraft [34,35]. Shape memory alloys are metallic alloys capable of transforming their crystalline structures between two phases to deform and recover their shapes through heating and cooling, which is suitable for morphing actuators because of their properties, including high-actuation energy densities and large recoverable strains [35]. A morphing wing under subsonic cruise flight conditions was developed in [36] using a flexible skin and a group of SMA actuators, which reduced the fuel consumption. In [37], the authors designed an autonomous morphing helicopter rotor blade where the SMAs equipped on the blade modify the section camber according to ambient temperature. In [38], a super critic airfoil actuated by SMAs was developed, and the transonic aerodynamics are investigated. Nevertheless, to actively control the strain of SMA wires is not straightforward due to the temperature hysteresis exhibited during the phase transformation [39]. This gives rise to increasing investigations of learning-based control for SMAs. Reinforcement learning methods including Q-learning and Sarsa have been applied to adjust the strain of SMAs via resistive heating, where the hysteretic dynamics are modeled as hyperbolic tangent curves [40,41]. However, there is no explicit consideration for the hysteresis properties in the policy-learning procedure of these studies.

In this work, the morphing control for a deformable asymmetric airfoil based on deep reinforcement learning techniques is investigated. The contributions of this work are threefold, as follows. Firstly, the shape of the asymmetric airfoil is designed based on Bézier curves, and the morphing mechanism of such airfoils is modeled via SMA wires. Subsequently, a dynamic system between input voltages and airfoil shapes is developed, which characterizes the hysteresis behaviors in the manner of Markov decision processes. Finally, the morphing policy is constructed based on deep reinforcement learning approaches without accurate knowledge of system models, which adjusts the airfoils during flight procedures to track the optimal shapes in different flight conditions.

The rest of this paper is organized as follows. In Section 2, the morphing method is developed in three steps. In Section 2.1, the asymmetric airfoil is shaped using piece-wise Bézier curves. In Section 2.2, the morphing airfoil is modeled by SMAs, and a dynamic model is derived. In Section 2.3, a deep reinforcement learning-based morphing policy is developed. In Section 3, simulations are conducted to validate the proposed methods. In Section 4, this work is summarized.

2. Materials and Methods

2.1. Asymmetric Airfoil Shape Modeling

Well-known methods of curve synthesis for airfoil design and optimization include splines (e.g., B-spline and Bézier curves) [42,43], free-form deformation (FFD) [44], and class-shape transformations (CST) [45]. In this work, Bézier curves are selected for their straightforward design procedure and simple calculations. The shape of the asymmetric airfoil is parameterized via N control points, which are connected by Bézier curves [18,20]. To generate an untangled shape, the control points are distributed in an annulus with a predefined inner radius $R_1$ and outer radius $R_2$. Moreover, the annulus is partitioned into N sectors equally, in each of which a control point is placed. These points are sorted with respect to the azimuth and denoted as $p_i\in {\mathbb{R}}^2$ for $i=1,\cdots ,N$ in Cartesian coordinates. Each pair of adjacent points is augmented by two other points and then connected via a cubic Bézier curve. For each control point, $p_i$, an auxiliary angle, is calculated to determine the tangent to the curve at this point, which is given by

$${\theta }_i^{*}={\alpha }_i{\theta }_{i-1,i}+(1-{\alpha }_i){\theta }_{i,i+1}$$

(1)

where ${\theta }_{i,i+1}$ is the angle between point $p_i$ and $p_{i+1}$, and ${\alpha }_i\in [0,1]$ is an averaging parameter to modify the local smoothness of the curve. Then, two augmented points for the curve between $p_i$ and $p_{i+1}$ are calculated by

$$\begin{array}{cc}\hfill p_i^{*}& =p_i+{\eta }_i\parallel p_{i+1}-p_i\parallel ·e_i\hfill \\ \hfill p_i^{**}& =p_{i+1}-{\eta }_i\parallel p_{i+1}-p_i\parallel ·e_{i+1}\hfill \end{array}$$

(2)

where $e_i={\left[\cos \left({\theta }_i^{*}\right),\sin \left({\theta }_i^{*}\right)\right]}^{\top }$, and the scale parameter ${\eta }_i$ controls the local curvature. The curve connecting $p_i$ and $p_{i+1}$ is given by

$$b\left(t\right)={(1-t)}^3{p}_i+3{(1-t)}^{2}t{p}_i^{*}+3(1-t)t^2{p}_i^{**}+t^3{p}_{i+1},0\le t\le 1$$

(3)

An example of a valid shape is illustrated in Figure 1. The morphing wings in this work take NACA-2424 [46,47] as the baseline airfoil.

2.2. Dynamic System of Airfoil Morphing

Denote the polar coordinate of each point $p_i$ as ${\left\{r_i,{\theta }_i\right\}}_{i=1}^N$. Then the shape of the airfoil is fully determined by the radius $r_i$, the angle ${\theta }_i$ and the auxiliary angle parameters ${\alpha }_i$, ${\eta }_i$. In the flight procedure of the aircraft, we aim to morph the airfoil to maximize desired aerodynamic performances, such as the lift-to-drag ratio $C_l/C_d$, according to various conditions, including flight position, velocity or attack angle. With the combination of CFD solvers and optimization methods, the preferred airfoil shape at a given flight condition can be determined previously. Such optimized shapes serve as the reference or target airfoil shapes for the morphing task. Then, the problem of how to control the airfoil to achieve optimized shapes during the flight procedure where the flight condition varies is investigated.

Since there has been a variety of investigations on the position control of DC motors [48], we assume that the polar angle of each control point tracks the reference trajectory well via motors. Additionally, we assume that the auxiliary angle parameters can be adjusted rapidly and accurately. Therefore, in this work, we focus on the optimal morphing in the aspect of modifying the radii of control points.

Smart materials, especially shape memory alloys (SMA), are adopted to realize the airfoil morphing, where the radii of control points are modified via adjusting the length of SMA wires. Firstly, the dynamic model is constructed between the wire temperature and radii. An SMA wire changes its length through the crystal phase transformation between martensite and austenite according to the temperature. The transitions to martensite and austenite have different start and end temperatures, which leads to the hysteresis properties of the strain with respect to the temperature. Instead of common methods such as Preisach model and Krasnosel’skii–Pokrovskii model [49], the SMA hysteresis is characterized using hyperbolic tangent functions for their efficiency in computation and accuracy in curve fitting [41]. The strain is replaced by a radius factor ${\gamma }_i=(r_i-R_1)/(R_2-R_1)$ equivalently such that ${\gamma }_i\in \left[0,1\right]$. For heating and cooling starting with temperatures outside the transformation region, namely that the initial temperature is not between the end temperatures of the phase transformations, the hysteresis properties are modeled by the major hysteresis loops as

$$f_l^{major}\left(T\right)=\frac{h_0}{2}\tanh \left(\left(T-c_{tl}\right)c_b\right)+w\left(T-\frac{c_{tl}+c_{tr}}{2}\right)+\frac{h_0}{2}+c_s$$

(4)

$$f_r^{major}\left(T\right)=\frac{h_0}{2}\tanh \left(\left(T-c_{tr}\right)c_b\right)+w\left(T-\frac{c_{tl}+c_{tr}}{2}\right)+\frac{h_0}{2}+c_s$$

(5)

where T denotes the temperature, and $h_0$, $c_{tl}$, $c_{tr}$, $c_b$, w, $c_s$ parameterize the shape of the curves. Values of such parameters are chosen to fit the experimental data of SMA [40,41]. The radius factor $\gamma$ varies according to the $f_l^{major}$ curve when the temperature decreases and $f_r^{major}$ when the temperature increases. Furthermore, switching the temperature direction during the transformation procedure causes a reverse transformation starting from the current temperature and strain, which is not on the major loop of the reverse transformation. Such transforms are modeled using minor hysteresis loops, which are modeled by hyperbolic tangent curves with similar shapes as the major loops. The function of a rising minor loop is given as

$$f_r^{minor}(T,h)=\frac{h}{2}\tanh \left(\left(T-c_{tr}\right)c_b\right)+w\left(T-\frac{c_{tl}+c_{tr}}{2}\right)+h_0-\frac{h}{2}+c_s$$

(6)

where h is selected to ensure the intersection of the consecutive curves at the current point and is given by

$$h=g_r(h_{\mathrm{prev}},T)=\frac{h_{\mathrm{prev}}\left(\tanh \left(\left(T-c_{tl}\right)c_b\right)+1\right)-2h_0}{\tanh \left(\left(T-c_{tr}\right)c_b\right)-1}$$

(7)

and $h_{\mathrm{prev}}$ is the height parameter of the previous curve. Functions for the lowering curves are analogous as

$$f_l^{minor}(T,h)=\frac{h}{2}\tanh \left(\left(T-c_{tl}\right)c_b\right)+w\left(T-\frac{c_{tl}+c_{tr}}{2}\right)+\frac{h}{2}+c_s$$

(8)

and

$$h=g_l(h_{\mathrm{prev}},T)=\frac{h_{\mathrm{prev}}\left(\tanh \left(\left(T-c_{tr}\right)c_b\right)-1\right)+2h_0}{\tanh \left(\left(T-c_{tl}\right)c_b\right)+1}$$

(9)

An illustration of the transformation procedure is given in Figure 2.

After constructing the temperature-strain model, resist heating is used to actuate the SMA wires [50]. Given the applied voltage $v_i$, the temperature T follows the heat transfer model [39]

$$m_w{c}_w\dot{T_i}=\frac{v_i^2}{R_w}-h_w{A}_w\left(T_i-T_f\right) i=1,\cdots ,N$$

(10)

where $m_w$ is the mass per unit length of the SMA wire, $c_w$ is the specific heat, $R_w$ is the electrical resistance per unit length, $h_w$ is the heat exchange coefficient, $A_w$ is the wire circumferential area, and $T_f$ is the airflow temperature. Combining the temperature-strain relationship and (10), it is shown that the dynamic system between the radius and the input voltage are highly nonlinear due to the hysteresis characteristics. An illustration on the dynamics of SMA wires driven by a sinusoidal voltage input is given in Figure 3.

In this work, we tackle the morphing problem for the airfoil constructed by the SMA wires whose dynamics are given in (4)–(10). Note that the temperature-strain and voltage-temperature relationships modeled above are not directly accessible to our controller, but serve as the environment from which paths of the states can be sampled. We resort to deep reinforcement learning methods to design the morphing policy.

2.3. Reinforcement Learning based Morphing Control

Reinforcement learning (RL) methods are capable of learning a control policy from interactions between the given agent and environment, with no requirement on the knowledge of system models [51]. The learning procedures are based on Markov decision processes (MDPs), which are given by 4-tuples $\left\{S,A,R,P\right\}$, where S and A are the state space and action space containing all the states and actions, respectively, R is the reward function giving $r_k=R(s_k,a_k,s_{k+1})$ as rewards, and P is the transition function giving $P\left(s_{k+1}|s_k,a_k\right)$ as state transition probabilities. In this section, the airfoil morphing problem is solved in the RL framework, where we aim to find the optimal policy maximizing the expected total rewards.

Before choosing the states and actions for the morphing problem, the morphing system is investigated further to construct an MDP from it. Firstly, the voltage-temperature dynamics (10) of SMA wires is discretized via Euler methods as

$$T_{i,k}=T_{i,k-1}+\Delta t·\sigma \left(T_{i,k-1},v_{i,k-1}\right)$$

(11)

where $\Delta t>0$ is the discretizing time step, and

$$\sigma \left(T,v\right)=\frac{v^2}{m_w{c}_w{R}_w}-\frac{h_w{A}_w}{m_w{c}_w}\left(T-T_f\right)$$

(12)

With regard to the temperature-strain dynamics, since the minor loops converge with the major loops outside the SMA’s transformation region, we assume that the major loops determine the initial states of the wires, and minor loops dominate the dynamics during the morphing procedure. We denote the function (8) as $f_l(T,h)$ and (6) as $f_r(T,h)$ for simplicity, and introduce the signum function

$$\mathrm{sgn}\left(x\right)=\left\{\begin{array}{c}1,x\ge 0\hfill \\ 0,x<0\hfill \end{array}\right.$$

(13)

Then, $\mathrm{sgn}(T_k-T_{k-1})$ is used to discriminate the status of raising or lowering the temperature at time k. According to (11) and the fact that $\Delta t>0$, we describe the temperature direction at time k by

$${\chi }_{i,k}\triangleq \mathrm{sgn}\left(T_k-T_{k-1}\right)=\mathrm{sgn}\left(\sigma \left(T_{i,k-1},v_{i,k-1}\right)\right)$$

(14)

Note that the strain is dependent on the height parameter of the current loop. Recall from (7) and (9) that the value of the height parameter changes when the direction of temperature switches. Then, the time-varying parameter is determined by

$$\begin{array}{cc}\hfill h_{i,k}=& \left(1-{\chi }_{i,k-1:k}\right)·\left({\chi }_{i,k}{g}_r\left(T_{i,k},h_{i,k-1}\right)+\left(1-{\chi }_{i,k}\right)g_l\left(T_{i,k},h_{i,k-1}\right)\right)\hfill \\ +{\chi }_{i,k-1:k}{h}_{i,k-1}\hfill \end{array}$$

(15)

where

$${\chi }_{i,k-1:k}\triangleq \mathrm{sgn}\left(\sigma \left(T_{i,k-2},v_{i,k-2}\right)\sigma \left(T_{i,k-1},v_{i,k-1}\right)\right)$$

(16)

detects the reversal of temperature direction. Subsequently, the radius factor given the temperature and height parameter is calculated by choosing the rising or lowering loop according to the current temperature direction as

$${\gamma }_{i,k}={\chi }_{i,k}{f}_r\left(T_{i,k},h_{i,k}\right)+\left(1-{\chi }_{i,k}\right)f_l\left(T_{i,k},h_{i,k}\right)$$

(17)

Summarizing the relationships (11), (15) and (17), it seems appropriate to select T and h as states and v as actions to construct a RL environment for airfoil morphing. The radius factor can be treated as an observation since it is dependent on only the current voltage, temperature and height parameter. Since the length and temperature of SMA wires can be measured directly via sensors equipped on the airfoils, these values are assumed to be available. However, the height parameter is not a realistic physical characteristic but just a coefficient fitting the hyperbolic tangent curves to the actual SMA properties, which makes the measurement on h not available. Actually, we aim to adjust the position of control points to achieve a reference airfoil shape, such that we want to learn a policy on modifying the lengths (i.e., the radius factors). We observe from the loop functions (6) and (8) that, given the temperature T, the function of $\gamma$ with respect to h is bijective. Therefore, denoting the relationship (17) as ${\gamma }_{i,k}=f\left(h_{i,k},T_{i,k},T_{i,k-1},v_{i,k-1}\right)$, we can obtain an inverse function as

$$h_{i,k}=f^{-1}\left({\gamma }_{i,k},T_{i,k},T_{i,k-1},v_{i,k-1}\right)$$

(18)

Combining (15) and (18), the dynamics of radius factors are described as

$$\begin{array}{cc}\hfill {\gamma }_{i,k}=& f((1-{\chi }_{i,k-1:k})·[{\chi }_{i,k}{g}_r\left(T_{i,k},f^{-1}({\gamma }_{i,k-1},T_{i,k-1},T_{i,k-2},v_{i,k-2})\right)\hfill \\ & +\left(1-{\chi }_{i,k}\right)g_l\left(T_{i,k},f^{-1}({\gamma }_{i,k-1},T_{i,k-1},T_{i,k-2},v_{i,k-2})\right)]\hfill \\ & +{\chi }_{i,k-1:k}{f}^{-1}({\gamma }_{i,k-1},T_{i,k-1},T_{i,k-2},v_{i,k-2}),T_{i,k},T_{i,k-1},v_{i,k-1})\hfill \end{array}$$

(19)

Note that ${\chi }_{i,k-1:k}$ is dependent on the temperature and voltage at time $k-2$. This makes the transition function of $\gamma$ a second-order difference equation.

Proposition 1.

Given a second-order MDP with states $s^{\prime }\in {\mathbb{R}}^{n_{s^{\prime }}}$ and actions $a^{\prime }\in {\mathbb{R}}^{n_{a^{\prime }}}$ satisfying

$$p(s_k^{\prime }|s_{k-1}^{\prime },a_{k-1}^{\prime },s_{k-2}^{\prime },a_{k-2}^{\prime },\dots ,s_0,a_0)=p(s_k^{\prime }|s_{k-1}^{\prime },a_{k-1}^{\prime },s_{k-2}^{\prime },a_{k-2}^{\prime }) k\ge 2$$

(20)

select $s_k={\left[s_k^{\prime \top },s_{k-1}^{\prime \top }\right]}^{\top }$ and $a_k={\left[a_k^{\prime \top },a_{k-1}^{\prime \top }\right]}^{\top }$. Then, an MDP can be constructed by states s and actions a (with an additionally defined state transition function and reward function). Moreover, if the second-order MDP with $s^{\prime }$ and $a^{\prime }$ has a deterministic state transition function

$$s_k^{\prime }=f_s(s_{k-1}^{\prime },a_{k-1}^{\prime },s_{k-2}^{\prime },a_{k-2}^{\prime })$$

(21)

then the MDP constructed by s and a satisfies the transition function

$$s_k={\left[f_s^{\top }\left(\left[{\mathbf{I}}_{n_{s^{\prime }}},0\right]s_{k-1},\left[{\mathbf{I}}_{n_{a^{\prime }}},0\right]a_{k-1},\left[0,{\mathbf{I}}_{n_{s^{\prime }}}\right]s_{k-1},\left[0,{\mathbf{I}}_{n_{a^{\prime }}}\right]a_{k-1}\right),{\left(\left[0,{\mathbf{I}}_{n_{s^{\prime }}}\right]s_{k-1}\right)}^{\top }\right]}^{\top }$$

(22)

This proposition can be derived directly via the properties of MDPs. According to Equations (11), (19) and Proposition 1, it is reasonable to choose the radius factors $\left\{{\gamma }_{i,k-1},{\gamma }_{i,k}\right\}$ and temperatures $\left\{T_{i,k-1},T_{i,k}\right\}$ as states with input voltages $\left\{v_{i,k-1},v_{i,k}\right\}$ as actions. We restrict the states by $\gamma \in \left[0,1\right]$ and the actions by $v\in \left[0,V_{\max }\right]$.

The reference radius factors of the optimized airfoil shape under certain condition $c_k$ are denoted as ${\gamma }_{i,k}^{\mathrm{ref}}(c_k)$. When determining the reward, we expect the airfoil to morph to reference shapes accurately and rapidly. Therefore, a sparse reward function comprised of two components at time k is designed by

$$r_k^{\prime }=R\left(s_k,a_k,c_k\right)=\sum _{i=1}^N(r_{k,i}^{\mathrm{pos}}+r_{k,i}^{\mathrm{vol}})$$

(23)

where the position reward conveying the requirement of morphing accuracy is given by

$$r_{k,i}^{\mathrm{pos}}=\left\{\begin{array}{c}{r}_{\mathrm{p}},\hspace{1em}|{\gamma }_{i,k}-{\gamma }_{i,k}^{\mathrm{ref}}(c_k)|\le e_{\mathrm{thr}}\hfill \\ 0,\hspace{1em}|{\gamma }_{i,k}-{\gamma }_{i,k}^{\mathrm{ref}}(c_k)|>e_{\mathrm{thr}}\hfill \end{array}\right.$$

(24)

and the voltage reward aiming to increase the morphing speed is given by

$$r_{k,i}^{\mathrm{vol}}=\left\{\begin{array}{c}{r}_{\mathrm{v}},\hspace{1em}|{\gamma }_{i,k}-{\gamma }_{i,k}^{\mathrm{ref}}(c_k)|>e_{\mathrm{thr}} \mathrm{and} \mathrm{sgn}\left(v_{i,k}-v_{i,k-1}\right)=\mathrm{sgn}\left({\gamma }_{i,k}^{\mathrm{ref}}(c_k)-{\gamma }_{i,k}\right)\hfill \\ 0,\hspace{1em}\mathrm{otherwise}\hfill \end{array}\right.$$

(25)

where $r_{\mathrm{p}}>r_{\mathrm{v}}>0$ and $e_{\mathrm{thr}}>0$ are tunable hyperparameters. When the position error is small, the voltage reward is eliminated to mitigate the oscillation. The choice of adjacent voltages as actions permits the calculation of the voltage reward in the MDP framework. Furthermore, the total return to be maximized is given by

$$R\left(\tau \right)=\sum _{k=0}^K{r}_k^{\prime }$$

(26)

where $\tau =\left(s_0,a_0,c_0,s_1,a_1,c_1\dots \right)$ denotes the sequence of states, actions and conditions.

After establishing the MDP, we proceed to tackle the morphing task based on deep reinforcement learning techniques. Our learning method is designed based on the soft actor-critic (SAC) algorithm [52], which is an off-policy reinforcement learning method compatible with continuous state and action spaces. In the actor-critic framework, the agent learns to interact with the environment and obtain maximum rewards via training two types of neural networks iteratively. The first one is named a critic network and accepts current states and actions as input to approximate the action-value function, which serves as an evaluation of the current policy. The second one is denoted as an actor network and generates actions according to the system states and optional external inputs. After the training is converged, the policy is determined by the actor network and conducted for online executions, which in this work is the morphing task. In SAC, a stochastic policy is learned with additional entropy regularization in the rewards, which improves the ability of exploration and achieves faster convergence for a variety of control problems. According to the MDP of the morphing procedure, we choose the state $s_k$ to be $s_k=\left[s_{1,k-1}^{\prime },s_{1,k}^{\prime },\dots ,s_{N,k-1}^{\prime },s_{N,k}^{\prime }\right]$ and the action as $a_k=\left[a_{1,k-1}^{\prime },a_{1,k}^{\prime },\dots ,a_{N,k-1}^{\prime },a_{N,k}^{\prime }\right]$, where $s_{i,k}^{\prime }=\left[{\gamma }_{i,k},T_{i,k}\right]$ and $a_{i,k}^{\prime }=v_{i,k}$. Note that since the reference airfoil shapes guide the morphing, the flight condition should be incorporated for the generation and evaluation of the actions. We denote the distribution of the stochastic policy as $\pi \left(a\right|s,c)$.

With regard to the construction of the action-value function, instead of directly applying (23), we augment the reward with policy entropy as

$$r_k=R(s_k,c_k,a_k)=\sum _{i=1}^N\left(r_{k,i}^{\prime }+\beta H_{\pi }(s_{i,k},c_k)\right)$$

(27)

where

$$H_{\pi }(s_{i,k},c_k)={\mathbb{E}}_{a_{i,k}\sim \pi (·|s_{i,k},c_k)}\left[-\log \pi \left(a_{i,k}\right|s_{i,k},c_k)\right]$$

(28)

is the entropy representing the randomness of the policy, and $\beta >0$ is the trade-off coefficient. According to (26), we define a finite-horizon undiscounted return to be maximized. Nevertheless, a discount factor is applied when evaluating the value functions to focus on recent rewards, since the future reference shapes are not accessible at the current time step. Then, the action-value function with the regularized reward function is introduced as

$$Q^{\pi }(s,a,c,k)={\mathbb{E}}_{\tau \sim \pi }\left[-\sum _{t=k}^K\sum _{i=1}^N{\xi }^t{r}_{k,i}^{\prime }+\sum _{t=k+1}^K\sum _{i=1}^N{\xi }^t\beta H_{\pi }(s_{i,t},c_t)|s_k=s,a_k=a,c_k=c\right]$$

(29)

Afterwards, the Bellman equation for the action-value function is given as

$$\begin{array}{cc}\hfill Q^{\pi }& (s_k,a_k,c_k,k)=\hfill \\ & r_k+{\mathbb{E}}_{s_{k+1}\sim P,a_{k+1}\sim \pi }\left[\xi \left(Q^{\pi }(s_{k+1},a_{k+1},c_{k+1},k+1)-\beta \log \pi \left(a_{k+1}\right|s_{k+1},c_{k+1})\right)\right]\hfill \end{array}$$

(30)

and following the schedule of SAC [52], we approximate the expectation by

$$Q^{\pi }(s_k,a_k,c_k,k)\approx r_k+\xi \left(Q^{\pi }\left(s_{k+1},{\tilde{a}}_{k+1},c_{k+1},k+1\right)-\beta \log \pi \left({\tilde{a}}_{k+1}|s_{k+1},c_{k+1}\right)\right)$$

(31)

where $\left\{s_k,a_k,c_k,s_{k+1},c_{k+1}\right\}$ are sampled from replay buffers, and the next action ${\tilde{a}}_{k+1}$ is sampled from the current policy $\pi (·|s_{k+1},c_{k+1})$.

For the learning of the action-value functions, the double-Q trick is applied to avoid overestimation [53]. Two critic networks for Q functions are implemented as $Q_{{\varphi }_1}(s,a,c,k)$ and $Q_{{\varphi }_2}(s,a,c,k)$, where ${\varphi }_1$ and ${\varphi }_2$ are parameters. Additionally, for the stabilization of the training procedure, the target networks $Q_{{\varphi }_1^{\prime }}$ and $Q_{{\varphi }_2^{\prime }}$, which are copies of $Q_{{\varphi }_1}$ and $Q_{{\varphi }_2}$, are used and updated by polyak averaging after each time we update the main critic networks as

$${\varphi }^{\prime }\leftarrow \rho {\varphi }^{\prime }+(1-\rho )\varphi$$

(32)

where $\rho \in (0,1)$ is the update hyperparameter. Summarizing all these settings, the loss for critic networks is given by the mean squared Bellman error function as

$$L\left(\varphi \right)=\sum _{b\in B}\left(Q_{\varphi }(s_k,a_k,c_k,k)-y(s_{k+1},a_{k+1},c_{k+1},k+1)\right)$$

(33)

where B is the sampled batch with elements $b=\left\{s_k,a_k,c_k,r_k,k,s_{k+1},c_{k+1},k+1\right\}$ of the replay buffer, and where

$$y(s_{k+1},a_{k+1},c_{k+1},k+1)=\xi \left(Q_{{\varphi }^{\prime }}\left(s_{k+1},{\tilde{a}}_{k+1},c_{k+1},k+1\right)-\beta \log \pi \left({\tilde{a}}_{k+1}|s_{k+1},c_{k+1}\right)\right)$$

(34)

is calculated using target networks. Stochastic gradient descent is applied to update ${\varphi }_1$ and ${\varphi }_2$ with respect to the loss functions $L\left({\varphi }_1\right)$ and $L\left({\varphi }_2\right)$.

Subsequently, we aim to find the policy that maximizes the expected action-value function with respect to the actions. Denote the parameters of actor network as $\theta$. Since the action at time step k is composed of the input voltages at $k-1$ and k, the actor network should be designed as ${\pi }_{\theta }(s_k,c_k,a_{k-1})$, where an identity layer is applied to propagate the previous voltages. However, the distribution of new input voltages is still only dependent on $s_k$ and $c_k$. Therefore, we use ${\pi }_{\theta }\left(a_k\right|s_k,c_k)$ to denote the density of action $a_k$. Then, the value function to be optimized is given as

$$V^{{\pi }_{\theta }}(s_k,c_k,k)={\mathbb{E}}_{a_k\sim {\pi }_{\theta }}\left[Q^{{\pi }_{\theta }}(s_k,a_k,c_k,k)-\beta \log \left({\pi }_{\theta }\left(a_k\right|s_k,c_k)\right)\right]$$

(35)

The reparameterization trick is adopted here for the sake of the efficient computation of gradients [54]. We introduce a standard normal distributed variable $\zeta \sim p\left(\zeta \right)=\mathcal{N}\left(0,{\mathbf{I}}_{n_a}\right)$ and calculate the input voltage according to a deterministic squashing function as

$$v_k=\frac{V_{\max }}{2}\left(1+\tanh \left({\mu }_{\theta }(s_k,c_k)+{\sigma }_{\theta }(s_k,c_k)\circ {\zeta }_k\right)\right)$$

(36)

where ${\mu }_{\theta }(s,c)$ and ${\sigma }_{\theta }(s,c)$ are parameterized neural networks and ∘ denotes element-wise multiplication. Then, the expectation over $a_k\sim \pi$ can be converted to the expectation over the normal variable whose distribution is irrelevant to the states and net parameters. Additionally, the squashed Gaussian policy constrains the input voltages in $\left[0,V_{\max }\right]$. According to (36), the action given $s_k$, $c_k$ and ${\zeta }_k$ is written as $a_k=a_{\theta }(s_k,c_k,{\zeta }_k)$. Note that this is an invertible map between $a_k$ and ${\zeta }_k$. Therefore, we can compute the log-probabilities in closed form according to the change of the variable formula [55] as

$$\begin{array}{cc}\hfill \log \left({\pi }_{\theta }\left(a_k\right|s_k,c_k)\right)& =\log \left(p\left({\zeta }_k\right)\right)-\log \left|\det J_{a_{\theta }}\left({\zeta }_k\right)\right|\hfill \\ & =\log \left(p\left(a_{\theta }^{-1}(s_k,c_k,a_k)\right)\right)+\log \left|\det J_{a_{\theta }^{-1}}\left(a_k\right)\right|\hfill \end{array}$$

(37)

where $a_{\theta }^{-1}$ is the inverse function of $a_{\theta }$ given $s_k$ and $c_k$, and J denotes the Jacobian matrix.

With $Q^{{\pi }_{\theta }}$ approximated by the minimum of the two critic networks, the loss function for the actor network is obtained as

$$L\left(\theta \right)=\sum _{b\in B}\left(\underset{i=1,2}{\min }{Q}_{{\varphi }_i}(s_k,a_{\theta }(s_k,c_k,{\zeta }_k),c_k,k)-\beta \log ({\pi }_{\theta }\left(a_{\theta }(s_k,c_k,{\zeta }_k)\right|s_k,c_k)\right)$$

(38)

Then, we conduct the training by updating the actor and critic networks iteratively. The agent interacts with various randomly generated time-varying reference shape sequences to acquire training data, which faciliates the exploration and enables the policy to handle different morphing scenarios. When the training converges, we can use the actor network to calculate the required voltages and morph the airfoil to the reference shapes. Finally, the overall flowchart of the proposed morphing mechanism is given in Figure 4.

3. Results

In this section, a simulation is conducted to validate the proposed morphing method. The simulation is arranged in two stages. Firstly, we implement our method with random generated reference shapes and perform ablation studies to examine the superiorities brought by different parts of our algorithm. Subsequently, we apply the proposed method to track optimized airfoil shapes in different flight conditions and show the morphing procedures. The values of parameters in our simulation are given in

Table 1. Values of system parameters used in simulations.

Parameter	Value	Parameter	Value
$m_w$	$1.14\times {10}^{-4}$	$A_w$	$4.72\times {10}^{-4}$
$c_w$	837.4	$R_w$	50.8
$T_f$	20	$h_w$	120
H	0.995	$c_b$	0.147
$c_w$	$1.25\times {10}^{-5}$	$c_s$	0.001
$c_{tl}$	46	$c_{tr}$	65

[40,41].

3.1. Tracking Random Shapes

In this stage, the superiorities of our method are illustrated in a variety of perspectives, including the state/action selection, reward configuration and entropy regularization. Without loss of generality, piece-wise constant trajectories of the radius factor for one control point are generated to represent the reference shapes. Then, our method is compared with three different settings of RL algorithms. The proposed method is denoted as RLM-SAC.

• Second-order state/action versus first-order state/action
  In Section 2, second-order MDP is adopted to model the hysteresis characteristics of the morphing system. Therefore, we chose the states and actions as combinations of that in current step and previous step. We compared the performance with RL algorithms where the policy is generated according to only current states, and the value function was also evaluated with only current states and actions as inputs that are applied in existing investigations on controling SMA wires. We refer to this as RLM-FO.
• Sparse reward versus squared error reward
  We designed a sparse reward taking value in $\left\{0,1\right\}$, which is different from traditional RL-based morphing research. We compared that with the square error rewards, which is given by

  $$r_k^{\prime }=R\left(s_k,c_k\right)=-\sum _{i=1}^N{|{\gamma }_{i,k}-{\gamma }_{i,k}^{\mathrm{ref}}\left(c_k\right)|}^2$$

  (39)

  which is named RLM-SER.
• SAC versus DQN
  The entropy regularization improves the capability of exploration in our algorithm. A modified deep Q learning method was implemented as a comparison, where only the entropy loss was removed, and both the double-Q setting and reparameterization trick remained. We denote this as RLM-DQN.

All RL realizations were trained through 150 epochs, in each of which 5 episodes with 40 s of time and 200 time steps were executed. The critic network was constructed using multilayer perceptrons (MLP) of 3 hidden layers and 128 units per layer. The actor network adopted similar structures, where additional fully connected layers were attached to produce the mean and standard variations of the policy. The training was started with actions uniformly sampled from the valid action space bounded by $\left[0,V_{\max }\right]$ for 5000 steps to explore the state space sufficiently. Then, the networks were updated every step with a batch size of 200. A fixed learning rate was set as $0.002$, and other hyperparameters used in RL training are shown in

Table 2. Values of hyperparameters used in RL.

Parameter	Value	Parameter	Value
$\xi$	0.98	$\alpha$	0.2
$\rho$	0.995	$V_{\max }$	15
$r_{\mathrm{p}}$	1	$r_{\mathrm{v}}$	0.02
$e_{\mathrm{thr}}$	0.1

. The actions were generated from the stochastic policy in training phase but produced in a deterministic way according to the mean value of the actor network in the testing phase.

The results are shown in Figure 5, Figure 6, Figure 7 and Figure 8. In Figure 5 and Figure 6, we illustrate the rewards acquired in the training and test trajectories, respectively, during the training procedure. A Savitzky–Golay filter [56] was adopted to smooth the data such that the values and trends of the rewards are illutstrated more clearly. Note that the reward of RLM-SER was not included because of a different reward setting. After the training was finished, the algorithms were executed on 100 random generated test reference trajectories. We present the root-mean-squared error (RMSE) of the radius factor through the flight time in Figure 7. Some of the trajectories and the corresponding results produced by the four RL realizations are depicted in Figure 8 for an intuitive comparison.

From the results, we can see the benefits of each important component in our algorithm. Firstly, RL using square error reward totally fails to produce effective actions in our environment, which is shown from both the RMSE and example trajectories. RLM-FO acquires inferior performance compared with RLM-SAC, especially in the temperature-switching procedure, which can be validated by the middle sections of the example trajectories. This is a result of the fact that the hysteresis cannot be characterized well by first-order states and actions. Lastly, RLM-DQN obtains better performance than RLM-FO and RLM-SRE, and it achieves a similar reward to RLM-SAC at the end of training. However, from the illustrations of rewards in Figure 5 and Figure 6, it is shown that RLM-DQN converges much slower than RLM-SAC. This is due to the improvement of exploration capability provided by entropy regularization.

3.2. Morphing Procedure Simulation

In this stage, the trained actor network is applied to morph an airfoil controlled by four points in a given flight procedure with varying flight conditions. Since the focus of this work is morphing control, in each condition, the optimal shape is assumed to be solved in anticipation by shape optimization techniques and determined by the radii and angles of control points [14,18,20]. The average coefficient and scale coefficient are fixed as $\alpha =[0.12,0.4,0.4,0.12]$ and $\eta =[0.5,0.5,0.5,0.5]$. The trajectory of flight conditions and the corresponding parameters of optimal shapes are shown in Figure 9.

Since the position control technique of motors is relatively sufficiently developed, we assumed that the angles of the points are controlled by DC motors with accurately known linear dynamic models, which can drive the points to desired angles rapidly with subtle errors. Then we generate input voltages to heat the SMA wires and adjust the radius factors. The voltages, temperatures and radius factors of all points are summarized in Figure 10. It is shown that with constrained voltages, the wires can track the reference lengths well.

Furthermore, we illustrate the morphing procedures in Figure 11. We can see intuitively that with the proposed RLM-SAC method, the airfoil is capable of morphing into the optimized shape within about 3 s after encountering a new flight environment, which validates both the morphing accuracy and morphing speed.

Quantitative comparisons on the length factor differences and shape differences are given in

Table 3. Quantitative comparisons of different morphing methods.

	RLM-SAC	RLM-FO	RLM-SER	RLM-DQN
$e_{\mathrm{length}}$	14.04%	31.44%	323.6%	23.24%
$e_{\mathrm{shape}-\mathrm{avg}}$	0.0533	0.0645	0.2202	0.0601
$e_{\mathrm{shape}-\mathrm{end}}$	0.0019	0.0055	0.0494	0.0031

. The average length factor differences are calculated according to the reference length factor and actual length factor as

$$e_{\mathrm{length}}=\frac{1}{KN}\sum _{k=0}^K\sum _{i=1}^N\frac{|{\gamma }_{i,k}-{\gamma }_{i,k}^{\mathrm{ref}}|}{{\gamma }_{i,k}^{\mathrm{ref}}}$$

(40)

The difference between the actual and reference shapes are evaluated using distances between the control points. An average shape difference over all time is calculated as

$$e_{\mathrm{shape}-\mathrm{avg}}=\frac{1}{KN}\sum _{k=0}^K\sum _{i=1}^N\left|\right|p_{i,k}-p_{i,k}^{\mathrm{ref}}{\left|\right|}_2$$

(41)

where $\left|\right|p_i-p_i^{\mathrm{ref}}{\left|\right|}_2$ is the $L_2$ distance between the two Cartesian coordinates $p_{i,k}$ and $p_{i,k}^{\mathrm{ref}}$. Additionally, the steady shape error is evaluated by

$$e_{\mathrm{shape}-\mathrm{end}}=\frac{1}{4N}\sum _{t=1}^4\sum _{i=1}^N\left|\right|p_{i,k_t}-p_{i,k_t}^{\mathrm{ref}}{\left|\right|}_2$$

(42)

where $k_t$ denotes the end time step of each flight condition. It is shown that our method acquires the best performance on all metrics. The proposed RLM-SAC method provides an average 29.8% performance improvement over the second-best RLM-DQN method. The length differences and average shape differences, which are averaged over all time steps, demonstrate that our method can morph the airfoil into desired shapes more rapidly, while the steady shape difference validates the morphing accuracy. Lastly, in this work, we give the reference airfoil shape directly and focus on the morphing performance of the proposed method. It will be interesting and meaningful to combine the shape optimization task with the morphing control problem in the future.

4. Conclusions

In this work, a novel deep reinforcement learning-based morphing control method is proposed for an asymmetric morphing airfoil. The airfoil is designed via Bézier curves and is capable of morphing from a baseline shape to an asymmetric shape. The morphing mechanism is modeled via SMA wires, which adjust shape parameters, especially the radii of the control points. To actuate the SMA wires, resistive heating is performed, but the hysteresis characteristics between the SMA strain and temperature make the dynamic system nonlinear and non-Markovian, which brings difficulties to the design of the control algorithm and the RL framework. Therefore, hyperbolic tangent curves are adopted to model the strain-temperature relationship and derive a second-order MDP describing the system, which is then transformed into a valid MDP and provides guidance for the selection of states and actions. Based on the constructed MDP, we modify the SAC algorithm and develop an RL scheme where input voltages are generated to morph the airfoil instantaneously according to reference-optimized shapes. Lastly, ablation studies on random generated reference trajectories are conducted to demonstrate the benefits brought by different components of our RL implementations, and we perform simulations of morphing procedures to validate that our method is able to morph the airfoil into the optimized shapes rapidly and accurately. Future works include incorporating the aerodynamic performance optimization directly into the morphing control and exploiting learning-based morphing policies for more complicated bio-inspired morphing aircraft.