Stochastic Drift Counteraction Optimal Control of a Fuel Cell-Powered Small Unmanned Aerial Vehicle

Abstract

This paper investigates optimal power management of a fuel cell hybrid small unmanned aerial vehicle (sUAV) from the perspective of endurance (time of flight) maximization in a stochastic environment. Stochastic drift counteraction optimal control is exploited to obtain an optimal policy for power management that coordinates the operation of the fuel cell and battery to maximize the expected flight time while accounting for the limits on the rate of change of fuel cell power output and the orientation dependence of fuel cell efficiency. The proposed power management strategy accounts for known statistics in transitions of propeller power and climb angle during the mission, but does not require the exact preview of their time histories. The optimal control policy is generated offline using value iterations implemented in Cython, demonstrating an order of magnitude speedup as compared to MATLAB. It is also shown that the value iterations can be further sped up using a discount factor, but at the cost of decreased performance. Simulation results for a $1.5$ kg sUAV are reported that illustrate the optimal coordination between the fuel cell and the battery during aircraft maneuvers, including a turnpike in the battery state of charge ($SOC$) trajectory. As the fuel cell is not able to support fast changes in power output, the optimal policy is shown to charge the battery to the turnpike value if starting from a low initial $SOC$ value. If starting from a high $SOC$ value, the battery energy is used till a turnpike value of the $SOC$ is reached with further discharge delayed to later in the flight. For the specific scenarios and simulated sUAV parameters considered, the results indicate the capability of up to $2.7$ h of flight time.

1. Introduction

With the growing market for unmanned aerial vehicles (UAVs), a wide range of industries and organizations, including military, government, industrial, and recreational users, deploy this technology across the globe [1,2,3]. Among different types of UAVs, small unmanned aerial vehicles (sUAVs) [4] are attractive for military, aerial photography, and environmental monitoring applications due to their small size and flexible operation [5]. Considering the (i) hardware and weight constraints, (ii) limited onboard energy storage, and (iii) performance requirements for sUAVs, improving their endurance (maximizing their flight time) is of great importance for extending the duration of their missions, which could involve surveillance, search and rescue, disaster relief, traffic control, and precision agriculture, thereby motivating the development of novel propulsion systems and the implementation of optimal control policies for power and energy management. Among different propulsion systems for such sUAVs, a hybrid propulsion system consisting of a polymer electrolyte membrane fuel cell (PEMFC) and a battery has been proposed for long duration missions, e.g., in [6,7,8,9]. Other propulsion systems may incorporate energy harvesters, such as in [10]. In this paper we focus on novel approaches to the energy management of sUAVs through optimal coordination between the PEMFC and battery for the previously proposed fuel cell hybrid propulsion system.

Rule-based (e.g., thermostat-like on-off control [11]), dynamic programming-based [12], and model predictive control (MPC) [13] methods have been considered for the energy management of hybrid aircraft. As in automotive energy management applications [14], the use of simple rule-based strategies may not provide optimal performance, while the conventional formulations of MPC and dynamic programming do not directly address the flight time maximization objective. Furthermore, deterministic variants of MPC and dynamic programming may require an accurate preview of the propeller power and climb angle over a long horizon and are computationally demanding if optimization has to be performed online. Similarly, the Pontryagin maximum principle (PMP)-based guidance solutions [15] need accurate characterization of the flight environment.

In this paper, we consider a different approach to the problem of endurance maximization for a hybrid UAV with a polymer electrolyte membrane fuel cell (PEMFC) based on an application of stochastic drift counteraction optimal control (SDCOC) [16], which directly addresses the problem of maximizing the time to constraint violation in a stochastic environment. In our case, the objective is to maintain the vehicle flying for a maximum amount of time by coordinating the fuel cell and the battery to provide the requested propeller power subject to the limited amount of fuel and battery state of charge ($SOC$) onboard the vehicle. The transitions in aircraft climb angle and propeller power are modeled stochastically by a Markov chain with the transition probabilities determined from historical data representing typical missions of an sUAV. Then, a control policy that minimizes a cost functional reflective of expected time-to-violate constraints is determined offline through value iterations; this control policy is then deployed onboard for the online coordination of the fuel cell and the battery in the sUAV.

In a preliminary conference paper [17] by the second author of this paper, the application of SDCOC for power management of a hybrid sUAV with a direct methanol fuel cell (DMFC) was considered. While the DMFC is often considered as a suitable power source for ground vehicles [18] and has certain advantages, PEMFCs are more appealing for air mobility applications [6,7] due to their relatively lower operating temperature, allowing for a quick start-up [19], higher efficiency (up to 60% [18,20]) and power density, and higher safety due to the use of the solid electrolyte [18].

Different from [17], in this paper, we consider the application of SDCOC to the power management of a hybrid sUAV with a PEMFC rather than a DMFC. To accommodate a different fuel cell and an sUAV, the fuel cell model is changed and improvements are made to the models used to compute the propeller power and thrust, as well as the evolution of the $SOC$.

More importantly, the lack of the ability of the PEMFC to rapidly change its power output imposes a stringent operating constraint (rate limit on PEMFC power output), which was not treated in [17], but is treated in this paper. This rate limit increases the complexity of the problem as an extra state needs to be introduced in the model and handled in SDCOC, and it also changes the optimal policies and the optimal response of the system. For instance, as the fuel cell is not able to support fast changes in power output, the optimal policy is shown to charge the battery to a turnpike value if starting from a low initial state of charge value. If starting from a high $SOC$, the battery energy is used till a turnpike value of the state of charge is reached with further discharge delayed to a later phase of the flight. In either case, the high frequency chattering of fuel cell load demand power in [17], which cannot be supported by the PEMFC, is eliminated.

