Marie Truffot
email: marie.truffot@centralesupelec.fr

Axel Ivaldi

Antoine Renaud

Laurent Zimmer

Franck Richecoeur

Bistability of a modulated downwards oriented laminar flame

Introduction

In recent reacting systems, the fuel's multi-injection leads to a more uniform distribution while limiting long mixing chambers. For sustainable flames, a small central pilot is used to enable easier anchoring while preventing blow-off issues. However, this double-injection system has some drawbacks as potentially several flame shapes may be encountered [START_REF] Renaud | Bistable behaviour and thermo-acoustic instability triggering in a gas turbine model combustor[END_REF]. It is then of profound interest to understand the mechanisms for shape transition. Most of the flame transitions happen for turbulent flames that make complex the fundamental investigation of the dynamics. In the present study, we propose a laminar configuration that presents a bistable behavior similar to dynamics observed in high-power turbulent combustors. The burner's dynamic is investigated here with a two-step method. First, optical flow processing is applied to the 2D Mie scattering images to render accurate velocity fields minimizing the discrepancy compared to correlation algorithms [START_REF] Quénot | Particle image velocimetry with optical flow[END_REF]. Second, we extract the POD (Proper Orthogonal Decomposition) modes from these fields and plot the most energetics spatial modes and phase portraits in a configuration where multiple resonant frequencies stem from the harmonic excitation.

Experimental set-up

The overall dynamic of the Incopladis burner developed at the EM2C laboratory is studied. The burner is an inverted bunsen burner type with a heat release of around 800W. The experiment consists of a downward directed injection with an internal 20mm radius. A gaseous methane-air premix is injected with a global equivalence ratio of 0.91 and a typical Reynolds number around 1200, which corresponds to a mass flow of 0.29g/s for the air and 0.015g/s for the methane.

Figure 1. Inverted Bunsen burner set-up

The burner is protected from radiation with a shield and, water flow is used around the injection system to have a stable thermal condition. A servo proportional valve (Moog) is used to induce mass flow fluctuation by modulating the outlet with a sinusoidal signal with an offset of 5V and an amplitude of 5V. The valve is fully closed for a signal of 0V and fully opened for a signal of 10V and has a frequency band ranging from 0Hz to 150Hz. The modulation of the air mass flow is acoustic beyond 50Hz.

This burner is equipped with several diagnostics to allow the study of its dynamics. Initially, heat release measurements are tracked with a photomultiplier tube at a sampling frequency of 10 000Hz. Subsequently, a continuous laser is used to measure the instantaneous velocity fields of fresh gases with oil particles injected into the flow. The oil droplets have a typical size below 2µm, leading to a small Stokes number. The flow velocities are measured using 2D Mie scattering images for a forcing modulation rate of 70H and 85Hz. A highspeed camera takes the images at a rate of 10000Hz. The camera has a resolution of 896x848 pixel, and images are acquired at 90° from the laser sheet.

This study focuses on the flame responses characterization to a flow modulated with a sinusoidal excitation. An optical flow algorithm (derived from [START_REF] Sun | A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them[END_REF]) is used to extract the velocity from 2D Mie scattering images. Furthermore, velocity fields are post-processed using Proper Orthogonal Decomposition (POD).

Results

The lean premixed laminar flame operates first in a stationary regime then harmonic flow modulations induce a dynamic response. We want to show here that the flame' response at 70 and 85 Hz is drastically different and that in between these two frequencies, the flame features a bistable dynamic.

Figure 2. Photo-multiplier signal recorded during a fiveminute long sequence (blue). The moving average (red) shows the alternation of the flame emission intensity while the modulation frequency remains constant.

Indeed, by modulating the mass flow rate with frequencies between 70 and 85 Hz, the flame spontaneously alternates between two steady states A and B. This alternation occurs without external forcing or modification of any parameter of the experiment. The temporal evolution of the flame luminescence, recorded by the photo-multiplier (Figure 2), shows intermittent switches between two average values 𝐼 ! and 𝐼 " . Each state is reached several times during a sequence and can last several tens of seconds. The sequence presented here lasts 3 minutes; the moving average of this signal is plotted in red.

The signal's PDF confirms that the flame emission is mainly distributed on one set or the other. The flame alternates between two combustion configurations, the choice between these two states is not controlled by an experimental parameter, the transition from one state to the other is spontaneous. The residence times on each plateau are long enough to observe the average flame shape and measure the average velocity field associated with these two situations. We do not bring here elements explaining the origin of these spontaneous transitions, but we wish to take the time to describe the two possible states for these flames. These two states differ in terms of average velocity fields in the vicinity of the flame and the dynamics of flame response.

