Théo Mariotte
email: theo.mariotte@univ-lemans.fr

Anthony Larcher

Silvio Montrésor

Jean-Hugh Thomas

Microphone Array Channel Combination Algorithms for Overlapped Speech Detection

Keywords: overlapped speech detection, multi-microphone, distant speech, interpretability

Overlapped speech occurs when multiple speakers are simultaneously active. This may lead to severe performance degradation in automatic speech processing systems such as speaker diarization. Overlapped speech detection (OSD) aims at detecting time segments in which several speakers are simultaneously active. Recent deep neural network architectures have shown impressive results in the close-talk scenario. However, performance tends to deteriorate in the context of distant speech. Microphone arrays are often considered under these conditions to record signals including spatial information. This paper investigates the use of the self-attention channel combinator (SACC) system as a feature extractor for OSD. This model is also extended in the complex space (cSACC) to improve the interpretability of the approach. Results show that distant OSD performance with self-attentive models gets closer to the nearfield condition. A detailed analysis of the cSACC combinationweights is also conducted showing that the self-attention module focuses attention on the speakers' direction.

Introduction

Speaker diarization in the multi-party scenario is still a challenging task [START_REF] Ryant | The Second DIHARD Diarization Challenge: Dataset, Task, and Baselines[END_REF][START_REF] Ryant | The Third DI-HARD Diarization Challenge[END_REF][START_REF] Sell | Diarization is Hard: Some Experiences and Lessons Learned for the JHU Team in the Inaugural DIHARD Challenge[END_REF]. Diarization systems are subject to severe performance degradation when several speakers are overlapping, which may naturally occur in spontaneous speech. Overlapped speech detection (OSD) is thus needed for robust speaker diarization [START_REF] Garcia Perera | Speaker Detection in the Wild: Lessons Learned from JSALT 2019[END_REF] to process overlapped speech segments separately.

Recent advances in neural network based OSD have shown impressive results in near-field conditions compared to classical approaches [START_REF] Geiger | Detecting overlapping speech with long short-term memory recurrent neural networks[END_REF][START_REF] Andrei | Detecting overlapped speech on short timeframes using deep learning[END_REF][START_REF] Sajjan | Leveraging lstm models for overlap detection in multi-party meetings[END_REF][START_REF] Jung | Three-Class Overlapped Speech Detection Using a Convolutional Recurrent Neural Network[END_REF]. Notably, recurrent neural networks such as Long Short-Term Memory (LSTM) [START_REF] Bullock | Overlap-Aware Diarization: Resegmentation Using Neural End-to-End Overlapped Speech Detection[END_REF][START_REF] Bredin | Pyannote.Audio: Neural Building Blocks for Speaker Diarization[END_REF] reached state-of-theart near-field OSD performance. Few researches have however been produced on the distant speech scenario [START_REF] Chen | Multi-channel overlapped speech recognition with location guided speech extraction network[END_REF][START_REF] Cornell | Detecting and Counting Overlapping Speakers in Distant Speech Scenarios[END_REF][START_REF]Overlapped speech detection and speaker counting using distant microphone arrays[END_REF] while this configuration offers practical benefits by avoiding speakers to wear their own microphone.

On distant OSD, Cornell et al. [START_REF] Cornell | Detecting and Counting Overlapping Speakers in Distant Speech Scenarios[END_REF] proposed a Temporal Convolutional Network (TCN) [START_REF] Bai | An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling[END_REF] based OSD and speaker counting architecture. The authors also proposed a detailed benchmark of several OSD architectures [START_REF]Overlapped speech detection and speaker counting using distant microphone arrays[END_REF] and showed that distant OSD can be improved by fusing spatial features with acoustic features. Only handcrafted spatial features were investigated, which may not be as optimal as spatial features learned This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 101007666, the Agency is not responsible for this results or use that may be made of the information and from the French ANR Extensor (ANR-19-CE23-0001-01) in an end-to-end manner.

Spatial information available in multi-microphone signals can be exploited using spatial filtering. Many algorithms based on classical signal processing [START_REF] Benesty | Microphone Array Signal Processing[END_REF][START_REF] Souden | On optimal frequencydomain multichannel linear filtering for noise reduction[END_REF][START_REF] Applebaum | Adaptive arrays[END_REF][START_REF] Benesty | Design of Circular Differential Microphone Arrays[END_REF] or on neural networks [START_REF] Li | Neural network adaptive beamforming for robust multichannel speech recognition[END_REF][START_REF] Minhua | Frequency domain multi-channel acoustic modeling for distant speech recognition[END_REF][START_REF] Park | Robust multichannel speech recognition using frequency aligned network[END_REF] have been proposed in the literature. Signal-based approaches often require extracting explicit information about the input signal (e.g. speaker location, noise statistics). On neural approaches, Gong et al. [START_REF] Gong | Self-Attention Channel Combinator Frontend for End-to-End Multichannel Far-Field Speech Recognition[END_REF] proposed the Self Attention Channel Combinator (SACC) which learns how to combine channels using self-attention [START_REF] Vaswani | Attention is all you need[END_REF]. Combination-weights are learned on the magnitude of the multichannel Short-Time Fourier Transform (STFT) of the input signal. This approach has shown significant improvement in the context of distant Automatic Speech Recognition (ASR).

Building on these previous works, we propose a multichannel OSD system based on the SACC algorithm. SACC acts as a feature extractor from the raw multi-microphone input signal. Two types of OSD systems are implemented, respectively based on Bidirectional LSTM (BLSTM) and TCN sequence modeling. As combination-weights learned by SACC offer limited interpretation, we propose its extension in the complex space, cSACC. cSACC learns complex combination-weights to preserve the phase information in the STFT and get closer to standard spatial filtering algorithms (e.g. beamforming [START_REF] Benesty | Microphone Array Signal Processing[END_REF]). We also show that cSACC weights are easily interpretable and provide information on the spatial directions exploited by the OSD system. To the best of our knowledge, this work is the first application of SACC for the task of OSD and its extension to use complex weights is the first attempt in the domain.

The paper is organized as follows. In Section 2, we introduce the supervised OSD framework. The BLSTM and TCN architectures are also presented. In Section 3, we describe the SACC and the cSACC architectures and their integration into the OSD system. The dataset, experiments and results are presented in Sections 4 and 5. Finally, a detailed analysis of the cSACC combination-weights is conducted in Section 6. Conclusions and perspectives are drawn in Section 7.

Overlapped Speech Detection

Principle

We formulate OSD as a binary classification task. Let X = {X1, X2, . . . , XN } be a sequence of feature vectors, with N being the number of frames, and its associated sequence of binary labels y = {y1, y2, . . . , yN }. We wish to determine the optimal parameters θ of a model f (X, θ) to predict the sequence ŷ = {ŷ1, ŷ2, . . . , ŷN }, with ŷn ∈ {0, 1} being the predicted binary label of the n th frame.

The sequence of feature vectors X is extracted from the raw audio time signal x ∈ R M ×T with M being the number of microphones and T the number of samples. Feature extraction is performed by a function X = g(x) which may be handcrafted (e.g. MFCC, MVDR beamforming) or jointly optimized with the parametric model f in an end-to-end manner (e.g. SACC, cSACC). Similarly to [START_REF] Bredin | Pyannote.Audio: Neural Building Blocks for Speaker Diarization[END_REF][START_REF] Cornell | Detecting and Counting Overlapping Speakers in Distant Speech Scenarios[END_REF][START_REF]Overlapped speech detection and speaker counting using distant microphone arrays[END_REF], the feature extraction applies a downsampling operation to the input data to reduce the computational cost. The overall OSD pipeline is presented in Figure 1.

2.2. Sequence modeling and frame classification 2.2.1. BLSTM Sequence modeling is first performed using the BLSTM architecture as described in [START_REF] Bullock | Overlap-Aware Diarization: Resegmentation Using Neural End-to-End Overlapped Speech Detection[END_REF][START_REF] Bredin | Pyannote.Audio: Neural Building Blocks for Speaker Diarization[END_REF]. Two BLSTM layers composed of P = 256 cells are stacked. The resulting sequence is post processed using a three-layers feed forward network (FFN) with output sizes L1 = 128, L2 = 128 and L3 = 2 respectively. FFN layers are followed by tanh activation functions except for the last one. A softmax activation is applied to the output logits to compute classification probabilities.

TCN

TCN architecture has been proposed as an alternative to BLSTM and FFN for OSD [START_REF] Cornell | Detecting and Counting Overlapping Speakers in Distant Speech Scenarios[END_REF]. This architecture is composed of causal convolutional layers with residual connections. The key feature of the TCN is the dilated convolution, allowing to learn a large temporal context. We employed the same TCN architecture as [START_REF] Cornell | Detecting and Counting Overlapping Speakers in Distant Speech Scenarios[END_REF] composed of R = 5 residual convolutional blocks repeated P = 3 times 1 . Classification is performed by a 1-d convolutional layer followed by a softmax activation function to compute classification probabilities. 1 https://github.com/popcornell/OSDC