Additionally, in this paper, the value iterations are implemented in Cython rather than MATLAB, with an order of magnitude speedup as compared to the MATLAB implementation. As value iterations are frequently used to solve dynamic programming problems in different applications and Python is becoming increasingly popular, our results on the ten-fold speedup with Cython without a substantive increase of the code complexity are of reference value to other researchers considering the computational implementation of dynamic programming.

Furthermore, a discount factor is introduced into the cost function of SDCOC, and its impact on the convergence speed of the value iterations is illustrated. It is shown that this discount factor results in the faster convergence of value iterations, but the performance of the control policy (in terms of exit time) is decreased.

While the SDCOC theory was developed in [16], that reference did not address the fuel cell or sUAV application studied in this paper. Our approach to representing motor power demand and climb angle by a Markov chain with a finite number of states follows [21], which is the first (to the authors’ knowledge) paper proposing the use of stochastic dynamic programming for automotive powertrain control applications; that paper also did not address the fuel cell or sUAV application studied in the present paper, nor the drift counteraction problem formulation.

The remainder of this paper is organized as follows: Section 2 describes the sUAV sub-systems and their models. Section 3 presents an integrated model of the hybrid system and defines the problem in a form suitable for SDCOC. Section 4 summarizes SDCOC, and Section 5 reports the results. Finally, Section 6 presents concluding remarks.

2. Physical Description of the Systems and Model

An sUAV with a series hybrid propulsion system, shown in Figure 1, was chosen in which the power supplied by the battery and the power supplied by the PEMFC are combined to meet the propeller motor power demand. The PEMFC uses hydrogen as the fuel, which is stored in the tank, and air from the atmosphere. A fraction of the energy generated by the PEMFC can be used to charge the battery. The fuel cell pack and battery pack are sized large enough so that they are able to meet the sUAV’s mean power demand individually, should either one not be operating properly.

The model used in this paper for generating the SDCOC policies captures the battery’s $SOC$ dynamics, the fuel cell’s hydrogen rate dynamics, and the fuel cell load power dynamics. Thus, the states of this model are the $SOC$, the mass of hydrogen remaining in the gas tank, and the fuel cell load demand power. The motor power of the sUAV and climb angle are treated as operating variables, and the SDCOC controller determines changes in the fuel cell load demand power. This system level model has been implemented by combining component submodels and characterizations available from the literature; our methodology is generic and can accommodate changes in these component models.

2.1. sUAV Dynamics

A control-oriented dynamic model of the sUAV is used for SDCOC law development. The sUAV is constrained to a longitudinal flight path in a vertical plane [22].

Table 1. List of variables used in the sUAV model.

Variable	Description	Unit
v	Velocity of the sUAV	$\mathrm{m}/\mathrm{s}$
$\gamma$	Climb angle	$\mathrm{deg}$
T	Thrust force	$\mathrm{N}$
$\alpha$	Angle of attack	$\mathrm{deg}$
L	Lift force	$\mathrm{N}$
D	Drag force	$\mathrm{N}$
$C_L$	sUAV coefficient of lift	∖
$C_D$	sUAV coefficient of drag	∖
${\rho }_{air}$	Air density	$\mathrm{kg}/{\mathrm{m}}^3$
$P_{sUAV}$	Power required by the sUAV	$\mathrm{W}$
N	Angular speed of the electric motor	$\mathrm{RPM}$
$P_P$	Power generated by the propeller	$\mathrm{W}$
$P_{P,ideal}$	Ideal propeller power	$\mathrm{W}$
${\eta }_P$	Propulsive efficiency	∖
$U_M$	Electric motor driver’s input voltage	$\mathrm{V}$
$I_M$	Electric motor driver’s input current	$\mathrm{A}$
$P_M$	Elector motor driver’s input power	$\mathrm{W}$
${\eta }_M$	Motor efficiency	∖
$P_{FC,total}$	Total power of the fuel cell	$\mathrm{W}$
$P_{FC,load}$	Load demand power of the fuel cell	$\mathrm{W}$
$P_{aux}$	Power required by the auxiliaries	$\mathrm{W}$
$U_{FC}$	Single cell voltage	$\mathrm{V}$
$I_{FC}$	Single cell current	$\mathrm{A}$
$i_{FC}$	Single cell current density	$\mathrm{A}/{\mathrm{cm}}^2$
$U_{act}$	Activation polarization	$\mathrm{V}$
$U_{ohm}$	Ohmic losses	$\mathrm{V}$
$U_{conc}$	Concentration polarization	$\mathrm{V}$
$U_{OC}$	Equivalent open circuit voltage of a single fuel cell	$\mathrm{V}$
$R_{FC}^{\prime }$	Modified single fuel cell resistance	$\mathsf{\Omega }$
${\tilde{R}}_{FC}^{\prime }$	Variable defined in Equation (18)	$\mathsf{\Omega }·{\mathrm{cm}}^2$
$U_{B,OC}$	Open circuit voltage of the battery	$\mathrm{V}$
$SOC$	Battery’s state of charge	∖
$P_B$	Power of the battery	$\mathrm{W}$
$SO{C}_0$	Initial SOC	∖
$S_f$	Split fraction	∖
u	Control input	∖
$\mathsf{\Delta }{P}_{FC}$	Defined in Equation (25)	$\mathrm{W}$
$m_{FR}$	Mass of fuel remaining	$\mathrm{kg}$