Figure 3. Vertical velocity fields for a forcing frequency of 70Hz (left section: state A) and 85Hz (right section: state B)

State A corresponds to the flame observed when modulated at 70 Hz, it features a compact flame anchored to the lips of the injector, and the average velocity field is similar to that observed around a circular laminar burner, its average shape is conical as shown in the left part in Figure 3Erreur ! Source du renvoi introuvable.. State B corresponds to the flame observed when modulated at 85 Hz, the flame is lifted with a fresh gas area extended downstream (right part in Figure 3). The shape of the flame is open, and the baseline is flat. This configuration is original insofar as the flow remains laminar and not swirling.

When modulated at frequency equals to 70 Hz, the flame's dynamic response looks consistent with what the literature reports for laminar flames [START_REF] Bourehla | Appearance and Stability of a Laminar Conical Premixed Flame Subjected to an Acoustic Perturbation[END_REF]. The laminar flame is anchored to the burner lips, and the flow modulation gives rise to vortices in the injection plane that are convected by the flow, then deform the flame front periodically. When modulated at 85 Hz, the flame looks different. To observe the flame's response in both states, the velocity field is decomposed into orthogonal modes (POD). In both situations, five modes contain most of the kinetic energy. The spectra and spatial distributions of these five modes for the two flame states are presented in Figure 5, Figure 7 and Figure 8.

Figure 4. Phase Portrait of the POD mode 1 and 2 for a forcing modulation frequency of 70Hz (state A)

At 70Hz (state A), the spectra of modes 1 and 2 contain a single frequency, corresponding to the flame's response at the modulation frequency. The phase portrait of these two modes (Figure 4) shows that they are coherent in phase quadrature. They correspond to the propagation of the vortices created at the injector lips.

At 85 Hz (state B), modes 2 and 3 correspond to the same dynamics. Their spectra also show only the excitation frequency. Their phase portrait (Figure 6) shows the coherence and phase shift corresponding to the propagation along the vortices' flame front, as witnessed previously.

Figure 6. Phase Portrait of the POD mode 2 and 3 for a forcing modulation frequency of 85Hz (state B)

Concerning the other orthogonal modes, the 70 Hz modulation shows a non-identified stationary mode 3 while modes 4 and 5 spectra feature a single peak at twice the excitation frequency. The phase shift and coherence of these latter two modes indicate that this frequency doubling corresponds to half-size vortices propagating along the flame front (Figure 5, top line), as shown by these modes' spatial velocity distribution.

Figure 5. POD modes of the vertical velocity superposed with the streamline associated with the mode. First row: state A (modulation frequency of 70Hz). Second row: state B (modulation frequency of 85Hz) Figure 7. Power density spectrum of the first five most energetic POD modes for a forcing modulation frequency of 70Hz (state A)

Modes 1 and 4 observed during modulation at 85Hz are original in this configuration. These two modes have a denser and broader spectrum at low frequencies (Figure 8) but are not coherent with each other. They do not have a defined periodicity related to the modulation frequency. The spatial velocity distribution corresponding to these modes shows that they correspond to transverse movements (Figure 5, bottom line).

To understand the mode 1, we map the mode that reveals an asymmetry in the flow generating large diameter recirculations. This mode moves the flow transversely at low frequencies. Mode 4 shows an axial shear layer inducing a sinuous mode where the fluid is alternately deviated transversely. These transverse phenomena occur at a range of frequencies ten times smaller than the excitation frequency; they appear when the flame is detached (state B). In the first state, the flame is attached, its emission intensity fluctuates like the modulation. In the second state, the flame is detached and elongated, and two dynamics are superimposed: the first at the frequency of the excitation deforms the flame front, the second at a lower frequency generates sinuous modes keeping the flame lifted and causing the flame to oscillate.

For a given frequency range, the two states can be obtained alternately without any external control and without observing any apparent periodicity. Modulation of the mass flow rate appears to induce secondary instability at a lower frequency sufficiently intense to alter the average flame state and cause the flow to oscillate in the transverse direction.

Conclusions

Several recent studies have reported the existence of bi-stable swirled turbulent burners, i.e., burners with two modes of flame stabilization that can be achieved either by variation of an experimental parameter (staging, flow rate, equivalence ratio, and so) often presenting a hysteresis phenomenon or by spontaneous alternation depending on the dynamics of the system. This study shows that such a bi-stable system can appear in a pulsed laminar configuration where the Bunsen-like burner is oriented downwards. A lean methane/air flame is submitted to a flow modulation in a frequency range from 70 to 85 Hz. For a given frequency and during the same sequence, the flame alternately exhibits two different responses characterized by distinct average positions and dynamics. The study presented here relies on the velocity fields obtained by optical flow, the orthogonal mode processing of these fields, and the flame's spontaneous emission to describe the two states reached by the flame over time. The frequency and spatial dynamics are very different between the two states, confirming two stable states' presence reached alternately without evident periodicity. Obtaining a bistable regime in an elementary laminar configuration is the originality that makes it possible to perform accurate experimental measurements and modeling of these situations frequently encountered in more complex burners to understand the transitions.