Multichannel STFT Seq. Modeling Classification th in th out 0 1 Multichannel audio segment p(c = 1|x) Merge channels V K Q Scaled dot prod. Softmax Softmax V K Q Scaled dot prod. Softmax Softmax R MVN I MVN MVN + LogMel + j Modulus Feature extractor cSACC w R w I Xstft X R stf t X I stf t Xatt

Multichannel feature extraction

Combining and weighting channels coming from multiple microphones allows to select spatial information in a multimicrophone input signal (e.g. beamforming [START_REF] Benesty | Microphone Array Signal Processing[END_REF]). In order to automatically combine multiple channels, we apply the SACC method [START_REF] Gong | Self-Attention Channel Combinator Frontend for End-to-End Multichannel Far-Field Speech Recognition[END_REF] as a feature extractor for the OSD task. This method is also extended in the complex space (cSACC) to preserve all of the STFT information during the weight-estimation procedure.

Self-Attention Channel Combinator

Let X stf t ∈ C M ×N ×K be the multichannel STFT of the input signal x, where K is the number of frequency bins. The SACC algorithm computes combination-weights w ∈ R M ×N ×1 applied to each STFT channel before combining them. Those weights are determined from the log-magnitude of the STFT using self-attention [START_REF] Vaswani | Attention is all you need[END_REF]. Let q, k and v be three linear transformations, mapping the input STFT log-magnitude X log to the query and the key Q, K ∈ R M ×N ×D and the value V ∈ R M ×N ×1 . The combination weights are computed as follows:

w = softmax softmax QK T √ D V , (1)
with D being the output dimension of the linear layers mapping the STFT log-magnitude to Q and K. The last softmax activation constrains the weights to be within the interval [0, 1].

Mean and Variance Normalization (MVN) is applied on each frequency bin before feeding the self attention module to reduce the variation range of the input data. The combined STFT magnitude Xatt is finally obtained as the weighted sum of the different channels:

Xatt = M m=1 w ⊙ ∥X stf t ∥. (2)
Since channels combination can lead to a large variation in the data, MVN is applied to the combined STFT for each frequency bin, and converted to the log-mel scale using F = 64 filters as in [START_REF] Gong | Self-Attention Channel Combinator Frontend for End-to-End Multichannel Far-Field Speech Recognition[END_REF]. The process of the SACC algorithm is presented in Figure 1. The use of multiple self-attention heads [START_REF] Vaswani | Attention is all you need[END_REF] was investigated to compute combination-weights but did not bring significant improvement.

The SACC combination-weights can be visualized as a function of the speaker direction of arrival (DOA) [START_REF] Gong | Self-Attention Channel Combinator Frontend for End-to-End Multichannel Far-Field Speech Recognition[END_REF]. However, since weights belong to the real space, the lack of information about the phase limits their interpretation. In order to preserve the phase information of the STFT along the process and to better understand weights learned by SACC, we propose its extension in the complex space.

Complex Self-Attention Channel Combinator

We propose an extension of the SACC algorithm in the complex space (cSACC). This model learns complex weights directly from the incoming multichannel STFT X stf t . We separate the real and imaginary parts X R stf t and X I stf t respectively. Weights are then computed separately on the real part w R and the imaginary part w I using equation [START_REF] Ryant | The Second DIHARD Diarization Challenge: Dataset, Task, and Baselines[END_REF]. Two self-attention modules are thus needed. The complex combined STFT can be written as:

Xatt = M m=1 w R ⊙ X R stf t + j × w I ⊙ X I stf t , (3)
where j = √ -1. MVN is applied to to each frequency bin on both real and imaginary parts of the STFT before computing combination-weights. The magnitude of the complex weighted STFT is then converted to the log-mel scale using F = 64 filters before being fed to the sequence modeling network. Similarly to SACC, MVN is applied before log-mel conversion.

Experimental study 4.1. Dataset

Experiments are conducted on the AMI 2 meeting corpus [START_REF] Mccowan | The ami meeting corpus[END_REF] using two types of audio signals: close-talk speech material recorded with the speakers' headset and distant speech signals recorded by the microphone array (Array 1). It consists of a uniform circular array (UCA) composed of M = 8 omnidirectional microphones placed on the table during meetings. To train and evaluate our models, the AMI corpus is split into training, development and evaluation subsets, containing about 80 h, 10 h and 10 h of human-annotated audio signals respectively. The data partition follows the protocol proposed in [START_REF] Landini | Bayesian HMM clustering of x-vector sequences (VBx) in speaker diarization: Theory, implementation and analysis on standard tasks[END_REF] which guarantees different speakers across subsets. Overlapped speech binary labels are computed from the manual segment annotation provided. Audio signals are sampled at 16 kHz.

Baselines

Both SACC and cSACC feature extractors are compared to four different baselines. (i) Close talk MFCC: MFCC are extracted from close-talk recordings from the AMI headset mix. This model is referred to as the close-talk reference. (ii) Single Distant Microphone (SDM): MFCC are extracted from the signal coming from the first channel of the UCA. (iii) Channels sum: time signals coming from each microphone of the UCA are directly merged in the time domain without any weighting scheme. MFCC are then extracted. (iv) MVDR beamformer: signals are weighted and merged using the MVDR solution proposed in [START_REF] Souden | On optimal frequencydomain multichannel linear filtering for noise reduction[END_REF]. We use the recent MVDR implementation from the torchaudio toolkit [START_REF] Yang | Torchaudio: Building blocks for audio and speech processing[END_REF]. Since we are dealing with signals recorded under real conditions, speech and noise covariance matrices are estimated using the Coherent-to-Diffuse Ratio (CDR). The CDR is estimated under diffuse noise and unknown DOA conditions [START_REF] Schwarz | Coherent-to-Diffuse Power Ratio Estimation for Dereverberation[END_REF] (eq. 25). The STFT magnitude of the MVDR beamformed signal is then converted to mel-scale using F = 64 filters similarly to the SACC and cSACC approaches.

Training and evaluation procedures

OSD models are trained on 245k two-second segments randomly sampled from the AMI train set. Features are extracted on 25 ms sliding window with 10 ms shift. The learning rate is set to lr = 10 -3 since scheduling schemes did not bring any improvement. Binary cross-entropy is used as a training objective with stochastic gradient descent (SGD) optimizer [START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF]. To counteract class imbalance, 50% of the training segments are augmented on-the-fly by summing them to another randomly sampled training segment. Associated labels of each segment are also combined [START_REF] Bullock | Overlap-Aware Diarization: Resegmentation Using Neural End-to-End Overlapped Speech Detection[END_REF][START_REF] Bredin | Pyannote.Audio: Neural Building Blocks for Speaker Diarization[END_REF]. No other data augmentation procedure is applied to make the comparison between monochannel and multichannel approaches easier. Furthermore, multichannel data augmentation (e.g. additional reverberation) was investigated but did not bring significant improvement. This may be due to different locations of the speakers in the real acquisitions and in the room impulse response (RIR) simulations. After training, detection thresholds are tuned on the development set. Model performance is evaluated using F1-score and aver-2 https://groups.inf.ed.ac.uk/ami/corpus/ age precision (AP) [START_REF]Overlapped speech detection and speaker counting using distant microphone arrays[END_REF] and reported on both development and evaluation subset.

Results

Overlapped speech detection performance with the BLSTM architecture are presented in table 1. As expected, the use of a single distant microphone drastically degrades OSD performance with an absolute 23,5% loss on the evaluation F1-score. Direct sum of the channels offers similar performance. MVDR beamforming improves detection by an absolute 22% F1-score gain, probably because it better takes advantage of spatial information. Self-attention based channel combination also reduces the gap between close talk and distant OSD. The SACC and cSACC methods reach about 23% and 19% absolute improvement on the F1-score respectively compared to the SDM configuration. The same behavior is observed on the AP scores with about 19% and 13% improvement in the distant speech scenario. The SACC model slightly outperforms MVDR without requiring to extract explicit information from the multi-microphone input signal (e.g. noise statistics, speakers location). The selfattention mechanism automatically extract useful spatial information for distant OSD from the multichannel audio data. The extension of the SACC approach in the complex space does not improve the results compared to SACC. However, this formulation allows to better interpret the combination-weights learned by the model, as shown in Section 6. Table 2 presents OSD performance with each feature extractor using the TCN architecture. The use of TCN globally improves detection compared to the BLSTM architecture. For example, AP is improved by 7,6% in the close-talk scenario and by 12% in the SDM scenario. Furthermore, SACC and cSACC still improve detection performance compared to SDM, reaching similar performance as MVDR beamforming. The score gap between SACC and cSACC still appears with TCN sequence modeling. The difference in performance with cSACC could be explained by the fact that the real and imaginary parts of the weight are learned separately. To tackle this, complexvalued deep neural networks [START_REF] Lee | Fully complex deep neural network for phase-incorporating monaural source separation[END_REF][START_REF] Choi | Phase-aware speech enhancement with deep complex u-net[END_REF] or learnable analytic filterbank [START_REF] Pariente | Filterbank design for end-to-end speech separation[END_REF] are going to be investigated. Preliminary work has also shown that the use of a different optimizer could improve the performance with cSACC architecture.