defines the notations for the variables used in the model.

Table A1. Parameters used in the sUAV model.

Variable	Description	Value	Unit
m	Mass of the sUAV	1.5	$\mathrm{kg}$
g	Gravitational acceleration	9.81	$\mathrm{m}/{\mathrm{s}}^2$
$S_{ref}$	Wing area	0.09	${\mathrm{m}}^2$
$C_{D0}$	sUAV coefficient of drag at $\alpha =0$	0.1038	∖
K	Coefficient in Equation (3)	0.0637	∖
$d_P$	Diameter of the propeller	0.24	$\mathrm{m}$
J	Advance ratio	0.37	∖
$R_M$	Motor resistance	0.105	$\mathsf{\Omega }$
$I_{M,0}$	Motor current at zero load	1.3	$\mathrm{A}$
$K_V$	Motor speed constant	1490	$\mathrm{RPM}/\mathrm{V}$
$n_{FC}$	Number of single cells in series	12	∖
$P_0$	Bias power of the fuel cell	5	$\mathrm{W}$
${\kappa }_{aux}$	Coefficient in Equation (12)	0.05	$\mathrm{V}$
$A_{FC}$	Fuel cell area	200	${\mathrm{cm}}^2$
$\mathsf{\Delta }{g}_f$	Molar specific Gibbs free energy	237.3	$\mathrm{kJ}/\mathrm{mol}$
$n_e$	Number of ions passed in reaction	2	∖
F	Faraday constant	96,485	$\mathrm{C}/\mathrm{mol}$
$\overline{T}$	Temperature of the reaction	333.15	$\mathrm{K}$
${\alpha }_{FC}$	Charge transfer coefficient	0.5	∖
R	Universal gas constant	8.314	$\mathrm{J}/(\mathrm{mol}·\mathrm{K})$
$R_{FC}$	Ohmic resistance defined in Equation (15)	0.0024	$\mathsf{\Omega }$
$d_0$	Coefficient in Equation (16)	3e-5	$\mathrm{V}$
$d_1$	Coefficient in Equation (16)	8	${\mathrm{cm}}^2/\mathrm{A}$
$k_0$	Coefficient in Equation (18)	4	∖
$k_1$	Coefficient in Equation (18)	1	∖
$M_h$	Molecular weight of $H_2$	2	$\mathrm{g}/\mathrm{mol}$
$n_B$	Number of batteries in series	8	∖
$U_{B,min}$	Open circuit voltage when $SOC=0$	2.5	$\mathrm{V}$
$U_{B,max}$	Open circuit voltage when $SOC=1$	4.2	$\mathrm{V}$
$R_{B,int}$	Battery internal resistance	0.012	$\mathsf{\Omega }$
$C_B$	Standard discharge capacity	14400	$\mathrm{C}$
$I_{B,max}$	Maximum discharge current	35	$\mathrm{A}$

in Appendix A lists the model parameter values, partly based on [23,24,25]. The development of lightweight electric components (batteries, fuel cells, motors) for sUAVs is an active area of research; see, e.g., [26,27]. In our model, we assumed that such lightweight components are available to be consistent with the assumed sUAV’s weight.

Using a flat Earth coordinate system, the longitudinal equations of motion of the sUAV are given by:

$$\dot{v}=\frac{Tcos\left(\alpha \right)-D}{m}-gsin\left(\gamma \right),$$

(1)

$$\dot{\gamma }=\frac{Tsin\left(\alpha \right)+L}{mv}-\frac{gcos\left(\gamma \right)}{v},$$

(2)

where v is the velocity of the sUAV and $\gamma$ is the climb angle. The lift L and drag force D are characterized as:

$$L=\frac{1}{2}{\rho }_{air}{v}^2{S}_{ref}{C}_L,D=\frac{1}{2}{\rho }_{air}{v}^2{S}_{ref}{C}_D,$$

(3)

where $C_L=C_{L0}+C_{L,\alpha }\alpha ,C_D=C_{D0}+K{C_L}^2$. Neglecting vertical acceleration (i.e., with $L=mg$), solving Equations (1) and (2) yields the thrust required by the sUAV,

$$\begin{array}{cc}\hfill T=& {\left[{\left(m\dot{v}+mgsin\left(\gamma \right)+\frac{1}{2}{\rho }_{air}{v}^2{S}_{ref}{C}_{D0}+\frac{2K{m}^2{g}^2}{{\rho }_{air}{v}^2{S}_{ref}}\right)}^2+{(mv\dot{\gamma }+mgcos\left(\gamma \right)-mg)}^2\right]}^{\frac{1}{2}}.\hfill \end{array}$$

(4)

Here, ${\rho }_{air}$ is a function of altitude. The power required by the sUAV is then given by:

$$P_{sUAV}=Tv.$$

(5)

2.2. Propeller Model

The propeller model is used to relate the torque and angular velocity generated by the electric motor to the power required by the sUAV and the velocity of the sUAV, respectively [22]. With the propulsive efficiency given by ${\eta }_P$, the power required to drive the propeller is:

$$P_P=\frac{P_{sUAV}}{{\eta }_P}.$$

(6)

According to the disk actuator theory, the ideal propeller power is:

$$P_{P,ideal}=\frac{1}{2}Tv\left(1+\sqrt{1+\frac{8T}{\pi {\rho }_{air}{v}^2{d_P}^2}}\right).$$

In general, the actual power required would be about 15% greater than this [28], which means $P_P=1.15P_{P,ideal}$. Thus, ${\eta }_P$ can be calculated as:

$${\eta }_P=\frac{P_{sUAV}}{1.15P_{P,ideal}}=\frac{2}{1.15+1.15\sqrt{1+\frac{8T}{\pi {\rho }_{air}{v}^2{d_P}^2}}}.$$

(7)

Combing Equation (7) with (5) and (6) yields:

$$P_P=\frac{1.15P_{sUAV}}{2}+\frac{1.15P_{sUAV}}{2}\sqrt{1+\frac{8P_{sUAV}}{\pi {\rho }_{air}{v}^3{d_P}^2}}.$$

2.3. Electric Motor Model

Electric motors used in sUAV applications exhibit high speed and high torque, as well as high power-to-weight ratios [29]. Assuming the power factor is equal to unity and the magnetic losses can be neglected, the output power of the motor is given by:

$$P_P=(U_M-I_M{R}_M)(I_M-I_{M,0}).$$

(8)

The angular velocity of the motor in revolutions per minute (RPM) can be expressed as:

$$N=(U_M-R_M{I}_M)K_V,$$

(9)

which should be equal to the RPM of propeller $N=\frac{v}{J{d}_P}$. From Equations (8) and (9), the motor current, $I_M$, is:

$$I_M=\frac{P_P{K}_V}{N}+I_{M,0}.$$

The motor power and motor efficiency are given by, respectively,

$$P_M=U_M{I}_M,{\eta }_M=\frac{P_P}{P_M}.$$

2.4. Fuel Cell Model

A PEMFC system is the primary power source for the sUAV. The total power generated by the fuel cell stack is calculated as:

$$P_{FC,total}=n_{FC}{U}_{FC}{I}_{FC}.$$

(10)

This power must cover the load demand $P_{FC,load}$ and the power required for auxiliaries [18], $P_{aux}$,

$$P_{FC,total}=P_{FC,load}+P_{aux},$$

(11)

where $P_{aux}$ is the total power required for the compressor motor, the hydrogen circulation pump, the humidifier water circulation pump, the coolant pump, the cooling fan motor, and the bias power, $P_0$. After simplifications, $P_{aux}$ could be written as a function of the fuel cell current [30],

$$P_{aux}=P_0+n_{FC}{\kappa }_{FC}{I}_{FC}.$$

(12)

The fuel cell current is a function of the current density and the fuel cell area,

$$I_{FC}=i_{FC}{A}_{FC},$$

where $i_{FC}$ could be obtained by solving the equation,

$$U_{FC}=U_{rev}-U_{act}-U_{ohm}-U_{conc}.$$

(13)

The reversible cell potential $U_{rev}$ is related to the molar specific Gibbs free energy $\mathsf{\Delta }{g}_f$ and the number of ions passed in the reaction $n_e$ [24],

$$U_{rev}=\frac{\mathsf{\Delta }{g}_f}{n_{e}F}.$$

The activation polarization $U_{act}$ is a result of the energy required to initiate the reaction, which can be described by the semi-empirical Tafel equation [31,32,33],

$$U_{act}=c_0+c_{1}ln\left(i_{FC}\right),$$

where $c_0$ and $c_1$ depend on temperature. When the current density is small, this equation can be modified [34] as:

$$U_{act}=c_0(1-e^{-c_1{i}_{FC}}),$$

(14)

where $c_0=-5.8\times {10}^{-4}\overline{T}+0.5736$ and $c_1=\frac{R\overline{T}}{n_e{\alpha }_{FC}F}$.

The ohmic losses $U_{ohm}$ are due to the resistance to the flow of (i) ions in the membrane and in the catalyst layers and (ii) electrons through the electrodes [18],

$$U_{ohm}=i_{FC}{\tilde{R}}_{FC},$$

(15)

where ${\tilde{R}}_{FC}=R_{FC}{A}_{FC}$.

The concentration polarization $U_{conc}$, is given by

$$U_{conc}=d_0{e}^{d_1{i}_{FC}}.$$

(16)

With the parameters given in Appendix A, the polarization curve of a single cell is plotted in Figure 2. In reality, the current density could be controlled within a certain range. After excluding the very low current densities ($i_{FC}<0.1$ A/cm²), (13) could be linearized [34,35] as:

$$U_{FC}=U_{OC}-{\tilde{R}}_{FC}{i}_{FC},$$

(17)

where $U_{OC}$ is the voltage at which the linearized curve crosses the y-axis, which should not be confused with $U_{rev}$.

Unlike ground vehicles, the sUAV changes its orientation during the flight, which would change the inner resistance of the fuel cell by about five times [36] from horizontal to vertical. To this end, (17) is modified to account for this effect as:

$$U_{FC}=U_{OC}-{\tilde{R}}_{FC}^{\prime }{i}_{FC},$$

(18)

where ${\tilde{R}}_{FC}^{\prime }={\tilde{R}}_{FC}(1+k_{0}sin(k_1\left|\gamma \right|\left)\right)$. Combining Equation (18) with Equations (10)–(12) yields:

$$n_{FC}{R}_{FC}^{\prime }{I}_{FC}^2-(n_{FC}{U}_{OC}-n_{FC}{\kappa }_{aux})I_{FC}+P_{FC,load}+P_0=0,$$