Figure 8 .

 8 Figure 8. Power density spectrum of the first five most energetic POD modes for a forcing modulation frequency of 85Hz (state B)

Acknowledgements

The support of Safran Aircraft Engine and ANRT for the Ph.D. grants of M. Truffot and A. Ivaldi is gratefully acknowledged. The authors would also like to thank JM Dupays for helping in the management of the supercomputer ressources (Ruche Mesocentre). This work was performed using HPC resources from the "Mésocentre" computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/).


Member, IEEE Gennaro Notomista
email: gennaro.notomista@irisa.fr

Senior Member, IEEE Claudio Pacchierotti
email: claudio.pacchierotti@irisa.fr

Senior Member, IEEE Paolo Robuffo Giordano

Online Robot Trajectory Optimization for Persistent Environmental Monitoring

Keywords: Robotics, Energy systems, Constrained control

This paper presents a control-theoretic approach for trajectory optimization of mobile robots suitable for environmental monitoring. The method is based on the optimization of the Constructability Gramian in order to maximize the information collected while traversing a state trajectory. The maximization of the collected information is combined with energy constraints to define an optimization-based controller that achieves persistent environmental monitoring. The results of its application to the estimation of the concentration of a diffusing gas using a mobile robot are shown in simulation.

I. INTRODUCTION

T HE use of robotic systems for environmental monitoring is becoming more and more common thanks to their suitability to collect spatial and temporally disperse data [START_REF] Dunbabin | Robots for environmental monitoring: Significant advancements and applications[END_REF]. However, the longer the time horizon over which robots need to remain operational to collect environmental data, the more challenging their deployment becomes. In particular, energy management and control is paramount when autonomous, untethered robotic platforms are considered [START_REF] Sadler | Fundamentals of energy-constrained sensor network systems[END_REF].

In this paper, we propose an online trajectory optimization strategy for mobile robots suitable for persistent environmental monitoring tasks. The approach consists of two main components: (i) Trajectory optimization to maximize the gathered information, framed as an active sensing problem [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF], and (ii) Energy control using physics-based robot battery model implemented using control barrier functions [START_REF] Ames | Control barrier functions: Theory and applications[END_REF]. These two components are combined in a single optimization-based controller which computes both the optimized trajectory and the control input required by the robot to track it.

Robots deployed for persistent environmental monitoring have been considered in [START_REF] Smith | Persistent robotic tasks: Monitoring and sweeping in changing environments[END_REF], where a spatial information field has to be swept by the robots whose velocity is controlled over a predefined trajectory in space. The employment of autonomous agents for environmental monitoring has been recognized to be particularly useful for source localization applications [START_REF] Dhariwal | Bacterium-inspired robots for environmental monitoring[END_REF], [START_REF] Trincavelli | Towards environmental monitoring with mobile robots[END_REF]. In [START_REF] Khodayi-Mehr | Model-based active source identification in complex environments[END_REF], mobile sensors are controlled to move in order to maximize the information collected required to identify the source of an advective-diffusive environmental field model, implementing, as a matter of fact, an active sensing control strategy.

Active sensing has been studied under different perspectives (see, e.g., [START_REF] Mihaylova | A comparison of decision making criteria and optimization methods for active robotic sensing[END_REF]- [START_REF] Ouyang | Multi-robot active sensing of non-stationary gaussian process-based environmental phenomena[END_REF]). In [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF], the authors show how a norm of the so-called Constructability Gramian (CG) is an appropriate measure of the amount of information gathered along a given trajectory of the state of a system. Leveraging this insight, in this paper we develop a persistent environmental monitoring strategy where the estimation of environmental fields is to be achieved under limited energy resources, which characterize the mobile robotic platforms employed for such tasks. Moreover, compared to empirical methods devised to maximize robot deployment time [START_REF] He | Autonomous chemical-sensing aerial robot for urban/suburban environmental monitoring[END_REF] or minimize energy expenditure [START_REF] Martin | Long-term exploration & tours for energy constrained robots with online proprioceptive traversability estimation[END_REF], we do so in a provable fashion by employing control barrier functions to enforce energy constraints.

The main contributions of this paper with respect to existing approaches are: (i) Extending the control-theoretic method proposed in [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF] to online robot trajectory planning for environmental monitoring; (ii) Considering the limited energy availability on the robots using an accurate battery model; (iii) Combining (i) and (ii) in an optimization-based controller to achieve a computationally efficient persistent environmental monitoring algorithm.