Analysis

Hereafter, we conduct an analysis of the combination-weights learned by the cSACC algorithm. The analysis shows that the model is focusing on the speakers' direction when overlapped speech is detected. The response of standard spatial filters, the so-called beampattern, can be computed based on filter weights and the array geometry [START_REF] Benesty | Microphone Array Signal Processing[END_REF][START_REF] Benesty | Design of Circular Differential Microphone Arrays[END_REF]. The magnitude of the beampattern represents the gain applied in every θ azimuth angular direction. Hence, the beampattern computed on a set of cSACC combinationweights wn provides information about the directions in which the attention is focused.

Let ψm = 2π(m -1)/M be the m th microphone angular position and r the radius of the UCA. The beampattern is defined as [START_REF] Benesty | Design of Circular Differential Microphone Arrays[END_REF]:

Bn[wn, θ] = M m=1 w * m,n e j ω cos(θ-ψm) , (4)
with ω = ωr/c and c = 340 m/s being the speed of the sound.

The n subscript represents the frame-index since combinationweights are time-variant. The Time-Averaged Beampattern (TAB) can be computed from [START_REF] Garcia Perera | Speaker Detection in the Wild: Lessons Learned from JSALT 2019[END_REF]. It represents the average steering direction of the cSACC algorithm on a given time window of length NT . The TAB can be formulated as follows:

B = 1 NT N T -1 n=0
Bn[wn(ω, θs), θ].

(5)

Simulated utterances

The analysis of the cSACC combination-weights is performed on simulated data to analyze them in a controlled environment. First, clean utterances are taken from the Librispeech dataset [START_REF] Panayotov | Librispeech: An asr corpus based on public domain audio books[END_REF]. We simulate a L = 5 m × l = 4 m × h = 3 m room with T60 = 0.8 s reverberation time. RIRs are then computed between two sources and an 8-microphone UCA of r = 5 cm radius using the gpuRIR toolkit [START_REF] Diaz-Guerra | gpuRIR: A python library for room impulse response simulation with GPU acceleration[END_REF]. Finally, each clean utterance from each speaker is convolved with its own set of RIRs. Convolved signals from each location are delayed from each other and summed to generate the artificial distant overlapped speech.

cSACC combination-weights analysis

The cSACC combination-weights are analyzed within the best performing TCN based OSD architecture. An example of beampattern magnitude computed on those weights is proposed in Figure 2. The beampattern is presented as a heat-map function of time and DOA. DOA associated with higher amplitude indicates the steering direction of the cSACC algorithm. This figure shows that the model is switching attention between the angular direction of each speaker. Two examples of TAB for two active speakers are presented in Figure 3. The TAB is only computed on time segments where Figure 2: Beampattern magnitude computed from (4) as a function of time and DOA for two active speakers located in different angular locations (---). Frame-wise MVN has been applied as well as a 100 frames moving average for better readability.

the output probability of the model verifies p(c = 1|x) > 0.6. Hence, the main lobes of the TAB informs us on the angular directions where cSACC draws attention when overlapped speech is detected. The utterance-averaged SRP-PHAT [START_REF] Benesty | Microphone Array Signal Processing[END_REF] is also presented as a heatmap. It has been computed on a circular plane of two meters radius and allows to compare the cSACC steering directions to the distribution of the acoustic energy. Figure 3 shows that, over the utterance, cSACC model steers towards the direction of the two active speakers for these two examples. A quantitative study may be conducted to confirm these observations.

Conclusion

In this paper, we investigated the use of self-attention channel combinators as distant overlapped speech detection front-end. We showed that self attention-based algorithms reach state-ofthe-art performance without requiring handcrafted feature extraction. Furthermore, the original self-attention channel combinator (SACC) [START_REF] Gong | Self-Attention Channel Combinator Frontend for End-to-End Multichannel Far-Field Speech Recognition[END_REF] slightly outperforms MVDR beamformer without requiring to estimate speech signal statistics. The SACC model was also extended in the complex space to preserve all of the information in the STFT. Even if this approach does not improve detection performance, it allows a better interpretation of the learned combination-weights. Combinationweight analysis showed that the model seems to draw attention to the angular directions of the active speakers. Further work will be conducted on the use of self-attention to combine channels in a full diarization pipeline in order to evaluate the benefits of this kind of approach on the overall task. Other signal representations such as learnable filterbank [START_REF] Pariente | Filterbank design for end-to-end speech separation[END_REF] or pre-trained models (e.g. WavLM [START_REF] Chen | Wavlm: Large-scale selfsupervised pre-training for full stack speech processing[END_REF]) will be investigated as an alternative to the STFT.

Figure 1 :

 1 Figure 1: Flowchart of the OSD task with the proposed cSACC feature extractor. Sequence modeling is performed with BLSTM or TCN architecture. The bottom part of the figure presents the inference procedure: two detection thresholds are applied to the positive-class output of the model to extract a binary sequence.

Figure 3 :

 3 Figure 3: TAB computed from (5) (-) on two overlapped speech utterances compared to SRP-PHAT energy map and groundtruth (---) speakers location with two simultaneously active speakers. TAB is normalized between [0, 1] for better visualization.

Table 1 :

 1 OSD performance on the AMI meeting corpus for each feature extraction method used as BLSTM architecture frontend. Bold value indicates the best model in the distant speech scenario.

		F1-score (%) AP (%)
	AMI	Dev	Eval	Dev Eval
	Close talk MFCC	67,9	63,1	71,7 63,6
	Single channel MFCC 49,1	39,6	50,6 42,6
	Sum channels MFCC	49,1	31,5	42,4 34,2
	MVDR	62,7	60,1	67,0 59,6
	SACC	64,8	62,4	67,9 61,6
	cSACC	62,2	58,9	65,0 55,8

Table 2 :

 2 OSD performance on the AMI meeting corpus for each feature extraction method used as TCN architecture front-end.

		F1-score (%) AP (%)
	AMI	Dev	Eval	Dev Eval
	Close talk MFCC	71,6	68,1	76,7 71,2
	Single channel MFCC 62,0	48,5	59,4 54,6
	Sum channels MFCC	57,0	46,6	61,3 56,2
	MVDR	68,3	64,2	72,4	66,1
	SACC	68,2	64,2	72,7 64,3
	cSACC	65,0	60,8	68,7 60,1
	6.1. Beampattern


Yves Fotso Fotso
email: fotsofyves@yahoo.fr

Berge Tsanou

Samuel Bowong

Suzanne Touzeau

Frédéric Grognard

Yves Fotso Fotso

Optimal impulsive control of coffee berry borers in a berry age-structured epidemiological model

Keywords: Hypothenemus hampei, Plant-pest interactions, PDE dynamical model, Integrated pest management, Impulsive control theory, Maximum principle, Numerical simulations

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I INTRODUCTION

Coffee berry borer (CBB), Hypothenemus hampei (Coleoptera: Scolytidae) is the most damaging insect pest of coffee worldwide, affecting both the yield and quality of coffee products (Le [START_REF] Pelley | Pests of Coffee. Tropical science series[END_REF][START_REF] Damon | A review of the biology and control of the coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae)[END_REF], and causing more than US$500 million in damages annually [START_REF] Vega | Global project needed to tackle coffee crisis[END_REF]. The CBB life cycle is composed of five distinct stages: egg, larva, pre-pupal, pupal and adult [START_REF] Aristizàbal | Integrated pest management of coffee berry borer: Strategies from Latin America that could be useful for coffee farmers in Hawaii[END_REF]. It starts when an adult (fertilised) female emerges from an infested berry to seek a new host. It selects and bores a hole in a uninfested berry. It lay its eggs in internal galleries and remains inside the coffee berry after oviposition until it dies. After hatching, the larvae feed on the coffee seeds. Young females mate with their male siblings inside the berry. Male CBB do not fly and remain inside the berry [START_REF] Silva | Mating behavior of the coffee berry borer, Hypothenemus hampei (Ferrari) (Coleoptera: Curculionidae: Scolytinae)[END_REF]. Fertilised females emerge from the infested berry, and the loop is closed. The choice of a new host by colonising females depends mainly on the age of the coffee berry. CBB attack immature and mature coffee berries from three months after flowering up to harvest period, with a preference for older berries [START_REF] Rodríguez | A coffee agroecosystem model: II Dynamics of coffee berry borer[END_REF].