(19)

where $R_{FC}^{\prime }={\tilde{R}}_{FC}^{\prime }/A_{FC}$. Overall, $I_{FC}$ can be expressed as:

$$\begin{array}{cc}\hfill I_{FC}=& \frac{n_{FC}(U_{OC}-{\kappa }_{aux})-{\left[n_{FC}^2{(U_{OC}-{\kappa }_{aux})}^2-4n_{FC}{R}_{FC}^{\prime }(P_{FC,load}+P_0)\right]}^{\frac{1}{2}}}{2n_{FC}{R}_{FC}^{\prime }}.\hfill \end{array}$$

(20)

2.5. Battery Model

The battery model represents a pack of Model 21700 lithium polymer battery cells. The battery pack is assembled in such a way that the cells are connected in series. According to [37], the open-circuit voltage of the battery can be estimated as:

$$U_{B,OC}=SOC(U_{B,max}-U_{B,min})+U_{B,min}.$$

(21)

The battery power and the fuel cell load demand power sum up to provide the electrical power to the motor such that:

$$P_M=P_B+P_{FC,load}.$$

(22)

Further, the current drawn from the battery set is obtained by solving:

$$P_B=n_B(U_{B,OC}{I}_B-I_B^2{R}_{B,int}),$$

(23)

which should not exceed its maximum discharge current $I_{B,max}$.

The battery Coulombic efficiency in the battery model is assumed to be 100%. Thus, the $SOC$ is satisfied as:

$$SOC\left(t\right)=SO{C}_0-{\mathsf{\int }}_{t_0}^t\frac{I_B\left(t\right)}{C_B}dt.$$

(24)

where t, $t_0$, and $SO{C}_0$ are the current time, initial time, and initial $SOC$, respectively.

3. Hybrid System Model and Problem Formulation

3.1. Hybrid System

The fuel cell load demand power, which will be indicated as $P_{FC}$ in the following section, is the only variable under control. Due to the output characteristic of the PEMFC, the change of $P_{FC}$ is chosen to be 5% of its maximum power, which depends on $\gamma$ according to (18). The fuel cell load demand power dynamics are then:

$$P_{FC}\left(t_{n+1}\right)-P_{FC}\left(t_n\right)=u·\mathsf{\Delta }{P}_{FC},$$

(25)

where $u\in \{-1,0,1\}$, $\mathsf{\Delta }{P}_{FC}=5\%P_{FC,max}$, and $P_{FC}\left(t_n\right)$ is the fuel cell load demand power at $t=t_n$. Here, three different values of u correspond to decreasing, sustaining, or increasing $P_{FC}$. According to Equation (19), the maximum load cell power can be calculated as:

$$P_{FC,max}=\frac{4n_{FC}{R}_{FC}^{\prime }{P}_0-{(n_{FC}{U}_{OC}-n_{FC}{\kappa }_{aux})}^2}{-4n_{FC}{R}_{FC}^{\prime }}.$$

(26)

Using Equations (25) and (26), the final expression representing the fuel cell load demand power dynamics is given by:

$$\frac{d{P}_{FC}}{dt}=u·5\%\frac{4n_{FC}{R}_{FC}^{\prime }{P}_0-{(n_{FC}{U}_{OC}-n_{FC}{\kappa }_{aux})}^2}{-4n_{FC}{R}_{FC}^{\prime }}.$$

(27)

The $SOC$ dynamics are obtained by differentiating both sides of (24) with respect to time,

$$\frac{dSOC}{dt}=-\frac{I_B}{C_B}.$$

(28)

Combing (28) with (21) and (23) yields,

$$\frac{dSOC}{dt}=\frac{-n_B{U}_{B,OC}+\sqrt{{\left(n_B{U}_{B,OC}\right)}^2-4n_B{R}_{B,int}{P}_B}}{2n_B{R}_{B,int}{C}_B},$$

(29)

where $U_{B,OC}=SOC(U_{B,max}-U_{B,min})+U_{B,min}$.

The motor power and battery power are related by:

$$P_B=S_f{P}_M,$$

(30)

where $S_f$ is referred to as the split fraction, which could be calculated from (22) as:

$$S_f=\frac{P_M-P_{FC}}{P_M}.$$

Using Equations (29) and (30), the final expression representing the $SOC$ dynamics is given by:

$$\frac{dSOC}{dt}=\frac{-U_{B,OC}+\sqrt{{\left(U_{B,OC}\right)}^2-\frac{4R_{B,int}{S}_f{P}_M}{n_B}}}{2R_{B,int}{C}_B},$$

(31)

where the internal resistance $R_{B,int}$ and the battery capacity $C_B$ are assumed to be constant [38].

The mass of remaining fuel dynamics is obtained from Faraday’s law as:

$$\frac{d{m}_{FR}}{dt}=-\frac{n_{FC}{I}_{FC}}{n_{e}F}{M}_h,$$

(32)

where $I_{FC}$ is calculated from $P_{FC}$ as shown in (20).

Equations (25), (29), and (32) are the final form of the state equations used in this study, where the states of the system are $x=[SOC,M_{FR},P_{FC}]$, the control is $u\in \{-1,0,1\}$, the outputs of the system are $y=[S_f,P_B]$, and the operating variables are $w=[P_M,\gamma ]$. These operating variables are treated as measured disturbances in the model.