II. ACTIVE SENSING FOR ENVIRONMENTAL MONITORING

This section introduces the active sensing problem which will be extended to cover environmental monitoring scenarios. A control-theoretic approach based on the CG will be employed and combined with energy control to achieve persistent environmental monitoring. The solution proposed in this paper builds on the ideas originally proposed in [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF] which are briefly recalled in the next subsection for the reader's convenience.

A. Active Sensing

In active sensing tasks, robot trajectories are planned with the objective of maximizing the information gathered by the robots during motion. In order to quantify the gathered information, we start by defining models for robot motion and sensor measurement, as follows:

ẋ = f (x, u) + x y = m(x) + y . (1)
Here, f : R nx × R nu → R nx is a Lipschitz continuous vector field modeling the robot dynamics. The measurement model is given by m : R nx → R ny . In the context of persistent environmental monitoring considered in this paper, measurements y will be leveraged to produce an estimate of the state x using an extended Kalman filter. The system and measurement noise are modeled by zero-mean Gaussian distributions with covariance matrices Q and R, respectively, i.e. x ∼ N (0, Q) and y ∼ N (0, R).

In [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF], the authors propose a solution to the active sensing problem consisting in optimizing the state trajectory in order to maximize a norm of the CG for system [START_REF] Dunbabin | Robots for environmental monitoring: Significant advancements and applications[END_REF], which is defined over a time horizon [t s , t f ] as:

G c (t s , t f ) = t f ts Φ(τ, t f ) T C(τ ) T W (τ )C(τ )Φ(τ, t f )dτ, (2)
where Φ(τ, t f ) ∈ R nx×nx is the transition matrix of the dynamical system ẋ = f (x, u), C(τ ) = ∇ x m(x(τ )) T , and W is a weight matrix set to R -1 , the inverse of the covariance matrix of the measurement noise y [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF]. This way, a norm of the CG is a suitable metric for quantifying the information gathered while moving along the state trajectory x :

[t s , t f ] → R nx .
The infinite-dimensional trajectory optimization can be reduced to a finite-dimensional optimization problem by parameterizing the state trajectory x(t), t ∈ [t s , t f ] using a finite number of parameters x c ∈ R nx c . This can be done, for instance, by considering polynomial (such as cubic splines and B-splines) state trajectories. Then, the following gradientascent strategy can be employed to optimize the trajectory so to maximize the information gathered:

ẋc = K ∇ xc G c (t s , t f ), (3)
K being a positive gain. We refer to [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF] for the details about the evaluation of the gradient of the CG G c (t s , t f ) with respect to the trajectory parameters x c .

B. Environmental Monitoring as Active Sensing

In environmental monitoring tasks, robots are deployed to estimate an environmental field, such as temperature and gas concentration maps. In this paper, we assume the field of interest can be parameterized by a vector θ ∈ R n θ . Then, the robot and measurement model (1) can be extended as follows:

     ẋ = f (x, u) + x θ = g(θ, t) + θ y = m(x, θ) + y . ( 4
)
The dynamics of the environmental field-which need not be static-are modeled by the vector field g : R n θ × R → R n θ assumed to be Lipschitz continuous in the first argument, θ, and continuous in its second argument, t. This choice allows us to consider dynamically evolving environmental fields (such as gas concentration, solar light intensity, and air properties), and it is suitable to represent discretized solutions of partial differential equations modeling an environmental field. The formulation in (4) allows us to treat environmental monitoring tasks as instances of active sensing. Following what has been proposed in [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF], a trajectory for x can be planned in order to maximize the information gathered through the measurements y required to estimate θ.

Remark 1. In the model (4), the state is comprised of both x and θ. However, the control input u influences directly only the state of the robot, x, while θ evolves independently. Therefore, a state trajectory can be planned considering just the components x of the state of the system (4). It is worth, however, pointing out that θ is nevertheless taken into account in the state trajectory optimization through the gradient of G c (t s , t f ) which does depend on θ.

The formulation in (4) allows us to execute environmental monitoring tasks by optimizing the state trajectory. In the next section, we will present the battery charge control as a way of rendering environmental monitoring tasks persistent, so to allow robots to execute them over long time horizons.

III. BATTERY CHARGE CONTROL

Environmental monitoring tasks may take place over long time horizons depending on the nature of the phenomenon of interest that needs to be monitored. In this paper, we present an approach to perform energy-constrained environmental monitoring which consists of an optimization-based strategy for trajectory optimization constrained by the limited availability of energy that robots can carry around.