In response to CBB attacks, integrated pest management programs for coffee plantation have been developed, which involve several control strategies. Chemical control consists in applying synthetic insecticides designed to kill the CBB during the colonisation of young coffee berries [START_REF] Aristizàbal | Integrated pest management of coffee berry borer: Strategies from Latin America that could be useful for coffee farmers in Hawaii[END_REF]. Biological control is based on natural enemies or products derived from living organisms and aims at maintaining CBB at an acceptable level of harmfulness. The main enemies of CBB are parasitic hymenoptera, certain entomopathogenic fungus and nematodes [START_REF] Aristizàbal | Integrated pest management of coffee berry borer: Strategies from Latin America that could be useful for coffee farmers in Hawaii[END_REF][START_REF] Damon | A review of the biology and control of the coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae)[END_REF]. Cultural practices consist of various activities, such as sanitary harvesting, which is the elimination of residual berries present on the branches of trees and on the ground [START_REF] Damon | A review of the biology and control of the coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae)[END_REF][START_REF] Aristizàbal | Integrated pest management of coffee berry borer: Strategies from Latin America that could be useful for coffee farmers in Hawaii[END_REF]. Finally, trapping methods consist in using attractive traps to capture colonising females during their migration flights [START_REF] Messing | The coffee berry borer (Hypothenemus hampei) invades Hawaii: Preliminary investigations on trap response and alternate hosts[END_REF].

The purpose of this work is to optimise CBB control strategies consisting in spraying a biocontrol agent or bio-insecticide at discrete times. It is based on a dynamical model describing the interactions between CBB and coffee berries during a cropping season, which incorporates a berry age structure to account for CBB preference for older berries [START_REF] Fotso | Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case[END_REF]. To achieve this goal, we formulate and solve an optimal impulsive control problem whose objective is to maximise the profit, while minimising the CBB population for the next cropping season. We previously studied a similar problem in Fotso [START_REF] Fotso | Modelling and optimal strategy to control coffee berry borer[END_REF], but the model did not include the berry age structure and the control was continuous instead of impulsive. The present paper hence proposes two notable improvements compared to our previous work. Impulsive bio-insecticide spraying, in particular, is more relevant for field implementation.

The remainder of the paper has the following structure. In Section II, we formulate the control problem. In Section III, the first order necessary conditions for optimality are established. Finally, results are illustrated by numerical simulations in Section IV.

II THE MODEL AND CONTROL PROBLEM STATEMENT

Coffee berry -CBB interaction model

In this study, we consider the epidemiological model proposed in [START_REF] Fotso | Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case[END_REF] describing the infestation dynamics of coffee berries by CBB. Coffee berries are characterised by their age and their epidemic status: s(t, a) and i(t, a) are respectively the age-specific density of healthy and infested coffee berries at the time t and age a ∈ [0, a † ], where a † is the maximal berry age. The CBB population is divided in two groups, the colonising females denoted by y(t), which correspond to the flying fertilised females looking for their host, and infesting females that are laying eggs inside the berries. Although the sex ratio is largely female-biased [START_REF] Damon | A review of the biology and control of the coffee berry borer, Hypothenemus hampei (Coleoptera: Scolytidae)[END_REF], we assume that there are enough males in each generation to fertilise young females. The infestation dynamics of coffee berries by CBB are described by the following age-structured model:

         ∂ t s(t, a) + ∂ a s(t, a) = -β(a)f (B, y)s(t, a) -µ(a)s(t, a), ∂ t i(t, a) + ∂ a i(t, a) = β(a)f (B, y)i(t, a) -ν(a)i(t, a), ẏ(t) = -εf (B, y) βs(t, .) -µ y y(t) + φz(t), ż(t) = εf (B, y) βs(t, .) -µ z z(t), (1)
completed by the following boundary and initial conditions: s(t, 0) = g(t), i(t, 0) = 0, s(0, a) = s 0 (a), i(0, a) = i 0 (a).

(2)

In this model, new healthy berries are produced at time-dependent rate g(t). The colonising CBB females infest the healthy berries at rate β(a)f (B, y), where β(.) is the berry age-dependent infestation rate and f (., .) is an interaction function which depends on the total healthy berries B(t) = a † 0 s(t, a)da and colonising females y. This creates a berry transfer from healthy s(t, a) to infested berries i(t, a) and simultaneously, from colonising y(t) to infesting females z(t) (with βs(t, .) = a † 0 β(a)s(t, a)da). The scaling parameter ε corresponds to the number of colonising females per infested berry; usually, ε = 1 CBB/berry, since superparasitism is rarely observed in the plantation [START_REF] Rodríguez | A coffee agroecosystem model: II Dynamics of coffee berry borer[END_REF]. Infesting females lay eggs inside the coffee berries, which go through their development cycle until fertilised adult females emerge at rate φ.

Colonising and infesting CBB females undergo mortality at rate µ y and µ z ; healthy and infested berries at the rate µ(a) and ν(a). This model relies on the following assumptions.

Assumptions: Positivity and smoothness of the functions and parameters of system (1) 1. Parameters φ, ε, µ y , µ z and initial conditions y 0 , z 0 are nonnegative. 2. g(.) ∈ L ∞ + (0, ∞), β(.) ∈ L ∞ + (0, a † ); boundary conditions s 0 (.), i 0 (.) ∈ L 1 + (0, a † ) and are bounded. 3. Mortality rates µ(.), ν(.) ∈ L ∞ + (0, a † ), ν(a) ≥ µ(a); moreover, there exists a real number µ > 0 satisfying µ(a) ≥ µ for almost every a ∈ [0, a † ]. 4. Contact function f (., .) is bounded and C 1 -Lipschitz continuous for both arguments; moreover, f (B, y) decreases with B and increases with y, with f (B, 0) = 0 and for all y > 0, lim B→+∞ f (B, y)B is finite.

(L 1 + (I), .
) is the space of nonnegative measurable functions L 1 (I) equiped by the product norm and (L ∞ + (I), . ∞ ) the space of nonnegative functions L ∞ (I) over the set I ⊂ R.

We then include in this model the possibility to apply a bio-insecticide at discrete times in the coffee plantation. The bio-insecticide can for instance be based on the entomopathogenic fungus Beauveria bassiana [START_REF] Greco | Efficacy of Beauveria bassiana applications on coffee berry borer across an elevation gradient in Hawaii[END_REF]. It prevents berry infestation: when a colonising CBB bores its hole in a healthy berry and comes into contact with the bio-insecticide, it gets killed by the biocontrol agent. To implement this control in the model, we introduce a new state variable v(t), which corresponds to the bio-insecticide load in the plantation at time t. Note that variable i(t, a), representing infested berries, does not appear in the remaining equations of system (1) and is not affected by our control. Therefore, one can drop the i(t, a)-equation to control the system dynamics. We then obtain the following controlled system:

                   ∂ t s(t, a) + ∂ a s(t, a) = -β(a)q(v(t))f (B, y)s(t, a) + µ(a)s(t, a), ẏ(t) = -εf (B, y) βs(t, .) -µ y y(t) + φz(t), ż(t) = εq(v(t))f (B, y) βs(t, .) -µ z z(t), v(t) = -γv(t), for t = t n , ∆v(t n ) = v(t + n ) -v(t n ) = h n for t = t n , s(t, 0) = g(t), s(0, a) = s 0 (a), y(0) = y 0 , z(0) = z 0 , v(0 + ) = v 0 .
(3)

The control h n ≥ 0 is applied periodically in the plantation at discrete times t n = nτ , where n ∈ {0, 1, 2, . . . , N f } and τ is the application period, i.e. the time elapsed between two successive applications. The t + n notation depicts the instant just after t n , so v(t + n ) = lim j→0 + v(t n + j) is the bio-insecticide load instantly after the control application. The bio-insecticide load persists in the plantation and decays with rate γ. Its initial condition v(0 + ) is nonnegative. As in [START_REF] Fotso | Modelling and optimal strategy to control coffee berry borer[END_REF], function q(.) : R + -→ R + represents the modulation of the CBB infestation rate by the bio-insecticide as follows: q(v) = 1 -ξv v+k , where ξ ∈ (0, 1) denotes the maximal effectiveness of the bio-insecticide load and the positive constant k the load half saturation.

Optimal control problem

We now formulate the optimal control problem. The main goal of coffee growers is to produce good quality coffee berries, that can be sold at a high price, with low control costs, so as to maximise their profit. They also aim at reducing the CBB population at the end of the cropping season, to limit the risk for the next season. We assume that the harvest takes place at the end of the cropping season, denoted by t f . We also assume that (almost) all berries, and hence (almost) all infesting females that are inside berries, are picked at t f , so that only colonising females (y) remain in the plantation after harvest. Lastly, we assume that infested coffee berries have a negligible monetary value, so that the yield consists only of healthy berries.

We define the class of admissible controls as:

M := {h = (h n ) ∈ R N f +1 : 0 ≤ h n ≤ h max , ∀n ∈ {0, 1, 2, . . . N f }},
where h max represents the upper bound of control. The objective functional is defined for control h = (h n ) n∈{0,1,2,...N f } ∈ M as follows:

J (h) = a † 0 Θ(a)s(t f , a)da - N f n=0 C(h n ) -C y y(t f ), (4)
subject to the impulsive evolution System (3). The first term is the coffee berry yield, where function Θ(.) denotes the price of healthy berries, a bounded, continuous and increasing function of the berry age. The second term represents the control cost, where C(.) is a continuous function of the control h n . The last term is a penalty on the CBB population that remains in the plantation after harvest, weighted by constant C y .