Based on the above modeling assumptions and parameters in Table A1, the maximum fuel cell output power is 795 W at $\gamma =0$ deg, $496.14$ W at $\gamma =\pm 10$ deg, and $335.71$ W at $\gamma =\pm 20$ deg. The theoretical maximum power for the battery series (of eight batteries) is 2940 W, due to the limitation of the discharge current (35 A); the maximum power of the battery is 1176 W at any climb angle.

3.2. Problem Formulation

The forward Euler method is used in this paper to approximate the time derivatives. During each time segment $\mathsf{\Delta }t$, the motor power of the sUAV is $w_1$ and the climb angle is $w_2$.

The following updated equations approximately model the sUAV hybrid propulsion system:

$$\begin{array}{cc}\hfill SOC\left(t_{n+1}\right)& =SOC\left(t_n\right)+\frac{dSOC}{dt}\left(t_n\right)\mathsf{\Delta }t,\hfill \\ \hfill M_{FR}\left(t_{n+1}\right)& =M_{FR}\left(t_n\right)+\frac{d{M}_{FR}}{dt}\left(t_n\right)\mathsf{\Delta }t,\hfill \end{array}$$

where $SOC\left(t_n\right)$ and $M_{FR}\left(t_n\right)$ are the state of charge and the mass of hydrogen remaining at $t=t_n$.

The system is controlled by the change of the fuel cell load demand power $\mathsf{\Delta }{P}_{FC}$ at each discrete time instant. Thus, the fuel cell power is modeled as:

$$P_{FC}\left(t_{n+1}\right)=P_{FC}\left(t_n\right)+u\mathsf{\Delta }{P}_{FC}\left(t_n\right).$$

The motor power and climb angle are typically unknown a priori. In this paper, a Markov chain model is used to describe the evolution of $w_1$ and $w_2$ with the transition probabilities identified from the historical data. Once particular $w_1$ and $w_2$ values are encountered, a prediction of their probability distribution over the next time segment will be made using the Markov chain model.

The objective of the stochastic endurance maximization problem is to determine a control law that maximizes the time the sUAV can travel before the system states exit a prescribed set,

$$\begin{array}{cc}\hfill G=& \{(SOC,M_{FR},P_{FC}):SO{C}_{min}\le SOC\le SO{C}_{max},\hfill \\ \hfill & \phantom{=}{M}_{FR,min}\le M_{FR}\le M_{FR,max},\phantom{=}0\le P_{FC}\le P_{FC,max}\}.\hfill \end{array}$$

(33)

The constraints on the $SOC$ and $M_{FR}$ in Equation (33) reflect the minimum and maximum values of the battery state of charge and the mass of fuel, respectively. The constraints on $P_{FC}$ are reflective of the fact that the fuel cell load demand power cannot (i) exceed the maximum power of the fuel cell and (ii) be negative.

The optimal control policy developed in this paper through the application of DCOC specifies the change in fuel cell load power over one step, $\mathsf{\Delta }{P}_{FC}\left(t\right)=P_{FC}(t+1)-P_{FC}\left(t\right)$, as a function of $SOC\left(t\right)$, the mass of hydrogen fuel left, $M_{FR}\left(t\right)$, and the current fuel cell load power, $P_{FC}\left(t\right)$. The battery power complements fuel cell power in matching propeller requested power.

3.3. Markov Chain Modeling

A Markov chain model [39] is used to represent the evolution of $w$ (in our case, $w=[w_1,w_2]$). The transition probabilities of the Markov chain are defined as:

$$p_{ij}=prob\{w\left(t_{n+1}\right)\in W_j|w\left(t_n\right)\in W_i\},$$

(34)

where $W_i$ and $W_j$ ($i,j=1,\cdots ,N$) are cells partitioning the feasible range of the operating conditions. The state dependence of the transition probabilities adds flexibility in reflecting the typical motor power and climb angle profiles of an sUAV.

The $p_{ij}$’s can be obtained from the statistical analysis of the historical flight data,

$$p_{ij}=\frac{M_{ij}}{M_i},$$

(35)

where $M_{ij}$ is the total number of transitions from the cell $W_i$ to the cell $W_j$ (i.e., $\mathit{w}\left(t_n\right)\in W_i,\mathit{w}\left(t_{n+1}\right)\in W_j$), while $M_i$ is the total number of transitions from $W_i$ to any other cell, including $W_i$ [21].

4. Control Law Construction

Here, we adopt the SDCOC framework from [16], which is applied to a discrete-time model with the following form,

$$x\left(t_{n+1}\right)=f(x\left(t_n\right),u\left(t_n\right),w\left(t_n\right)),$$

(36)

where $x\left(t_n\right)$ is the state vector, $u\left(t_n\right)$ is the control vector, and $w\left(t_n\right)$ is the vector of operating variables, which is not known until the time instant $t_n$. The system has both control constraints and state constraints imposed as $u\left(t_n\right)\in U$ and $\{x\left(t_n\right),w\left(t_n\right)\}\in G$, respectively, where U and G are specified sets. A Markov chain with a finite number of states is used to represent transitions in $w\left(t_n\right)\in W=\{w^p:p\in P\}$. Here, P is the size of the grid for $w$. The transition probability from $w\left(t_n\right)=w^i\in W$ to $w\left(t_{n+1}\right)=w^j\in W$ is denoted by $p_{ij}$, expressed in (34). In a discounted variant of SDCOC, the objective is to determine a control function $u(x,w)$ such that, with $u\left(t_n\right)=u(x\left(t_n\right),w\left(t_n\right))$, a cost functional of the form,