In [START_REF] Notomista | Persistification of robotic tasks using control barrier functions[END_REF] and [START_REF] Notomista | Persistification of robotic tasks[END_REF], the persistification of generic robotic tasks is considered, through which the energy of the robots executing the tasks is controlled to remain about a lower bound. However, the considered battery model is conservative as it always takes into account the worst-case scenario of maximum discharge rate that the battery of the robot can undergo. In this section, we refine the persistification strategy by introducing a more faithful energy model that describes the dynamics of the state of charge (SOC) of a battery. Then, we introduce integral control barrier functions (I-CBFs) as a way of enforcing energy constraints. The SOC dynamic model will be used in conjunction with I-CBFs in order to ensure that the battery of the robots executing a task never depletes-i.e. the SOC never falls below a minimum threshold-while executing a task of interest, in this case the trajectory optimization task to realize environmental monitoring through active sensing.

A. State of Charge Dynamics

The state of charge is used to describe the remaining capacity of a battery, and it is therefore a very important parameter to take into account while designing energy-aware control strategies [START_REF] Piller | Methods for state-of-charge determination and their applications[END_REF]. Its accurate knowledge allows us to design algorithms to protect batteries from overdischarge and overcharge, increasing their life. However, batteries store chemical energy, which cannot be directly measured and has to be, therefore, estimated. The SOC estimation still remains a fundamental challenge, as it hinges on battery models which, in turn, depend on many parameters [START_REF] Krewer | Dynamic models of Li-Ion batteries for diagnosis and operation: A review and perspective[END_REF]. A survey of methods which have been proposed to estimate the SOC of batteries can be found in [START_REF] Chang | The state of charge estimating methods for battery: A review[END_REF], [START_REF] Murnane | A closer look at state of charge (SOC) and state of health (SOH) estimation techniques for batteries[END_REF]. In the following, we present the basic SOC dynamic model (found, e.g., in [START_REF] Zhao | Least-squares based Coulomb counting method and its application for state-of-charge (SOC) estimation in electric vehicles[END_REF]) which can be employed to render the environmental monitoring task persistent, by preventing robot batteries from depleting.

Let e(t) denote the value of the SOC of the robot battery at time t, C its rated capacity, η its Coulombic efficiency, and i(t) the current flowing from the battery-i(t) > 0 if the current is flowing out of the battery at time t, i(t) < 0 if it is flowing in. Then, the value of e(t) can be evaluated as follows:

e(t) = e(t 0 ) -C -1 t t0 ηi(t)dt, (5)
where, in general, the efficiency η can be a function of the current i(t), the temperature, the SOC itself, as well as the state of health (SOH) of the battery [START_REF] Zhao | Least-squares based Coulomb counting method and its application for state-of-charge (SOC) estimation in electric vehicles[END_REF].

In a large variety of robotic applications, it is reasonable to assume that the current supplied by the battery is proportional to the magnitude of the control input u supplied to the robot. This holds, for instance, in cases where electric motors and actuators-for which exerted torque is proportional to absorbed current-are employed for motion and locomotion purposes [START_REF] Toliyat | Handbook of electric motors[END_REF]. Therefore, we can model the relation between current and robot input as i(t) = α i u 2 , where α i is a class K function.

Starting from (5), and using this assumption, yields the following SOC dynamic model for the battery of robot i:

ė = -ηC -1 α i u 2 . (6)

B. Integral Control Barrier Functions

The model of the battery SOC dynamics in (6) depends on the control input u, and in general it does so in a noncontrol-affine fashion. As such, a straightforward application of CBF-based approaches would result in differential constraints which are not affine in u [START_REF] Ames | Control barrier functions: Theory and applications[END_REF]. As will be shown in the next section, affine input constraints are desirable in order to synthesize robot controllers under real-time constraints at high frequencies commonly used in the control loops of modern robotic platforms. To circumvent this issue, control dependent control barrier functions [START_REF] Huang | Guaranteed vehicle safety control using control-dependent barrier functions[END_REF] or the more general formulation of integral control barrier functions [START_REF] Ames | Integral control barrier functions for dynamically defined control laws[END_REF], defined below, can be employed.

Definition 2 ( [24]

). For the system ẋ = f (x, u), with corresponding safe set S ⊂ R nx × R nu defined as the 0superlevel set of a function h : R nx × R nu → R with 0 a regular value:

S = {(x, u) ∈ R nx × R nu : h(x, u) ≥ 0}.
Then, h is an integral control barrier function (I-CBF) if for any (x, u) ∈ R nx × R nu and t ≥ 0, p(x, u) = 0 ⇒ d(x, u, t) ≤ 0, where the expressions of p(x, u) and d(x, u, t) are given in [START_REF] Ames | Integral control barrier functions for dynamically defined control laws[END_REF].