Our purpose is to maximise the objective functional over the admissible class of controls, i.e. seek h = (h n ) n∈{0,1,2,...N f } that belongs to M so that:

J (h ) = max h∈M J (h). (5)

III OPTIMALITY CONDITIONS

We first establish the existence of an optimal solution to problem (5) in Theorem 1 (see proof in Appendix 1.1) and then characterise the solution using the maximum principle method.

Theorem 1:

There exists an optimal control application h = (h n ) n∈{0,1,2,...N f } that belongs to M, which maximises the objective functional J (.) subject to system (3).

Using the framework of Zabczyk (1995) and Yong et Zhang (1992), we derive the optimal control from a combination of the state and the adjoint variables. We determine the adjoint equations by first introducing the sensitivity functions. Let us denote p = (s, y, z, v) and define the solution map: h -→ p(h). The sensitivity functions are defined by the Gâteaux derivatives

(λ s , λ y , λ z , λ v ) = lim →0 + -1 p(h + h) -p(h) , for h ∈ R N f +1
. Consequently, the sensitivity functions λ s (t, a), λ y (t), λ z (t) and λ v (t), corresponding to the state variables s(t, a), y(t), z(t) and v(t), satisfy the following equations:

               ∂ t λ s + ∂ a λ s = -q(v)[β(a)f (B, y)λ s + βsf B (B, y) λ s + βsf y (B, y)λ y ] + q v (v)f (B, y)βsλ v -µλ s , λy = φλ z -εf (B, y) βλ s -ε βs f y (B, y)λ y -ε βs f B (B, y) λ s -µ y λ y , λz = q(v)[εf (B, y) βλ s + ε βs f y (B, y)λ y + ε βs f B (B, y) λ s ] -q v (v)εf (B, y) βs λ v -µ z λ z , (6)
and:

λv = -γλ v for t = t n , λ v (t + n ) = λ v (t n ) + hn for t = t n , n = 0, 1, 2, . . . N f . (7)
Functions f y (., .), f B (., .) and q v (.) represent the partial derivatives of functions f (., .) and q(.) with respect to their arguments y, B and v. The sensitivity functions respect the following initial and boundary conditions:

λ s (t, 0) = 0, λ s (0, a) = 0, λ y (0) = λ z (0) = λ v (0) = 0. (8)
Next, we introduce Φ s (t, a), Φ y (t), Φ z (t) and Φ v (t), the adjoint variables corresponding to the state variables s(t, a), y(t), z(t) and v(t) respectively. The adjoint equations are derived by using the adjoint operator associated with the sensitivity equations (6,7), together with appropriate tranversality and boundary conditions. The adjoint equations are given by:

               ∂ t Φ s + ∂ a Φ s = q(v)[βf (B, y)(Φ s -εΦ z ) + f B (B, y)( βsΦ s -ε βs Φ z )] + [εβf (B, y) + ε βs f B (B, y)] Φ y + µΦ s , Φy = q(v)f y (B, y) [-ε βs Φ z + βsΦ s ] + ε βs f y (B, y)Φ y + µ y Φ y , Φz = -φΦ y + µ z Φ z , Φv = q v (v)f (B, y) [-βsΦ s + ε βs Φ z ] + γΦ v , (9)
and:

Φ v (t + n ) = Φ v (t n ) for n = 0, 1, 2, . . . N f , (10)
with the following tranversality conditions:

Φ s (t f , a) = Θ(a), Φ s (t, a † ) = 0, Φ y (t f ) = -C y , Φ y (t f ) = Φ v (t f ) = 0. ( 11
)
We obtain the following result about the characterisation of the optimal control strategy.

Theorem 2:

For any optimal control application h

= (h n ) n∈{0,1,2,...N f } ∈ M, then the Gâteaux derivative of J (.) is DhJ (h) = N f n=0 [Φ v (t + n ) -C (h n )] hn
, where C (.) represents the derivation of function C(.) with respect to its argument hn and h ∈ V h , with the set

V h := { h ∈ R N f +1 : ∃ ∈ [0, 1]; h + h ∈ M}.
The method used to prove Theorem 2 is the impulsive maximum principle (see proof in Appendix 1.2), which gives the first order necessary optimality conditions (Yong et Zhang, 1992;Zabczyk, 1995), but does not provide an efficient way to compute the optimal control. We now suppose that the control cost is C(h n ) = ηh 2 n , so that the objective functional (4) is quadratic in control. For this particular cost function, Theorem 3 provides a means for computing the optimal control strategy (see proof in Appendix 1.3).

Theorem 3:

There exists an optimal control h = (h n ) n∈{0,1,2,...N f } ∈ M, with corresponding optimal states s (t, a), y (t), z (t) and v (t), maximising the objective functional J (.) defined in (4). Moreover, the characterisation of the optimal control h is given by:

h n = max 0, min Φ v (t + n ) 2η , h max , ∀n ∈ {0, 1, 2, . . . N f }. ( 12
)

IV NUMERICAL SIMULATIONS

We present numerical simulations illustrating the effect of the optimal control strategy on the infestation process of coffee berries by CBB in the plantation. Our numerical approach is based on an extension of the forward-backward sweep method [START_REF] Leander | Optimal control of continuous systems with impulse controls[END_REF], originally proposed to solve the optimal control of continuous ordinary differential equations with impulsive controls. The procedure consists of the following steps: first, the state equations in (3) are solved using the forward semi-implicit finite difference in time and backward difference in age, with an initial guess for the control variable. Second, the adjoint equations in ( 9) are solved by backward semi-implicit finite difference in time and forward difference in age, using the solutions of the state equations. Next, the control is updated with the new values of the state and adjoint solutions given by Theorem 3. The algorithm is repeated until the states and control converge.

In the simulations, a constant berry production rate g(t) ≡ g and an age-independent berry mortality rate µ(a) ≡ µ are chosen. Furthermore, the infestation rate is defined by the following function:

β(a) = β min 0 ≤ a < a β , β min + β a (1 -e -k β (a-a β ) ) a β ≤ a ≤ a † , (13)
its average value being β = 1 a † a † 0 β(a)da. The CBB-berry interaction function f is modelled by f (B, y) = y y+αB+1 . More details on these functions are given in [START_REF] Fotso | Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case[END_REF]. To apply Theorem 3, we use a quadratic control cost C(h n ) = ηh 2 n . Moreover, we choose a sigmoid function to model the price of healthy coffee berries according to berry age:

Θ(a) = Θa m a m + a m Θ , ( 14
)
where Θ is the asymptotic price of healthy coffee berries, a Θ is the age at which berries are at half asymptotic price and m ∈ N is the Hill constant. Note that with a finite age bounded by a † , the maximum price that mature berries can reach is

Θa m † a m † +a m Θ .
Moreover, we use the following initial conditions. At the beginning of the cropping season, flowering has not started yet. We also assume that there are only colonising females, as infesting females were eliminated from the plantation when berries were picked during the preceding harvest. So there are initially neither coffee berries nor infesting females, i.e. s 0 (a) = 0 for all a ∈ [0, a † ] and z(0) = 0. The number of initial colonising females is set at y(0) = 10 4 females. We also assume that the control is absent in the plantation so the initial load is v(0) = 0.

We consider two periods for the bio-insecticide application in the plantation during the cropping season: τ = 7 days (1 week, red curves) in and τ = 14 days (2 weeks, magenta curves) in Figure 1. Panels (a-d) show the dynamics of the state variables, panel (e) the optimal bioinsecticide application; panel (f) represents the effect of bio-insecticide load on the infestation rate, that is σ(v) = 1 -q(v) = ξv k+v .

Figure 1: Simulation of system (3) without control (dashed blue curves) and with weekly (plain red curves) and every 2 weeks (plain magenta curves) optimal bio-insecticide applications. The healthy berry growth without pest is also represented in (a) (dash-dotted black curve).

In Figure 1, trajectories of colonising and infesting females increase, both in the uncontrolled (dashed blue curves) and controlled (red and magenta curves) case. Due to the CBB preference for older coffee berries, this growth starts slowly. In all cases, the bio-insecticide (panel (e)) is applied at its maximal value at the beginning of the cropping season and decreases progressively. Its effect (panel (f)) rapidly increases but remains below its maximum value (ξ = 0.8), which denotes a trade-off between the control effectiveness and its cost. The colonising (panel (b)) and infesting (panel (c)) female trajectories are much lower when the optimal control is applied (red and magenta plain curves) than in the case without control (dashed blue curve). Consequently, the healthy coffee berries (panel (a)) are notably higher with control at the end of cropping season. The control is more effective when applied weekly than every two weeks (Figure 1), both in terms of healthy berries and colonising females at the end of the cropping season. However, this result could change if the control cost function C in the objective functional (4) was modified, for instance if the bio-insecticide cost η was higher or if C included fixed labour costs for each application.