$$J^{x_0,w_0,u}={\mathbb{E}}_{x_0,w_0}\left[\sum _{t=0}^{{\tau }^{x_0,w_0,u}\left(G\right)-1}{\delta }^t·1\right],$$

(37)

is maximized. Here, ${\tau }^{x_0,w_0,u}\left(G\right)\in {\mathbb{Z}}^{+}$ represents the first time instant when the trajectories of $x\left(t_n\right)$ and $w\left(t_n\right)$, which are denoted by $\{x^u,w^u\}$ and result from applying the control $u\left(t_n\right)=u(x\left(t_n\right),w\left(t_n\right))$ with values in the set U, exit the prescribed compact set G. $\delta$ is a discount factor [40]. For $\delta =1$, (37) maximizes the exit time, i.e., the time till the prescribed constraints become violated. The use of the discount factor $0<\delta <1$ facilitates faster convergence of the value iterations. Note that $\{x^u,w^u\}$ is a random process, ${\tau }^{x_0,w_0,u}\left(G\right)$ is a random variable, and ${\mathbb{E}}_{x_0,w_0}[·]$ denotes the conditional expectation given the initial values of $x$ and $w$.

To solve (37), the value iterations approach is used, which produces a sequence of value function approximations, $V_n$, at specified grid-points $x\in \{x^k:k\in K\}$,

$$\begin{array}{cc}\hfill V_0(x^k,w^i)& \equiv 0,\hfill \\ \hfill V_n(x^k,w^i)& =\underset{u^m,m\in M}{max}\left\{\sum _{j\in J}{F}_{n-1}(f(x^k,u^m,w^i),w^j)·p_{ij}·{\delta }^t+1\right\},\hfill \end{array}$$

where $u\in \{u^m:m\in M\}$ is a specified grid for $u$. Here, K and M are the size of the grid for $x$ and $u$, respectively. In each iteration, once the values of $V_{n-1}$ at the grid-points have been determined, linear or cubic interpolation is employed to approximate $V_{n-1}(f(x^k,u^m,w^i),w^j)$ as $F_{n-1}(x,w^i)=\mathrm{Interpolate}\left[V_{n-1}\right](x,w^i)$, if $(x,w^i)$$\in G$, and $F_{n-1}(x,w^i)=0$, if $(x,w^i)\notin G$. A termination criterion of the form $|V_n(x,w^i)-V_{n-1}(x,w^i)|\le ϵ$ for all $x\in \{x^k:k\in K\}$ and $i\in P$, where $ϵ>0$ is sufficiently small, is used.

Once an approximation of the value function, $V_{*}$, is available, an optimal control law is determined as:

$$\begin{array}{cc}\hfill u_{*}(x,w^i)\in \{& u:V_{*}(x,w^i)-\sum _{j\in J}{V}_{*}(f(x,u,w^i),w^j)·p_{ij}·\delta -1\le ϵ\}.\hfill \end{array}$$

5. Control Law Computations and Results

5.1. sUAV Configuration and Model Parameters

The model was parameterized for a $1.5$$\mathrm{kg}$ sUAV [23] that can be used for aerial photography and environmental monitoring applications. The minimum and maximum $SOC$ values were set to $SO{C}_{min}=0.2$ and $SO{C}_{max}=0.8$. The minimum and maximum values of $M_{FR}$ were set as $M_{F{R}_{min}}=2$$\mathrm{g}$ and $M_{F{R}_{max}}=9$$\mathrm{g}$. For the value iterations, the $SOC$ grid was chosen with a step size of $0.05$, and the $M_{FR}$ grid was chosen with a step size of $0.5$$\mathrm{g}$. The grid for the control variable u was set as $\{-1,0,1\}$.

The transition probabilities for the operating variables (motor power and climb angle) were obtained from the time histories of the sUAV motor power and climb angle using (35) and assuming a time step $\mathsf{\Delta }t=1$ s. These time histories were based on a scenario in which an sUAV follows a moving ground vehicle that sUAV operators are interested in monitoring. In this scenario, the ground vehicle, and consequently the sUAV, is assumed to be traveling with the velocity profile defined by concatenating the EPA Highway Cycle [41] nine times. For the sUAV, the speed profile is modified to remain above the stall speed while avoiding extreme acceleration values.

The climb angle time history, shown in Figure 3, was obtained from the Google Earth elevation profile for a path from Monroe, West Virginia, to Princeton, West Virginia, with the help of GPS visualizing software [42]. See [43] for the assessment of the accuracy of such extracted profiles.

Figure 4 provides the time histories of the sUAV motor power calculated based on the equations in Section 2.3. The trajectories in Figure 3 and Figure 4 were used to compute the transition probabilities.

5.2. Control Law Computation

Cython was used for control law computations as it is more efficient than MATLAB in handling nested for loops and two-dimensional interpolation. In our numerical experiments with dynamic programming, Cython was about 10 times faster than MATLAB.

To further speed up the value iterations, a discount factor was introduced. When testing the effect of the discount factor on the optimal policy, a zero climb angle ($\gamma =0$) was assumed, which means that the only operating variable was the motor power.

Table 2. Average exit time for different discount factor.

$\mathit{\delta }$	Number of Iteration	Computing Time (min)	Exit Time with 20% Initial $\mathit{SOC}$ (s)	Exit Time with 80% Initial $\mathit{SOC}$ (s)
0.99	2258	830.02	6358.44	9742.99
0.97	753	100.69	6276.18	9716.86
0.95	448	58.27	6221.37	9640.24
0.93	317	39.22	6186.50	9610.22
0.91	244	30.55	6159.65	9602.55