Analogously to CBFs, I-CBFs are amenable to ensure the forward invariance of subsets of the product space R nx ×R nu . Therefore, they can be used to ensure that the SOC of robot batteries never falls below a desired minimum threshold. In other words, we would like the set of values of SOC, e, above a desired minimum threshold to be forward invariant, which is equivalent to saying that we would like e to be always greater than or equal to a desired minimum threshold. With this above definition of I-CBFs, Theorem 1 in [START_REF] Ames | Integral control barrier functions for dynamically defined control laws[END_REF] gives sufficient conditions for the forward invariance of sets defined through an I-CBF.

In order to apply I-CBF-based control, the considered system ẋ = f (x, u) has to be dynamically extended as follows [START_REF] Notomista | Long-duration robot autonomy: From control algorithms to robot design[END_REF]:

     ẋ = f (x, u) ė = -ηC -1 α i u 2 u = v, (7)
where now v is the control input, and u becomes part of the state. With these dynamically extended dynamics, the following I-CBF can be used to keep the SOC value above a minimum threshold, denoted by e min :

h x (x, e) = e -e min -α c ( p(x) -p c ) , ( 8
)
where e is the value of the SOC of the battery of robot, p(x) denotes the position of the robot, p c is the location of a charging station, and α c is a monotonically increasing function. The term α c ( p(x) -p c ), then, represents the amount of energy required to reach the charging station. The I-CBF h x has relative degree 2 with respect to the new control input v. Therefore, following the approach in [START_REF] Notomista | Persistification of robotic tasks[END_REF], we can define an following auxiliary I-CBF, h u , which explicitly depends on u, and thus has relative degree 1 with respect to v:

h u (x, e, u) = ḣx (x, e, u) + γ x (h x (x, e)) , (9)
where γ x is a class K function. In the next section, it will be shown how to utilize the defined I-CBF in order to control the SOC of a robot battery while executing an environmental monitoring task through by solving the active sensing problem discussed in the previous section.

IV. PERSISTENT ENVIRONMENTAL MONITORING

Combining the model required for active sensing tasks with the one that describes the dynamics of battery SOC, we can state the full model of a robot with state x, SOC e, control input u, dynamically extended control input v, and measurement model m, moving in an environment characterized by an environmental field parameterized by θ, as follows:

               ẋ = f (x, u) ė = -ηC -1 α i u 2 u = v θ = g(θ, t) y = m(x, θ). (10)
We are now ready to state the optimization-based controller used to realize persistent environmental monitoring, in the form of the following optimization nonlinear program (NLP):

Main NLP (11) minimize v, ẋc,δ v 2 + κδ 2 (11a) subject to ∇ xc V T sens ẋc ≥ φ(x c ) -δ (11b) v = ψ ( ẋc ) (11c) ḣu (x, e, u, v) ≥ -γ u (h u (x, e, u)), (11d)
where the optimization variables are the dynamically extended control input, v, the parameters x c describing the state trajectory x to optimize, and a slack variable δ, whose purpose will be explained in the following.

As discussed in Section II, environmental monitoring seen as an active sensing task can be executed by optimizing the parameters of the state trajectory, x c , to plan a path that maximizes a norm of the CG, which, in this work, is chosen to be the trace. The advantage of the trace over other norms will be clear in the next section and is discussed in more detail in Remark 5. Thus, letting V sens (x c ) = trace (G c (t s , t f ))

2 , and choosing a positive definite function φ, the inequality

∇ xc V T sens ẋc ≥ φ(x c
) results in an increase of the trace of the CG, as desired.

The energy constraints are encoded by the inequality (11d) obtained starting from the auxiliary I-CBF h u defined in [START_REF] Mihaylova | A comparison of decision making criteria and optimization methods for active robotic sensing[END_REF]. The slack variable δ is used to prioritize the objective of keeping the SOC above a minimum threshold over the environmental monitoring task. This way, the latter will be executed as long as the former is satisfied. And, if the energy in the robot battery is depleting, the main active sensing task will be sacrificed just by the amount δ. Moreover, the slack variable δ makes the feasibility of the NLP (11) hinge on constraint (11d) only, whose satisfaction depends on the function h x , as well as the robot and energy dynamics.

Finally, ψ is a function relating the derivative of the trajectory parameters, ẋc , to the dynamically extended input v. In Section V, an example of such a mapping will be given. We conclude this section by stating the main result of the paper which allows for a persistent implementation of environmental monitoring tasks executed through active sensing trajectory optimization.

Proposition 3. If a solution to [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF] exists for all times, then the dynamically extended control input v allows the robot modeled by [START_REF] Chung | On a decentralized active sensing strategy using mobile sensor platforms in a network[END_REF] to execute persistent environmental monitoring.