Figure 2: Age distribution of the healthy coffee berries and their price (yield) at end of the simulation (t = t f ) with optimal bio-insecticide application every week (plain red curves) or every 2 weeks (plain magenta curves), without control (dashed blue curves), and without pest (dash-dotted black curves). Cases without control and with various constant age-independent infestation rates β are also represented (blue shaded area delimited by β min and β max = β min + β a , plain blue curve for average value β).

Figure 2 illustrates the evolution of the final berry density and price of coffee berries as functions of the berry age at harvest for bio-insecticide applications every week (red curves) and every 2 weeks (magenta curves). We observe a strong decrease of the mature coffee berries (≥ 90 days) without and with control (left panel) due to the berry age preference of CBB. However, the price distribution of berries increases with berry age but remains considerably lower than in the case without CBB (right panel). Without control (blue curves), one can observe the strong impact of the infestation rate (blue shaded area), as well as the effect of CBB preference for older berries: with the same average value for the infestation rate, the final berry density is higher with (dashed curve) than without (plain curve) age preference. Indeed, in the former case, young berries have a better chance to remain healthy.

V CONCLUSION

In view of the extent of damages caused by CBB in coffee plantations, we aimed in this work at controlling the infestation dynamics of coffee berries by CBB. The control is based on the application of a bio-insecticide at discrete and periodic times, to prevent CBB from infesting healthy berries. The controlled model we used is a semi-discrete system with a continuous berry age structure to represent the CBB preference for older berries. An optimisation problem was formulated in order to maximise the yield of healthy coffee berries, while minimising the cost of control, as well as the remaining CBB population for the next cropping season. We showed the existence of an optimal control and gave its characterisation using the maximum principle. Numerical simulations confirmed that the application of the bio-insecticide effectively controls the CBB and considerably increases the yield and profit at the end of cropping season. Weekly applications were found to be more efficient than applications every two weeks, both in terms of CBB population reduction and profit maximisation.

des Grands Lacs "Mathématiques, Informatique, Biosciences et Géosciences l'Environnement" (AUF, French Embassy in Cameroon, IRD).

Yong J., Zhang P. (1992). Necessary conditions of optimal impulse controls for distributed parameter Bulletin of the Australian Mathematical Society 45(2), 305-326. doi:10.1017/S0004972700030173.

Zabczyk J. (1995). Mathematical control theory. Birkhäuser, Boston.

A APPENDIX

Proof of Theorem 1

Let h m = max{h n , ∀n = 0, 1, 2, . . . N f }, we have the inequality v(t + n ) ≤ v(t n ) + h m . Let w(t) be the solution of the impulsive system

ẇ(t) = -γw(t) for t = t n , ∆w(t n ) = + n ) -w(t n ) = h m , for t = t n , n = 0, 1, 2 . . .
Then, for t ∈ (t + n , t n+1 ), we have from the first equation of system above w(t) = w(t + n )e -γ(t-tn) . Therefore we get w(t + n+1 ) = w(t + n )e -γτ + h m . As e -γτ < 1 , we deduce that the sequence (w(t + n )) n converges to a fixed point, that is w(t + n ) → hm 1-e -γτ . Finally, we obtain the following periodic solution of w(t) = hme -γ(t-tn)

1-e -γτ

≤ hm 1-e -γτ . Therefore, from the comparison principle [START_REF] Lakshmikantham | Theory of Impulsive Differential Equations[END_REF], we get v(t) ≤ w(t) ≤ hm 1-e -γτ . Hence the state v is bounded. Since function q verifies 1 -ξ ≤ q(v) ≤ 1, then by the comparison principle and using the same approach as in [START_REF] Fotso | Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case[END_REF], the subsystem defined by the s(t, a), y(t) and z(t) variables has a unique nonnegative and bounded solution.

Since the state variables of system (3) are bounded, the objective functional J (.) is finite. So it is possible to define d = sup h∈M J (h) and thus, there is a maximising sequence (h j ) j∈N so that the sequence (J (h j )) j∈N converges to d. Since the set M is compact, then there exists a subsequence still denoted by (h j ) j∈N that converges to h = (h n ) n∈{0,1,2,...N f } so that h belongs to M. As the objective functional J (.) is continuous, then it follows that J (h ) = d. This achieves the proof.

Proof of Theorem 2

Let Q := [0, t f ] × [0, a † ]. For any given h ∈ V h , there is > 0 small enough, such that h := h + h ∈ M. Let s (t, a), y (t), z (t), and v (t) be the solution of system (3) corresponding to the control strategy h . Then, the directional derivative of objective functional J (.) is given by

DhJ (h) = lim →0 + -1 [J (h + h) -J (h)] = a † 0 Θ(a)λ s (t f , a)da - N f n=0 C (h n ) hn -C y λ y (t f ). (15)
Multiplying the equations of the sensitivity system (6) by the adjoint state variables Φ s (t, a), Φ y (t), Φ z (t) and Φ v (t) respectively, and also the equations ( 9) by the sensitivities λ v (t, a), λ y (t), λ z (t) and λ v (t) respectively. We derive the following relation:

Q λ s [∂ t Φ s + ∂ a Φ s ] dadt + t f 0 λ y Φy dt + t f 0 λ z Φz dt + t f 0 λ v Φv dt + Q Φ s [∂ t λ s + ∂ a λ s ] dadt + t f 0 Φ y λy dt + t f 0 Φ z λz dt + t f 0 Φ v λv dt = 0. (16)
We have the following relation:

t f 0 v λv dt = N f -1 n=0 tn+1 tn Φ v λv dt + t f t N f Φ v λv dt = N f -1 n=0 [Φ v λ v ] tn+1 t + n - N f -1 n=0 tn+1 tn Φv λ v dt + [Φ v λ v ] t f t + N f - t f t N f Φv λ v dt = N f -1 n=0 Φ v (t n+1 )λ v (t n+1 ) - N f n=0 Φ v (t + n )λ v (t + n ) - t f 0 Φv λ v dt (since Φ v (t f ) = 0) = N f -1 n=0 Φ v (t n+1 )λ v (t n+1 ) - N f n=0 Φ v (t + n )[λ v (t n ) + hn ] - t f 0 Φv λ v dt = - N f n=0 Φ v (t + n ) hn - t f 0 Φv λ v dt (since Φ v (t n ) = Φ v (t + n ) and Φ v (t 0 ) = 0). (17)
Integrating the equation ( 16) and taking into account the initial and boundary conditions of equations ( 6)-( 9) and also the relation ( 17), we get the relation:

a † 0 Θ(a)λ s (t f , a)da -C y λ y (t f ) - N f n=0 Φ v (t + n ) hn = 0. (18)
Therfore,

a † 0 Θ(a)λ s (t f , a)da - N f n=0 C (h n ) hn -C y λ y (t f ) = - N f n=0 C (h n ) hn + N f n=0 Φ v (t + n ) hn . (19)
Hence DhJ (h) =

N f n=0 [Φ v (t + n ) -C (h n )]
hn for all h ∈ V h . This achieves the proof.

Proof of Theorem 3

Let h = (h n ) n∈{0,1,2,...,N f } ∈ M be the optimal impulsive control which maximises the objective functional J (.) and let s (t, a), y (t), z (t) and v (t) be the corresponding state variables. Then, for an arbitrary but fixed h ∈ V h and there is > 0 small enough such that J (h + h) ≤ J (h ). It is derived from Theorem 2 that

DhJ (h ) = N f n=0 [Φ v (t + n ) -2ηh n ] hn ≤ 0.
Since hn is arbitrary for n ∈ {0, 1, 2, . . . , N f }, which implies that

h n = Φ v (t + n ) 2η
. By taking the lower and upper bounds into account, we obtain the optimal impulsive control (12). This achieves the proof.

All parameter values are given in Table 1.

Table 1: Model and control parameter values. Most parameter values are based on biological collected in the literature; more information is available in [START_REF] Fotso | Modelling and optimal strategy to control coffee berry borer[END_REF][START_REF] Fotso | Modelling and control of coffee berry borer infestation[END_REF][START_REF] Fotso | Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case[END_REF]

  .

	Symbol Description	Value
	t f	Duration of a cropping season	250 days
	a †	maximum age of coffee berry	250 days
	g	Production rate of new coffee berries	1200 berries.day -1
	µ	Natural mortality rate of healthy coffee berries	0.002 day -1
	ε	Colonising CBB per berry (scaling factor)	1 female.berry -1
	β(a) Infestation function (13):	day -1
	β min	minimum infestation rate	0.004 day -1
	β a	age-dependent extra infestation rate	0.036 day -1
	k β	infestation coefficient	0.035 day -1
	a β	infestation threshold age	90 days
	α	CBB-berry interaction constant	0.7 female.berry -1
	φ	Emergence rate of new colonising females	2 day -1
	µ y	Natural mortality rate of colonising females	1/20 day -1
	µ z	Natural mortality rate of infesting females	1/27 day -1
	ξ	Bio-insecticide load maximal effectiveness	0.8
	k	Bio-insecticide load half saturation constant	200 g.day -1
	γ	Bio-insecticide decay rate	1/50 day -1
	η	Bio-insecticide cost	0.002 $.g -2
	τ	Bio-insecticide application period	7 or 14 days
	h max Maximal bio-insecticide application	180 g
	Θ(a) Coffee berry price function (14):	$.berry -1
	Θ	coffee berry asymptotic price	0.025 $.berry -1
	a Θ	berry age at half asymptotic price	120 days
	m	Hill constant	7 -
	C y	Cost of remaining colonising females	10 -4 $.female -1