shows the average exit time based on 100 random simulations for discount factors from $0.91$ to $0.99$. The stopping criterion was chosen with $ϵ={10}^{-10}$ for all $\delta$. Computations were performed on a Hasee K780G-i7 laptop with a CORE i7-4710MQ (2.5–3.5 GHz) processor and 24 GB of RAM.Note that the number of iterations and the computing time decrease as the discount factor decreases, but so does the exit time. The discount factor $\delta =0.95$ was ultimately chosen as a compromise between value iteration convergence speed and solution accuracy. Figure 5 shows that the value iterations with a discount factor of $\delta =0.95$ converge much faster than those with $\delta =1$.

5.3. Endurance Maximization Results

We used $ϵ={10}^{-10}$ in the stopping criterion for the value iterations. Figure 6 illustrates the resulting control policy. Note that when $SOC$ is low, the control policy calls for an increase in $P_{FC}$ to charge the battery. This is reasonable given that the fuel cell cannot alone respond rapidly to fast changes in motor power request, and hence, the battery has to be charged to do so.

The simulation results are given for three cases in Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18. The first case (Scenario I) corresponds to a higher initial $SOC$, and the second case (Scenario II) considers a lower initial $SOC$. The third scenario is for a mid-range initial $SOC$ and is used to confirm the $SOC$ behavior observed in the first two scenarios. In all cases, the initial fuel mass and initial fuel cell power are the same: $M_{FR,0}$ = 6 $\mathrm{g}$ and $P_{FC,0}$ = 0 $\mathrm{W}$. The dashed lines in Figure 8, Figure 12, and Figure 16 indicate the constraints mentioned in Section 5.1. The spikes of power in Figure 7, Figure 11, and Figure 15 correspond to the time instants when the sUAV starts to accelerate while the positive and negative spikes of climb angle represent the time when the sUAV starts to climb or descend.

Figure 7, Figure 8, Figure 9 and Figure 10 illustrate the closed-loop response for the first simulation scenario. The initial $SOC$ is $0.8$, and it decreases rapidly until it reaches a value of about $0.5$. Then, it stays near this value between 2000 and 5000 $\mathrm{s}$. Finally, when the mass of hydrogen reaches a relatively low value, the $SOC$ starts to decrease and continues to decrease until the constraints are violated. The fuel cell load demand power keeps a relatively low value during the whole flight, and the mean value of the split fraction is negative during the 2000 to 5000 $\mathrm{s}$ time interval, which is the period when the $SOC$ is kept at about $0.5$.

Figure 11, Figure 12 and Figure 13 illustrate the closed-loop response for the second simulation scenario. The initial $SOC$ is $0.2$. The battery is charged until it reaches a value of about $0.5$ to enable the battery to sustain rapid propeller power fluctuations. Then, the $SOC$ stays near that value of $0.5$ between 500 and 1500 $\mathrm{s}$. Finally, when the mass of hydrogen reaches a relatively low value, the $SOC$ starts to decrease and continues to decrease until the constraints are violated. The fuel cell load demand power increases rapidly at first to charge the battery, then it keeps a relatively low value during the rest of the flight. The mean value of the split fraction is negative from the beginning to about 1500 $\mathrm{s}$, which is the period when the battery is charged from $SOC=0.2$ to about $0.5$.

According to the results from Scenarios I and II, a turnpike behavior of the battery $SOC$ is observed, with the $SOC$ converging to about $0.5$ and staying at that value for a while before decaying. To confirm this turnpike behavior, we additionally considered the responses with the developed policy to the initial SOC of $0.6$. These are shown in Figure 15, Figure 16, Figure 17 and Figure 18. The exit times for Scenarios I, II, and III were, respectively, 9890 $\mathrm{s}$, 6019 $\mathrm{s}$, and 8478 $\mathrm{s}$.

6. Conclusions

This paper considers an endurance maximization problem for a small unmanned aerial vehicle (sUAV) with a hybrid propulsion system consisting of a polymer electrolyte fuel cell and a battery, both driving an electric motor connected to a propeller. A stochastic drift counteraction optimal control (SDCOC) approach is employed to develop control policies for optimally coordinating the fuel cell and the battery while enforcing the constraints on the fuel cell power output rate of change. Cython is used to implement value iterations and demonstrated an order of magnitude speedup versus MATLAB without increasing the code complexity, due to its efficiency in handling nested for-loops. Additionally, the use of a discount factor is shown to significantly speed up the value iterations at the price of decreased performance. The results illustrate the effectiveness of the SDCOC strategy in regulating the charging behavior of the battery by the fuel cell to provide the capability to respond to rapidly varying motor power demand.

The proposed approach based on SDCOC is particularly suitable for handling stochastic disturbances and can be applied to sUAVs exposed to headwinds with the headwind modeled as a stochastic disturbance. Accounting for such wind disturbances, extensions to include thermal dynamics, systematic and comprehensive comparison with other energy management approaches and propulsion system choices, the systematic study of the robustness to model uncertainties, as well as actual flight experiments represent directions for continuing research. In particular, our study of the discount factor impact on the computation time and exit time suggests that the flight time is sensitive to the choice of the energy management strategy; our approach based on SDCOC is optimal in the sense of maximizing expected flight time within a stochastically modeled environment.