Proof. Let v , x c , and δ be the solution of [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF]. Then, the following inequality holds: ḣu (x, e, u, v ) ≥ -γ u (h u (x, e, u)).

(

) 12
Theorem 1 in [START_REF] Ames | Integral control barrier functions for dynamically defined control laws[END_REF] ensures that h u (x(t), u(t), e(t)) ≥ 0, ∀t ≥ 0. By the definition of h u in ( 9), a positive h u implies that ḣx (x, e, u) + γ x (h x (x, e)) ≥ 0. Then, by Theorem 3 in [START_REF] Ames | Control barrier functions: Theory and applications[END_REF], h x (x(t), e(t)) ≥ 0, ∀t ≥ 0. Thus, by ( 8), e ≥ e min + α c ( p(x) -p c ). Hence, since α c ( p(x) -p c ) quantifies the energy required to reach the charging station located at p c , at the charging station e ≥ e min . As a result, the SOC never falls below the desired lower threshold e min .

Setting the gradient of the Lagrangian associated with the problem [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF] to zero, it immediately follows that δ ∝ κ -1 . Therefore, as ∇ xc V T sens ẋc = Vsens ≥ φ(x c ) -δ, V sens increases as long as φ(x c ) > δ, which can be rendered arbitrarily small by increasing the parameter κ. This condition corresponds to the execution of the environmental monitoring task by maximizing the information gathered along the state trajectory parameterized by x c .

V. SIMULATIONS

This section shows an application of persistent monitoring in a planar environment where a gas is spreading from a source, whose location has to be identified. The environmental field is the gas concentration modeled by the following diffusive dynamics [START_REF] Crank | The mathematics of diffusion[END_REF]:

m : (x, θ) → θ 4 e -(p 1 (x)-θ 1 ) 2 +(p 2 (x)-θ 2 ) 2 θ 3 (13)
where θ 1 θ 2 T is the source location to be identified, θ 3 and θ 4 determine the spread of the gas concentration, and

[p 1 (x), p 2 (x)] T = p(x)
is the robot position. The dynamics of θ are as follows:

θ = 0 0 K 1 -K 2 t -3 2 T ( 14
)
where K 1 and K 2 have been set to 0.005 and 0.01, respectively. It is worth pointing out that the model [START_REF] He | Autonomous chemical-sensing aerial robot for urban/suburban environmental monitoring[END_REF][START_REF] Martin | Long-term exploration & tours for energy constrained robots with online proprioceptive traversability estimation[END_REF] represents a time-varying gas concentration model in which the unknown parameters θ(t) = const.

The robot executing the persistent environmental monitoring task has been modeled using single integrator dynamics ẋ = u, x = x 1 x 2 T ∈ R 2 . In the energy model ( 6), the function α i has been chosen to be the identity function, i.e. α i (s) = s. Moreover, following the derivation in [START_REF] Notomista | Persistification of robotic tasks using control barrier functions[END_REF], the energy required to reach the charging station located at p c is evaluated as follows:

α c ( p(x) -p c ) = log ( p(x) -p c ) -log(d max ), (15)
where d max is the maximum distance from the charging station at which the robot can recharge its battery. As an example, consider wireless charging technologies: in these settings, d max represents the footprint of the wireless charging station.

The state trajectory to be optimized has been parameterized as a cubic spline using 4 parameters per each dimension of the state. Then, the parameter vector x c is given as follows:

x c = p x1 s p x1 f t x1 s t x1 f p x2 s p x2 f t x2 s t x2 f T (16)
v = ψ( ẋc ) = 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 ẋc (17)
As a result, ( 11) is a convex quadratic program in v, ẋc , and δ, and, as such, it can be efficiently solved online. Algorithm 1 summarizes the control strategy employed to obtain the results presented in the following two subsections.

Remark 4 (Multi-rate estimation and planning). To improve estimation performance, the state estimation using the EKF in Algorithm 1 can be executed at a higher rate compared to the trajectory optimization implemented by solving [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF]. In particular, in the simulations shown in this section, the EKF runs at a rate 10 times faster than the trajectory optimization. Execute input v * solution of [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF] 4:

Collect measurement y 5:

Estimate x using EKF on system (4) 6:

Update trajectory using ẋ c solution of [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF] 7:

Re-parameterize trajectory 8: end while Remark 5 (Trajectory re-parameterization). The parameters x c of the trajectory are updated at each iteration of the control loop using the optimal ẋc , solution of the optimization problem [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF]. Thus, the trajectory evolves over time while the robot is tracking it using the input v = ψ( ẋc ). Letting ∆t be the time horizon used in the definition of G c -i.e., t f = t s +∆t in (2)and dt the duration of the control loop, Step 7 in Algorithm 1 ensures that at, iteration k + 1, the trajectory in the interval [t s -dt, t f -dt] is identical to the trajectory in the interval [t s , t f ] optimized at iteration k. This step is justified by the fact that the trace of G c satisfies the Bellman principle, as observed in [START_REF] Salaris | Online optimal perception-aware trajectory generation[END_REF].

Remark 6 (Trajectory length and input bounds). As the trace of G c defined over the time interval [t s , t f ] is a measure of the information gathered along the state trajectory from x(t s ) to x(t f ), the optimal ẋc may lead to an unbounded increase of the length of trajectory, requiring control inputs v, and therefore u, which exceed the input bounds which exist on a physical robotic platform. Therefore, in practice, enforcing maximum-length constraints on the trajectory may be required.

In the case of the cubic spline trajectory parameterization considered in this section, this is simply achieved by upperbounding the vector of parameters x c .

A. Environmental monitoring

This section shows the simulation results of the execution of environmental monitoring as an active sensing task implemented using (3), and we compare them with an EKF-based estimation without trajectory optimization. Figure 1 shows the trajectories followed by the robot while optimizing the trace of G c (solid line), and without this optimization (dashed line). As a result of this trajectory optimization, Fig. 2 reports the graphs of the estimated parameters θ (Fig. 2a) as well as the trace of G -1 c (Fig. 2b) during the course of the simulation. As can be seen, by optimizing the trajectory followed by the robot, the parameters θ converge to the true ones, depicted as dashed lines in Fig. 2a. Moreover, Fig. 2b clearly shows the advantage brought by the trajectory optimization in terms of the trace of G -1 c -lower values indicate more information collected. The solid line is obtained by updating the trajectory as in (3), while the dashed line corresponds to the case without trajectory optimization.

B. Persistent Environmental monitoring

In this section, energy constraints are included in the scenario considered before. The system dynamics are the same, however, in this case, the trajectory is updated using ẋc solution of (11) rather than by running the gradient ascent algorithm (3). Moreover, a charging station is located at p c = 0, 0 T , where the robot is driven to by the input v solution of [START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF] when in need of energy.

The results are reported in Figures 3 and4. In particular, the former shows the trajectory followed by the robot (left column) and the energy in the battery of the robot (right column) for each charging cycle (corresponding to each row). Figures 4a and4b show the estimated parameters θ and the trace of G -1 c . Compared to Figures 2a and2b, a slight performance detriment can be noticed as a result of the energy constraints that drive the robot away from the optimal trajectory. This corresponds to values of the slack variable δ which significantly deviate from 0 to allow energy constraints to be satisfied.

VI. CONCLUSIONS

In this paper, we presented a control-theoretic method for online trajectory planning suitable to be combined with energy control and execute persistent environmental monitoring tasks. The trajectory optimization and the energy constraints are considered holistically in a single optimization problem which, under suitable design choices, becomes a convex quadratic program, and thus can be efficiently solved in an online fashion. We showed the results of multiple simulated scenarios of persistent environmental monitoring applied to the estimation of a gas concentration field that follows diffusive dynamics. In order to illustrate the effectiveness of the proposed approach, we compared it to the case where trajectory optimization or energy constraints are not considered. 4b, respectively. In Fig. 4a, ground truth values are depicted as dashed lines. For reference, in Fig. 4b, the performance of the system without energy constraints-shown in Fig. 2b as a solid line-is depicted as a dashed line.

Algorithm 1 environmental monitoring 1 :

 11 Persistent Initialize trajectory parameters x c 2: while true do 3:

Fig. 1 .

 1 Fig. 1. Trajectory of the robot executing environmental monitoring obtained by optimizing the state trajectory (solid) and without optimization (dashed).

  With (solid) and without (dashed) trajectory optimization.

Fig. 2 .

 2 Fig.2. Environmental monitoring simulation. The estimation of the field parameters θ and the trace of the inverse of the CG are reported in Figures2a and 2b, respectively. In Fig.2b, lower values correspond to more information collected.

Fig. 3 .Fig. 4 .

 34 Fig.3. Trajectory (left column) and energy (right column) of the robot executing persistent environmental monitoring by optimizing the state trajectory and tracking it using the solution of[START_REF] Ny | On trajectory optimization for active sensing in gaussian process models[END_REF]. A charging station located at [0, 0] T is depicted as a gray circle in the figures on the left column. The first, second, and third discharging-charging cycles are shown in Figures3a-3b, 3c-3d, and 3e-3f, respectively.

This works was sponsored by the ANR project ANR-20-CHIA-0017 "MULTISHARED"