ACKNOWLEDGEMENTS

This work was supported by EPITAG, an Inria associated team part of the LIRIMA (https:// team.inria.fr/epitag/), as well as the Collège doctoral régional de l'Afrique Centrale et


domain: spi phys

text:
Ali Maziz

Cedric Plesse

Caroline Soyer

Eric Cattan

Frederic Vidal

High speed electromechanical response of ionic microactuators

Keywords: conducting IPNs, PEDOT, microactuators, microsystems, fast ionic motion

This paper presents the synthesis and characterization of thin and ultra-fast conducting polymer microactuators which can operate in the open air. Compared to all previous existing electronic conducting polymer based microactuators, this approach deals with the synthesis of robust interpenetrating polymer networks (IPNs) combined with a spincoating technique in order to tune and drastically reduce the thickness of conducting IPN microactuators using a so-called "trilayer" configuration. Patterning of electroactive materials has been performed with existing technologies, such as standard photolithography and dry etching. The smallest air-operating microbeam actuator based on conducting polymer is then described with dimensions as low as 160x30x6 µm 3 . Under electrical stimulation the translations of small ion motion into bending deformations are used as tools to demonstrate that small ion vibrations can still occur at frequency as several hundreds of Hz. Conducting IPN microactuators are then promising candidates to develop new MEMS combining downscaling, softness, low driving voltage, and fast response speed.

INTRODUCTION

Electronic conducting polymers (ECPs) are now considered as a very important class of materials which exhibit interesting electrical and optical properties [START_REF] Paoli | Conductive polymer blends and composites[END_REF]. They have also attracted considerable attention because of potential dimensional changes generated by the ion expulsion/inclusion motions during oxidation and reduction processes [START_REF] Baughman | Conducting Polymer Artificial Muscles[END_REF][START_REF] Smela | Conjugated polymer actuators for biomedical applications[END_REF][START_REF] Otero | Artificial muscles with tactile sensitivity[END_REF][START_REF] Ding | High performance conducting polymer actuators utilizing a tubular geometry and helical wire interconnects[END_REF]. As a result, ECPs might be used as the active material in actuators for interesting applications such as robotics, prosthetics and micro-valves [START_REF] Baughman | Conducting Polymer Artificial Muscles[END_REF]. Actuators designed for operation in open air are usually built as three-layer system with an internal layer composed of a solid polymer electrolyte (SPE) sandwiched between two ECP layers (ECP//SPE//ECP). The relative differential expansion between conducting polymer layers results in bending [START_REF] Baughman | Conducting Polymer Artificial Muscles[END_REF].

In those systems, the diffusion phenomena of ions and solvents into the ECP is a rate-limiting step in the actuation performances as with every ionic EAP [START_REF] Mirfakhrai | Polymer Artificial Muscles[END_REF], making them considered to be slow devices, typically in the range of 1 Hz which is not sufficient for many applications. In order to increase speed, the ionic conductivity has to be increased and/or the distance between electroactive layers has to be reduced, The miniaturization of tri-layer actuators operating in air has been described using a laser ablation technique [START_REF] Alici | Conducting polymer microactuators operating in air[END_REF] or a combination of photolithography and reactive ion etching (RIE) [START_REF] Khaldi | Conducting interpenetrating polymer network sized to fabricate microactuators[END_REF] for fabricating sub-millimeter sized devices. In such system the SPE layer plays a key role in the ultimate performances of these electrochemomechanical devices and the decrease of its thickness is the main limiting factor hindering further miniaturization of the conducting polymer actuators operating in air. Indeed, they have to provide suitable mechanical properties while keeping ionic conductivity as high as possible. This drawback has been recently solved by the synthesis of ultrathin conducting interpenetrating polymer network (IPN) microactuators using the spincoating method and patterning through microsystem techniques. The SPE is synthesized as an IPN based on a PEO network for ionic conductivity with improvement of mechanical properties by the interpenetration of a high molecular mass elastomer, nitrile butadiene rubber (NBR), within the PEO network. After interpenetration of two poly(3,4ethylenedioxythiophene) (PEDOT) electrodes and incorporation of an ionic liquid as electrolyte, large displacements under +/-4V above 900 Hz corresponding to the fundamental eigenfrequency have then been demonstration for a 690x45x12 µm 3 microbeam actuator [START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF].

In the present work, the smallest air-operating microbeams actuator based on ECP will be described (160x30x6 µm 3 ). The synthesis of 6 µm thick SPE will be presented as well as the interpenetration of PEDOT electroactive electrodes. Patterning of this pseudo-trilayer structure will be discussed as a function of the chemical nature of the polymer partners. Finally electromechanical response under electrical stimulation will be presented.

EXPERIMENTAL SECTION

Synthesis of Poly(ethylene oxide)/nitrile butadiene rubber thin film. PEO/NBR 50/50 IPN films are obtained by an in-situ polymerization according to a modified method of previously described procedure [START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF][START_REF] Festin | Robust solid polymer electrolyte for conducting IPN actuators[END_REF][START_REF] Festin | Electro-active Interpenetrating Polymer Networks actuators and strain sensors. Fabrication, position control and sensing properties[END_REF]. First the linear NBR solution is prepared by dissolving NBR in cyclohexanone (1:5 weight ratio). After complete dissolution, DCP (6 wt% vs NBR) as cross-linker was then added to the mixture. PEO network precursors (50 wt% vs NBR) are introduced by adding PEGM monomer (75 wt% of the PEO network) and PEGDM crosslinker (25 wt% of the PEO network). DCPD initiator (3 wt% PEO precursors) was then introduced and the mixture was stirred for 30 min under vacuum. When homogenization was achieved, aqueous PVA solution at 300 g.L -1 was firstly spincoated onto a glass slide (3000 rpm, 3000 rpm.s-1 for 30 s, corresponding to 2.2 µm thick layer) and dried at 80°C. Afterwards, the homogeneous mixture of PEO and NBR was spincoated. The thickness of the polymer films depends on the selected spinner speed, time and solution viscosity. The substrate is then kept at 50°C for 4 h and then postcured for 2 h at 80°C. The temperature is then increased for 60 min at 160°C for the crosslinking step of NBR. The lift-off of the PEO/NBR IPN is carried out in water by dissolving the PVA sacrificial layer. The resulting PEO/NBR IPNs were then dried 4 h at 70°C under vacuum. Freestanding PEO/NBR IPNs are thus obtained.

Trilayer actuators preparation. The methodology for the preparation of conducting IPN actuators has been described previously [START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF]. IPN films are swollen with EDOT vapor under reduced pressure for given lengths of time. The swollen films were then immersed in an iron III trichloride FeCl 3 aqueous oxidative solution (1.5 mol.L -1 ) for an oxidative polymerization of EDOT. This synthetic pathway ensures a non-homogeneous distribution of the electronic conducting polymer throughout the thickness of the PEO/NBR IPN. After polymerization, the resulting films were washed with methanol, i.e. excess of FeCl 3 is removed, and dried under vacuum at 40°C for 1 h. The edges of the thin film were then trimmed off. The chemical synthesis was optimized to avoid short circuit between the two ECP electrodes. Very thin C-IPN films of 3.5 to 14 µm are obtained. Prior to electrochemomechanical characterization, the actuators were immersed in EMITFSI for at least 72 hours until saturation.

Microactuators patterning process. Patterning using standard photolithography and a dry etching technique were carried out according to a procedure described elsewhere [START_REF] Khaldi | Conducting interpenetrating polymer network sized to fabricate microactuators[END_REF][START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF]. Firstly, C-IPN films were fixed on silicon substrates with a PVA adhesive layer (from a spincoated aqueous solution at 1500 rpm/750 rpm.s -1 for 15s, 5 μm thick layer was obtained). The positive photoresist SPR220-7µm (Microchem Corporation) is spincoated at 1000 rpm for 40s, with a ramp rate of 750 rpm.s -1 . The softbake for SPR220 was 115°C on a contact hotplate with slow ramp up and down to avoid formation of cracks. The resulting ~ 17 µm thick photoresist is then patterned. Photoresist was exposed through a mask to an energy dose of 400 mJ.cm -2 during 40 s with a light at 365 nm. A hold time of at least 180 min is required to complete the photoreaction before the next step. The sample is then soaked in the corresponding photoresist developer (Microchem Corporation MF-326) for 10 min in order to obtain a patterned protective layer on top of the C-IPN. In a second step, dry etching of C-IPN was achieved using reactive ion etching. The C-IPN actuators were etched with a mixture of O 2 /CF 4 gas (90/10) at 300 W and 200 mT. Finally, the remaining photoresist is removed with acetone and the C-IPN lift-off was carried out by dissolution of PVA sacrificial layer in water. Ionic liquid incorporation was then performed producing self-standing microactuators.

Electrochemomechanical characterizations. The frequency-dependent tip displacements of the microactuators were measured using a laser displacement sensor (Keyence laser LK-G32). From the resulting displacement signal, strain difference between the electrodes was calculated according to Sugino et al. [START_REF] Sugino | Actuator properties of the complexes composed by carbon nanotube and ionic liquid: The effects of additives[END_REF]:

∆ = (1)
Where Δε is the strain difference, D is half of the peak to peak displacement, h is the thickness of the actuator (measured with SEM microscope after swelling in electrolyte), and L is the distance from the clamped end of the actuator to the projection of the laser.

The output force was measured using a FT-S Microforce Sensing Probes from FemtoTools. Microforce sensors are capable of measuring forces from millinewtons (10 -3 N) down to several nanonewtons (10 -9 N) along the sensors' probe axis. Both compression and tension forces can be measured. The measured force is calculated by the equation F = Sensor Gain * (Vout and Vout is whereas a co calibration da here. The acting as the PEO that can be manipulated many times without suffering from cracking. PEO/NBR SPE can be obtained with tunable thicknesses from 30 µm to 600 nm with very good reproducibility.

Synthes

In a second step conducting IPN were prepared according to the procedure described elsewhere [START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF][START_REF] Festin | Robust solid polymer electrolyte for conducting IPN actuators[END_REF]. Typically, EDOT monomer is incorporated into the SPE by gas-phase swelling under reduced pressure until the swelling ratio reached the desired value of 150% of the initial SPE mass. Next, samples were dipped for different length of time into an aqueous iron III trichloride (FeCl 3 ) solution at 40 °C for the oxidative chemical polymerization of EDOT. Water is chosen as the solvent for the oxidative polymerization of EDOT due to the low EDOT solubility in water, i.e. 1.5 10 -2 mol L -1 . The polymerization rate and the diffusion of EDOT out of the film into the aqueous solution is slow and allows the penetration of FeCl 3 in the matrix before a significant amount of EDOT leaves the matrix. This synthetic pathway ensures a non-homogeneous distribution of the ECP throughout the thickness of the PEO/NBR IPN: two PEDOT layers are then interpenetrated into the both faces of the SPE leading to the desired pseudo-trilayer configuration. Three different thicknesses have been obtained: 3.5, 8 and 14 µm. The electrolyte necessary to the redox process of the conducting polymer is subsequently introduced in the electroactive material through a swelling step. The chosen electrolyte is the ionic liquid N,N'-ethylmethylimidazolium bis(trifluoromethylsulfonyl)imide (EMITFSI). After the swelling step the total thickness of the samples is increasing to 6, 12 and 19 µm.

Then the samples can be macroscopically dimensioned with scissors or microsized using photolithography and reactive ion etching. In both cases microbeam geometries have been chosen since mechanical models are well established.

Determination of the Young modulus of macroscopic bending beam actuators

The Young modulus of the beams has been measured on macroscopic samples according to two approaches. First the Young modulus has been estimated by measuring the resonant frequency of 19µm thick beams with different free lengths L ranging from 2.5 cm to 4 mm (Figure2A and 2B) under +/-3V electrical stimulation. Resonant frequencies ranging from 6.4 to 60 Hz have been obtained respectively and are increasing linearly with 1/L 2 . According to equation ( 2) the Young modulus of the macrobeams has been estimated to be 178 MPa which is in the range of previous values [START_REF] Maziz | Demonstrating kHz Frequency Actuation for Conducting Polymer Microactuators[END_REF].

= √

(2)

Where α n is the resonant frequency mode (1.875 for fundamental resonant frequency) and ρ is the mass density (1700 kg/m 3 ).

The Young modulus has been also determined using the Hooke's law. A displacement is applied to the macrobeam (L=2.5mm, width=0.75mm, thickness=19µm) with a sensor probe and the resulting force is measured. The stiffness of the beam is then used to calculate the Young modulus according to equation (3):

= ( 3
)
Where k is the stiffness of the beam, L the length, w the width and h the thickness. A Young modulus of 170 MPa is obtained. Then the two methods provide similar values for the microbeam modulus.

Patterni

The ultrathin materials are performed w conducting IP compared wi the dry etchin Etching rate of the PEO single network and NBR single network are 1 µm/min and 0.5µm/min respectively. However no etching step can be measured for PEDOT in these conditions as it has been described previously. The PEO/NBR IPN presents an etching rate around 0.4 µm/min, close to the one of the rubber network. This result can be explained by the continuity of the NBR network in the whole SPE, then being the limiting factor of the etching process. Interestingly PEO/NBR/PEDOT conducting IPNs present etching rate higher than those of individual components. Etching rate is ranging between 0.8 and 1.05 µm/min for the three different thicknesses of the electroactive materials. This behavior has already been described previously [START_REF] Khaldi | Conducting interpenetrating polymer network sized to fabricate microactuators[END_REF] and has been attributed to a "self-degradation" process in the presence of PEO and PEDOT.

In these conditions selectivities ranging from 1.15 to 1.5 are obtained, indicating that the Megaposit SPR 220 7.0 photoresist was etched more slowly than the actuator. As a consequence the conducting IPNs can be fully patterned before the protective photoresist layer is fully etched. Then microfabrication of conducting IPN microactuators can be performed with precise control of the geometry for the three different thicknesses.

Actuation of the conducting IPN microactuator

To turn microbeam into microactuator EMITFSI ionic liquid is incorporated through a swelling step of the material. As illustrated on figure 3 the microbeam actuator as small as 160x30x6µm 3 can be obtained which is to the best of our knowledge the smallest air-operating microbeam actuator based on ECP. Microbeam actuators are then stimulated electrically with square wave potential (+/-2V) at frequencies ranging from 1Hz to 500Hz. Interestingly, the actuator deformation remains almost constant (90% of the initial amplitude at 1Hz ) up

CONCLUSION

This paper reports on the fabrication of the smallest air-operating microbeam actuator based on conducting polymer. This result was obtained by the synthesis of ultrathin solid polymer electrolyte using a spin coating technique. The solid polymer electrolyte was prepared with an IPN architecture and was synthesized from poly(ethylene oxide) (PEO) and nitrile butadiene rubber (NBR) networks. The combination of ionic conductivity of the PEO partner along with the mechanical properties of NBR network led to flexible SPE as thin as 3.5µm, suitable with microsystem processes. In a second step, the conducting IPN actuator was prepared by oxidative polymerization of 3,4-ethylenedioxythiophene (EDOT) inside the SPE in order to insure inhomogeneous dispersion of PEDOT throughout the PEO/NBR thickness. To micro-size these materials dry etching was used and led to microbeam actuator as small as 160x30x3.5µm 3 . After incorporating an ionic liquid as electrolyte, the swollen microbeam (160x30x6µm 3 ) was stimulated under square wave potential of +/-2V. Actuation was demonstrated up to 500 Hz indicating that ionic motions are still present in the membrane at such frequency and can induce oscillatory bending movement of this soft ionic microactuator.

  Figure 1 su configuration electrolyte m present good of the materi poor mechan [10-11,13-14 combining a (NBR) will e

Figure 1 :

 1 Figure 1:

  Figure 2: (B) Fund of +/-3V measurem dimensio

Figure 3 :

 3 Figure 3: SEM micrograph of conducting IPN microbeam (160x30x6 µm 3 )

Figure 4

 4 Figure 4: potential.

Table 1

 1

		the resonant freq quency of the be eams
		of 1/L 2 under e electrical stimula ation
		imental setup f for Young mod dulus
		posed by sensor r probe tip. Actu uator
		ng photolithog graphic metho ods. The
	at 300 W and 2 200 mT. The R RIE was
	su ummarizes the e etching rate s of the
	d each individu ual component. . Etching rates	are also
	at will be used	to protect the	actuator locall y during

Table 1 .

 1 Etching rate and selectivity of conducting IPN, IPN and individual components in a O 2 /CF 4 (90/10) plasma, 300W, 200 mTorr.

	Material	Nature	Etching time (min)	Etching rate (µm/min)	Selectivity a
	PEO-NBR-PEDOT 14 µm	conducting IPN	5	1.05	1.5
	PEO-NBR-PEDOT 8 µm	conducting IPN	5	0.85	1.2
	PEO-NBR-PEDOT 3.5 µm	conducting IPN	5	0.8	1.15
	PEO-NBR	IPN	5	0.4	0.6
	PEO	saturated polymer	4	1	1.4
	NBR	unsaturated polymer	4	0.5	0.7
	PEDOT	conjugated polymer	/	n.m. b	n.m b
	a ratio between the etching rate of the material and the etching rate of the Megaposit SPR 220 7.0 photoresist in the O 2 /CF 4 (90/10)
	plasma at 300W, 200 mTorr
	b non-measurable