Dr Basma Abdi
email: basma.abdi@aphp.fr

Keywords: SARS-CoV-2, COVID-19, seroprevalence, HIV, humoral immune response, risk factors

Objectives

We studied the kinetics of SARS-CoV-2 antibodies and identified the associated risk factors with the presence of these antibodies.

Methods

In this longitudinal cohort study, we included all PLWHIV. Samples were tested for IgG targeting the nucleocapsid (N) and all positive samples were tested for IgG targeting the spike (S) and IgA targeting S.

Results

1,901 PLWHIV were included: 64.4% male, median age: 53 years (44-60). At inclusion, 254 (13.4%, 95% CI 11.9, 15.0) had positive IgG anti-N. Levels of IgG anti-N, IgG anti-S and IgA anti-S were 3.95 (standard error: 0.16), 199.4 BAU/mL (21.7) and 3.14 (0.21), respectively. Over one year, levels of IgG anti-N and anti-S decreased significantly (-2.83 p<0.0001 and -94.9 BAU/mL p= 0.0010, respectively), while IgA level increased significantly (+2.97 p=0.0032). Multivariable analysis showed that Sub-Saharan African patients were more likely to have positive IgG anti-N in comparison with other patients (OR:4.78 (95% CI 3.39, 6.73), p<0.0001), while active smoking was a protective factor (OR: 0.57 (95% CI 0.36, 0.90), p=0.0176).

Conclusion

We showed in these PLWHIV a high SARS-CoV-2 seroprevalence. A higher seroprevalence was observed in sub-Saharan Africa patients, probably due to socioeconomic status and a lower seroprevalence was observed in active smokers.

As reported in this journal, a recent study in South Africa suggest that HIV is not a risk factor for moderate or severe COVID-19 disease neither is it a risk factor for mortality [START_REF] Venturas | Comparison of outcomes in HIVpositive and HIV-negative patients with COVID-19[END_REF] However other studies described that HIV could be associated with a similar or a higher risk of acquiring COVID-19 and/or worse outcomes [START_REF] Lee | COVID-19 in People Living with HIV: A Systematic Review and Meta-Analysis[END_REF][START_REF] Brown Lillian | The interplay between HIV and COVID-19: summary of the data and responses to date[END_REF] . To date, no longitudinal studies have been conducted in PLWHIV. Thus, we aimed to determine the SARS-CoV-2 seroprevalence in our PLWHIV and to identify factors potentially associated with COVID-19 infection, and then to evaluate the kinetics of anti-SARS-CoV-2 antibodies.

Our study is a longitudinal prospective cohort conducted between April 2020 and September 2021. All HIV-1 patients followed in the Pitié-Salpétriêre hospital were invited to participate in this study. Residual plasma samples obtained from plasma HIV-RNA viral load measurements were used to perform SARS-COV-2 serologies. IgG was measured using the Abbott Alinity instrument. It is a chemiluminescent microparticle immunoassay for semi-qualitative detection of IgG against nucleoprotein (N) and quantitative detection of IgG against the receptor-binding domain (RBD) of spike (S) protein. IgA against the S1 domain of the S protein was measured using enzyme-linked immunosorbent assays (ELISA, EuroImmun). At inclusion, all plasma samples were screened for IgG anti-N and all samples that were confirmed positive for IgG anti-N were tested to detect IgG anti-S and IgA anti-S. Patients with positive serology were evaluated at 6 and 12 months (M).

Univariable and multivariable logistic regression models were used to assess factors associated with the risk of a positive serology at baseline. We used a multiple imputation approach with Fully Conditional Specification method to fill in missing data. The change from baseline in antibodies levels overtime were compared using mixed models for repeated measures with random effects and unstructured covariance matrix.

A total of 1,901 PLWHIV were enrolled in the study. 64.4% of them were male with median age of 53 years (44-60). Only 57 patients reported previous COVID-19 infection without any complications. Among the participants, 26.6% were active smokers and 38.3% were from sub-Saharan Africa.

At inclusion, 254 patients were seropositive, corresponding to a seroprevalence rate in PLWHIV of 13.4% (95% IC 11,9%, 15%). Median age was 50 years (43-57), 53.5% men, 72.7% Sub-Saharan African, only 37 previous COVID-19 infection and 10.8% were active smokers. The main characteristics of the study population are summarized in table 1.

Among seropositive patients, 88.2% and 64.1% had positive IgG anti-S and IgA anti-S respectively at baseline. The mean levels of antibody concentrations were 3.95 (Standard error (SE) 0.16) for IgG anti-N, 199.4 BAU/mL (SE 28.3) for IgG anti-S and 3.14 (SE 0.21) for IgA anti-S. At M6, 51.9%, 87.3% and 75.4% patients had positive IgG anti-N, IgG anti-S and IgA anti-S respectively. At M12, 35.2%, 87.6%, and 81.2% patients had positive IgG anti-N, IgG anti-S and IgA anti-S respectively. Over one year, levels of IgG anti-N and anti-S decreased significantly (-2.83 (SE 0.20) p<0.0001 and -94.9 BAU/mL (SE 28.3) p<0.0001 respectively), while IgA anti-S level increased significantly (+2.97 (SE 0.95) p<0.0001).

Univariable et multivariable logistic regression analyses were performed to assess independent factors associated with positive serology at baseline (Table 2). Fourteen factors were retained for the multivariable analysis showed that the geographical origin and smoking were independently associated with positive SARS-CoV-2 antibodies. Sub-Saharan African patients were more likely to have positive IgG anti-N in comparison with patients originating from France and other countries (OR: 4.78 [95% CI 3.39;6.73], p<0.0001), while active smoking was a protective factor (OR: 0.57 [95% CI 0.36; 0.90], p=0.0176).

To our knowledge, this is the first study evaluating the seroprevalence and assessing the kinetics of SARS-CoV-2 antibodies during one year in the HIV population. Our findings show a higher seroprevalence of SARS-CoV-2 in PLWHIV in comparison to that reported in general population in France in the same period [START_REF] Fabrice | Antibody status and cumulative incidence of SARS-CoV-2 infection among adults in three regions of France following the first lockdown and associated risk factors: a multicohort study[END_REF] . This result could be explained by social and behavioural determinants of health associated with COVID-19 transmission in different communities especially in PLWHIV. Indeed, we found a higher seroprevalence of SARS-CoV-2 in African Sub-Saharan HIV patients, which may reflect social inequalities in health and healthcare in France for people of sub-Saharan African origin.

We showed that levels of IgG anti-N and IgG anti-S decreased significantly while levels of IgA anti-S increased significantly over one year.

Our results reinforce previous studies of evolution of antibody immunity to SARS-CoV-2 showing the decrease of antibody levels with time 5- of 6 to 9 months [START_REF] Ortega Natalia | Seven-month kinetics of SARS-CoV-2 antibodies and role of pre-existing antibodies to human coronaviruses[END_REF] . Antibody kinetics suggest that HIV patients don't exhibit an efficient immune response in case of virus re-exposure.

However, direct conclusions about protective immunity cannot be made only on the basis of humoral immunity. Other investigations in memory B and T cells are needed.

We showed also that active smoking was associated with a lower rate of IgG anti-N antibodies. Our result is in line with other studies showing that active smoking was associated with a lower rate SARS-CoV-2 antibodies and supporting the role of nicotine as protective for SARS-CoV-2 infection [START_REF] Fabrice | Antibody status and cumulative incidence of SARS-CoV-2 infection among adults in three regions of France following the first lockdown and associated risk factors: a multicohort study[END_REF][START_REF] Konstantinos | Systematic review of the prevalence of current smoking among hospitalized COVID-19 patients in China: could nicotine be a therapeutic option?[END_REF] . Many authors hypothesize that this protective role is associated with the nicotine regulation of angiotensin-converting enzyme-2 receptor expression which is involved in SARS-CoV-2 entry 8 . However, data on whether COVID-19 has a greater incidence in non-smokers is still contradictory and the causal role of tobacco in lung cancer and chronic obstructive pulmonary disease should not encourage smoking to limit the risk of developing COVID-19.

In conclusion, the higher seroprevalence observed in sub-Saharan Africa patients highlights the need of an implementation of health and prevention system taking care of vulnerable people especially PLWHIV. More investigations are needed to understand the association between smoking and COVID-19.

Table 1 : Patients characteristics at baseline

 1

	1901	53 (44-60)	254	(43-57)

Table 2 : Factors associated with the risk of a positive serology SARS-CoV-2 (IgG anti-N) at baseline.

 2 Univariable and multivariable logistic regression models were used to assess factors associated with the risk of a positive serology (IgG anti-N) at baseline. Multiple imputation approach with Fully Conditional Specification method was used to fill in missing data. Analyses were run on each of the 10 data sets, including the imputed values, and the results were combined with Rubin's rules. Variables with a univariate p-value <0.20 were included in the multivariable logistic regression model. The significance level of the p-value in the multivariable model was set at 0.05.

		Yes	115	14 (12.3%)	0.9 (0.51, 1.61)
	Bictegravir	1901 No	1540	196 (12.8%)	1	0.1249	1	0.3779

Acknowledgements

We thank the clinical research team (Aida Zarhrate-Ghoul, Soumeya Hammal).

Transparency declarations

None to declare

Funding

This work was supported by the ANRS-MIE and the EMERGEN consortium.

Authors contribution

Conceptualization, BA., VP., A.G.M., LA. and VC; Methodology, ALN., LA., BA., VP. and A.G.M.; Software, ALN. and LA.; Validation, BA., VP., A.G.M., LA. and VC.; Formal Analysis, ALN. and LA,.; Investigation, SD., RP., KZ., SM., ET., CS., MAV., CK., LS., RT.,MW. and SS. ; Data Curation, BA., ALN., LA.; Writing -Original Draft Preparation, BA., VP., ALN. and LA.; Writing -Review & Editing, all authors reviewed and accepted the final version of the article.; Supervision, VP., A.G.M., VC. and LA..; Project Administration, BA., VP. and YD.; Funding Acquisition, A.G.M., VC. and VP.


Carine Jauberthie
email: carine.jauberthie@lifl.fr

Lilianne Denis-Vidal

Ghislaine Joly-Blanchard
email: ghislaine.joly-blanchard@utc.fr

Identifiability and Estimation of Aircraft Parameters and Delays

Keywords: Nonlinear systems, Time-delay systems, Identifiability, Differential algebra, Estimation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

System identification based on physical laws often involves parameter estimation. Before performing estimation problem, it is necessary to investigate its identifiability. It is a mathematical and a priori problem. This paper is concerned with the identifiability of a model given by aerospace domain describing aircraft dynamics. This model is derived from the general equations of motion:

         m dV dt = F e , dC dt = M e , (1)
where dV dt denotes the acceleration of the gravity center, C the kinetic moment, F e the sum of the forces acting on the aircraft and M e the moment of these forces. In order to improve the model accuracy, the input which consists here in the turbulence, is supposed to act at three different points of the fuselage. Let σ i denote the shift operator associated to τ i > 0 and defined for some function v(t) by:

(

σ i v)(t) := v(t -τ i ), i = 1, .., 3, (2)
and let also τ 3 = τ 2τ 1 . The projection of equations [START_REF] Belkoura | Identifiability analysis of linear delay-differential systems[END_REF] on the aerodynamic reference frame of the aircraft yields the following nonlinear and retarded equations [START_REF] Coton | Validation de modèles de représentation du comportement des aérodynes en rafales verticales[END_REF]:

p                        V = k 0 sin(θ -α) + V 2 (k 1 + k 2 α), V ( θ -α) = k 0 cos(θ -α) + k 3 q V + V 2 [k 4 + p 1 (α + u) + p 2 (α + σ 1 u)+ p 3 (α + σ 2 u) + p 3 p 4 (σ 3 α + σ 2 u)], q = k 5 q V + V 2 [k 6 + p 5 (α + u)+ p 6 (α + σ 1 u) + p 7 (α + σ 2 u)+ p 4 p 7 (σ 3 α + σ 2 u)], θ = q, (3)
where V denotes the speed of the aircraft, α the angle of attack, θ the pitch angle, q the pitch rate, and u the input. The set of parameters k i , i = 0, ..., 6 is assumed a priori known (here for instance k 0 = g, the gravity coefficient) and doesn't need to be estimated. The state vector is defined by:

x := (V, α, q, θ) T

and is also assumed available from measurements.

In this paper we are interested in the identifiability of the set:

p = (p, τ 1 , τ 2 ) with p = (p 1 , ..., p 7 ). (5)
The vector p is assumed in some subset Ω of R 9 , and the main problem in our case comes from the introduction of the delays in the set of unknown parameters.

The identifiability of linear delay-differential systems described by:

   ẋ(t) = r i=0 (A i x(t -τ i ) + B i u(t -τ i )) , 0 = τ 0 < τ 1 < . . . < τ r , x(s) = x 0 (s) , s ∈ [-τ r , 0] , (6) 1
has been analyzed by [START_REF] Belkoura | Identifiability analysis of linear delay-differential systems[END_REF], [START_REF] Nakagiri | Unique identification of coefficient matrices, time delays and initial functions of functional differential equations[END_REF], [START_REF] Orlov | On identifiability of linear time delay systems[END_REF] and [START_REF] Verduyn | Parameter identifiability of differential delay equations[END_REF].

In this case, the set of unknown parameters is {A 0 , . . . , A r , B 0 , . . . , B r , τ 1 , . . . , τ r }. To our knowledge, there is no existing method for solving the identifiability problem for general nonlinear and retarded systems with a priori unknown delays. Unlike the free delay case where the properties one can obtain from the linearized system still hold for the original plant in the vicinity of the linearization point, the extension of such result to the retarded case is still an open problem. The aim of this paper is to show how a particular approximation derived from the original nonlinear and retarded system may result in different identifiability conclusions and can be used to estimate the parameters and delays of the model. The presentation is organized as follows: An algebraic approach is developed in Section 2, where the nonlinearities are maintained and a particular input (a practicable input) is used with Padé approximation on the retarded terms of state. In Section 3 the parameters and delays of the model are estimated using this identifiable approximated model. Section 4 presents some numerical results. Most of the computations, concerning identifiability, related to the aircraft models are implemented in Maple, a symbolic computation language and in Matlab for the estimation procedures.

2 Algebraic approach

Input-output approach for testing Identifiability

The input-output approach can be used to test identifiability of some nonlinear systems. This method is based on differential algebra [START_REF] Kolchin | Differential algebra and algebraic groups[END_REF] and consists in rewritting, when it is possible, the nonlinear system as a differential polynomial system that will be completed with ṗ = 0 where p is the vector of parameters of the system. The resulting system Γ p can be described by the following polynomial system:

Γ p    R( ẋ, x, u, u, p) = 0, S(x, y, p) = 0, ṗ = 0. ( 7
)
The notion of identifiability is strongly connected to observability. In the 90s Fliess and Diop propose a new approach of nonlinear observability and identifiability based on differential algebra [START_REF] Fliess | Nonlinear observability, identifiability, and persistent trajectories[END_REF]. The initial conditions are ignored as they are in the model Γ p . A solution of Γ p is a quadruplet of functions (x, y, u, p) which satisfies all the equations of the model. Thus, the solution of these equations may not be unique and some solutions may be degenerate. Therefore the set of non-degenerate solutions x(p, u), ȳ(p, u), corresponding to every possible initial condition, have to be involved in the identifiability definition. Here we adopt the definition introduced in [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF].

Definition 2.1 The model Γ p is globally identifiable with respect to Ω at p if for any p * ∈ Ω, p * = p there exists a control u, such that ȳ(p, u) = ∅ and

ȳ(p, u) ∩ ȳ(p, u) = ∅.
The previous definition is well in line with the usual formulation [START_REF] Walter | Identifiability of state space models[END_REF] which considers initial conditions: If ȳ(p, u) = ȳ(p * , u) then p = p * . In most models there exist atypical points in Ω where the model is unidentifiable. Therefore, the previous definition can also be generically extended so that:

Definition 2.2 Γ p is said to be globally structurally identifiable if it is globally identifiable at all p ∈ Ω except at the points of a subset of zero measure in Ω.

In this approach I is the radical of the differential ideal generated by the equations of Γ p and the ranking:

[p] ≺ [y, u, u] ≺ [x] (8)
is chosen in order to eliminate the state variables.

In general, I should be written as the intersection of regular differential ideals which admit a characteristic presentation. A characteristic presentation [START_REF] Boulier | Computing representations for the radicals of a finitely generated differential ideals[END_REF] is a set of polynomials which is a canonical representant of the ideal and it gives the exhaustive summary of the system which allows us to obtain identifiability results. Generally several characteristic presentations are obtained but only one gives general input-output polynomials and general results of identifiability. The others correspond to particular cases of value of parameters or degenerate solutions. The method is validated by checking the independence of some monomials in y, u and their derivatives occuring in the input-output polynomials.

The computation are achieved with an algorithm implemented in Maple [START_REF] Denis-Vidal | An algorithm to test identifiability of non-linear systems[END_REF]. It is based on the Rosenfeld-Groebner algorithm which has been realized and implemented by F. Boulier in the package Diffalg [START_REF] Boulier | Computing representations for the radicals of a finitely generated differential ideals[END_REF]. Now, when the delays are unknown, the presented approach is unusable unless some modifications. It has been shown in [START_REF] Jauberthie | Aircraft Parameter and Delay Identifiability[END_REF], by using this approach, that among possible approximations some expansion in Taylor series of the retarded terms leads to an unidentifiable model when the use of Padé approximants to approximate the retarded terms gives a global identifiable model. In the next section another idea is considered. To analyze identifiability and also to estimate the parameters and delays of the model ( 3) only an approximation of the retarded state variable is done and a specific input is used.

Identifiability analysis

The delay expression in frequency domain is given by an exponential function which is approached with Padé approximants at first order by:

e -sτ3 ≈ (1 - sτ 3 2 )/(1 + sτ 3 2
).

In the time domain, this yields the differential equation:

(σ 3 v)(t) + τ 3 2 (σ 3 v)(t) = v(t) - τ 3 2 v(t). (9)
With this approximation, one has to define a new state variable z 3 (t), which consists in the unknown term σ 3 α. The input u is given by:

u(t) = 5.73e kt l 1 -l 0 (e kt + l 0 )(e kt + l 1 ) , (10)
where l 0 , l 1 , k are known constants. It is clear from (3) that the identifiablity of p (as well as that of p) only depends on the second and third equations.

Let us define z(t) = e kt thus ż(t) = kz(t) and F 1 (x, ẋ), F 2 (x, ẋ) such that:

   V 2 F 1 (x, ẋ) = V ( θ -α) -k 0 cos(θ -α) -k 3 qV -V 2 k 4 , V 2 F 2 (x, ẋ) = q -k 5 qV -V 2 k 6 .
(11) On the other hand:

               u(t) = 5.73(l 1 -l 0 )z (z + l 0 )(z + l 1 ) , σ 1 u(t) = 5.73a 1 (l 1 -l 0 )z (z + a 1 l 0 )(z + a 1 l 1 ) , σ 2 u(t) = 5.73a 2 (l 1 -l 0 )z (z + a 2 l 0 )(z + a 2 l 1 ) , ( 12
)
where a i = e kτi , i = 1, 2 are parameters to identify (replacing τ i , i = 1, 2).

Finally the original problem is rewritten in the following rational form [START_REF] Jauberthie | Aircraft Parameter and Delay Identifiability[END_REF] where the set of polynomials is given by:

D p                    w 1 = p 1 (α + u) + p 2 (α + σ 1 u)+ p 3 (α + σ 2 u) + p 3 p 4 (σ 2 u + z 3 ), w 2 = p 5 (α + u) + p 6 (α + σ 1 u)+ p 7 (α + σ 2 u) + p 4 p 7 (σ 2 u + z 3 ), τ 3 ż3 = 2(-z 3 + α) -τ 3 α, y 1 = w 1 , y 2 = w 2 , y 3 = α, y 4 = α, y 5 = z 3 , ṗ = 0. ( 13
)
The approach described in the preceding section 2.1 gives the general characteristic presentation which contains two input-output polynomials which are too large to be written here. The exhaustive summary of the parameters in the set p is analyzed with the same algorithm which concludes to the structural global identifiability of D p .

3 Parameter and delay estimation

Introduction

The aim of this section is to estimate parameters and delays of the model ( 3) by using a classical least-square objective function given by:

J(p) = N i=1 (z(t i ) -y(t i , p)) R -1 (z(t i ) -y(t i , p)), (14)
where R is the measurement noise covariance matrix given by:

    25.10 -4 0 0 0 0 0.04 0 0 0 0 0.04 0 0 0 0 0.04     .
The vector z represents the data. The measurement noise is assumed to be white gaussian with zero mean. The initial vector p 0 = (p 0 , τ 10 , τ 20 ) is assumed to be known. The cost function is minimized with respect to the unknown parameters and delays. This problem is solved by a Quasi-Newton method which is implemented in the toolbox optimization of Matlab 5.3. It leads to an estimation denoted by pls the soobtained parameter and delay vector.

Estimation results

The following study was conducted in simulation with the input u(t) given by [START_REF] Malek-Zavarei | Time-delay systems, Analysis, Optimization and Applications[END_REF] with l 0 = e k t1 and l 1 = e k t2 where k is equal to 100, t1 = 0.3 seconds and t2 = 0.6 seconds.

The columns of the Table 1 successively give the true parameter vector p, the initial parameter vector p 0 and the estimated parameter vector pls . The last column gives the relative errors as indicated.

Relation between the approximated system and original system

The adequacy between the approximate system described previously and the discretized original system has been verified. The obtained trajectories (figures 5, 6 ,7, 8) with the both systems are very closed: in fact they are indistinguishable relatively to the measure error of all states. The trajectories presented in the figures 5, 6, 7, 8 are obtained by the following way: p = pls , the discretized original system is given by full line and the approximate system described previously is given by dotted line.

Conclusion

This paper presents an approach for the parameter and delay identifiability of aircraft dynamics which consists in a nonlinear and retarded system. In case of unknown delays, the classical methods to test identifiability can't be applied. By using a practicable input and an approximation of the retarded state based on the Padé approximants, a structurally identifiable system is obtained. This last system is used to estimate parameters and delays of the original nonlinear and retarded model. This procedure gives satisfactory results for the delay estimation and some parameter estimation. The accuracy of the others could be improved by using optimal input design.

Figure 2 :

 2 Figure 2: Speed reconstruction with pls .

Figure 3 :

 3 Figure 3: Pitch angle reconstruction with p 0 .

Figure 4 :

 4 Figure 4: Pitch angle reconstruction with pls .

Figure 5 :

 5 Figure 5: Speed reconstruction (approximate system and discretized original system).

Figure 6 :

 6 Figure 6: Angle of attack reconstruction (approximate system and discretized original system).

Figure 7 :

 7 Figure 7: Pitch rate reconstruction (approximate system and discretized original system).

Figure 8 :

 8 Figure 8: Pitch angle reconstruction (approximate system and discretized original system).

Table 1 :

 1 Estimation results.

	Parameter	p	p 0	pls	| pls -p| | p|
	p 1	0.2440	0.2806	0.4009 0.6430
	p 2	4.6260	5.3199	5.0489 0.0914
	p 3	0.8450	0.9717	0.5950 0.2959
	p 4	0.2550	0.2932	0.2324 0.0886
	p 5	0.3720	0.4278	0.4448 0.1957
	p 6	0.2440	0.2806	0.1235 0.4943
	p 7	-3.0360 -3.4914 -2.9239 0.0370
	τ 1	0.0120	0.0138	0.0126 0.0500
	τ 2	0.0440	0.0506	0.0451 0.0250
	The results presented Table 1 show a good improv-
	ment in the estimation of parameters p 2 , p 4 , p 7
	and delays. The estimation of the other parame-
	ters could be improved by using an optimal input
	design which consists in searching a suitable choice
	of experimental conditions. Thus a significant in-
	crease in accuracy of the parameter estimation can
	be obtained.
	The trajectories presented in the figures 1, 2, 3, 4
	are obtained by the following way: p = p (full line)
	and either p = p 0 (dotted line) (figures 1, 3) or
	p = pls (dotted line) (figures 2, 4). They show a
	good trajectory reconstruction and point out the
	efficiency of the approximated method.


A Villemazet

A Durand-Petiteville

V Cadenat

Autonomous Navigation Strategy for Orchards Relying on Sensor-based Nonlinear Model Predictive Control

This paper deals with autonomous navigation through orchards. It proposes a strategy based on a sensorbased model predictive control law coupled with a spiral-based framework allowing to deal in a similar manner with the two main tasks: the in-row navigation and headland maneuvers. The Model-based Predictive Control scheme allows to deal with stability, boundaries of the actuators, and offers a solution to switch between the different tasks. Simulation results based on ROS and Gazebo are presented at the end of the paper.

I. INTRODUCTION

In this paper, we present a navigation strategy relying on sensor-based model predictive control to make a robot autonomously navigate in an orchard. This work belongs to the movement seeking at dealing with the need of doubling the worldwide food production by 2050 [START_REF] Foley | Solutions for a cultivated planet[END_REF] via the mechanization of the agricultural tools [START_REF] Reid | The impact of mechanization on agriculture[END_REF]. More specifically, one focuses on the autonomous navigation challenge present in the mowing, spraying, pruning, and harvesting processes. To achieve these tasks, the robot has to autonomously drive from the entrance of an alley to its exit, navigate in the headlands to reach the entrance of the next row, and repeat these two steps to cover the whole area of interest. While navigation systems in open-field usually rely on GNSS receivers, this approach may fail in the orchards because of the tree canopy which may block the localization signal [START_REF] Li | Review of research on agricultural vehicle autonomous guidance[END_REF]. For this reason, orchard navigation strategies mainly rely on embedded exteroceptive sensors. Thus, the main approach to drive in the alleys consists in detecting the tree rows, computing geometrical lines in the robot's coordinate frame and using them for guidance. There exists numerous related works and a large spectrum of sensors has been considered, such as vision or Lidar. Some of them are briefly presented hereafter. In [START_REF] Radcliffe | Machine vision for orchard navigation[END_REF] and [START_REF] Opiyo | Medial axis-based machine-vision system for orchard robot navigation[END_REF], a monocular camera is used as the main sensor. In the first one, a sky-based approach allows to compute the vehicle heading, while in the second one, a medial axis-based machine vision is designed to calculate the path to follow. In [START_REF] Durand-Petiteville | Tree detection with low-cost three-dimensional sensors for autonomous navigation in orchards[END_REF], the tree positions are computed by detecting their shadows in a point cloud provided by a RGB-D camera. Finally, in [START_REF] Gai | Using a depth camera for crop row detection and mapping for under-canopy navigation of agricultural robotic vehicle[END_REF], a Time-of-flight camera allows to compute an occupancy grid of the alley. Lidars have also been used. In [START_REF] Blok | Robot navigation in orchards with localization based on particle filter and kalman filter[END_REF], the localization performances of a particle filter and a Kalman filter are compared, while in [START_REF] Danton | Development of a spraying robot for precision agriculture: An edge following approach[END_REF] lines, segments and circles are extracted to detect the tree rows. Both works rely on 2D Lidar data. Despite the variety adrien.durandpetiteville [at]ufpe.br of the approaches, none of the mentioned works deal with the headland navigation problem. Thus, robust, accurate and repeatable turning at the end of a row using relative positioning information with respect to trees is very difficult [START_REF] Vougioukas | Agricultural robotics[END_REF]. A first method proposed in [START_REF] Subramanian | Autonomous vehicle turning in the headlands of citrus groves[END_REF] relies on dead reckoning with slip compensation to perform the u-turn. [START_REF] Bayar | Localization and control of an autonomous orchard vehicle[END_REF] follows a similar approach to execute maneuvers, despite the use of a 2D Lidar, a camera and a tree-detection algorithm to detect the end of the row. Another solution consists in introducing easily distinguishable artificial landmarks at tree rows extremities in order to detect the end of the current row, the entrance of the next one, and localize the robot during turning, using once again dead reckoning [START_REF] Zhang | Mapping orchards for autonomous navigation[END_REF]. Finally, it was proposed to perform the headland maneuvers (u-turn in [START_REF] Durand-Petiteville | Design of a sensor-based controller performing u-turn to navigate in orchards[END_REF] and Ω-shape turn in [START_REF] Flecher | Nonlinear output feedback for autonomous u-turn maneuvers of a robot in orchard headlands[END_REF]) by following a spiral centered on the last tree of the row. Its position is obtained thanks to a 2D Lidar and the spiral is tracked using a non-linear output feedback controller. These approaches thus allow to navigate in the headlands using exteroceptive data, thus increasing the robustness and accuracy of the maneuvers.

In this paper, we aim at continuing the works presented in [START_REF] Durand-Petiteville | Design of a sensor-based controller performing u-turn to navigate in orchards[END_REF] and [START_REF] Flecher | Nonlinear output feedback for autonomous u-turn maneuvers of a robot in orchard headlands[END_REF]. First, we propose to extend the spiral tracking approach to the line following problem. Thus, both in-row and headland navigation rely on the same framework: a spiral tracking problem where the centers of the spirals are computed based on the position of the surrounding trees. Next, the nonlinear output feedback controller is replaced by a Non-linear Model-based Predictive Controller (NMPC) to perform the spiral tracking. NMPC offers two advantages. First, it becomes possible to take into account several constraints. In this work, we include a terminal constraint to deal with the stability of the spiral tracking, which was proven to be only local in [START_REF] Durand-Petiteville | Design of a sensor-based controller performing u-turn to navigate in orchards[END_REF] and [START_REF] Flecher | Nonlinear output feedback for autonomous u-turn maneuvers of a robot in orchard headlands[END_REF]. Moreover, we bound the control inputs to take into account the actuator limitations. In the future, additional constraints could be considered to deal with the presence of obstacles or the field of view of the sensors. Second, a NMPC scheme represents an interesting framework to deal with a sequence of tasks. Finally, we present a task switching method allowing to progressively modify the reference values in the prediction horizon. Thus, it becomes possible to switch from one task to another without creating discontinuities in the control law.

The paper is organized as follows. In the next section, we present the robotic system and the orchard environment. Next, we present the proposed navigation framework and detail the data extraction, the NMPC control and the task switching method. Finally, simulation results are presented.

II. MODELLING

This section is devoted to the problem modelling. It presents the considered robotic system together with the orchard and introduces the required elements to build the navigation strategy. In this work we consider the Agilex Hunter 2.0 car-like robotic platform (see Fig. 1(a)). It is equipped with a 270 • field of view laser range-finder positioned at its front allowing to detect the surroundings. To model this system, let us define the following frames (see Fig. 1(b)): F w = (O w , x w , y w , z w ) fixed to the world, F r = (O r , x r , y r , z r ) attached to the robot, and F l = (O l , x l , y l , z l ) linked to the laser. The robot pose is given by χ(t) = [x(t), y(t), θ(t)] where x(t) and y(t) are the coordinates of O r in F w , whereas θ(t) is the angle between x w and x r . The control vector is defined as U(t) = [υ(t), γ(t)] where υ(t) is the linear velocity along x r and γ(t) the angular position of the steering angle. For such a system, with L the distance between the front and rear wheel axes, the kinematic model is given by:

A. Robot modelling

     ẋ(t) = υ(t) cos(θ(t)) ẏ(t) = υ(t) sin(θ(t)) θ(t) = υ(t) L tan(γ(t)) (1)
B. Orchard description An orchard is shaped by the hand of man. Therefore, it is an environment with a known general structure. Indeed, as shown in Fig. 2, the trees are planted along parallel straight lines, forming rows where the robot can navigate. The rows width and the distance between two consecutive trees are considered to be constant over the whole orchard and are supposed to be roughly known. To navigate through the orchard, the robot has to follow a row, switch to the next one and repeat these two motions until the area is covered. In this work, it is proposed that the first task is performed by following a row central line defined by a set of points, named inner-row points. To realize the row switching, the robot has to perform maneuvers in the headlands (the free spaces in blue on Fig. 2). To do so, it is proposed to define this task using the positions of the trees located at the rows extremity, named pivot trees. Both these latter and the inner-row points will be deduced from the laserrange data, leading to a full sensor-based approach and thus a better execution robustness. These aspects are detailed in the next section.

III. NAVIGATION FRAMEWORK

This section presents the core of our contribution. It successively describes the proposed sensor-based navigation strategy, the extraction of the required information from the sensory data and the design of a suitable controller. To design the orchard navigation strategy, two cases are considered: (i) in-row navigation, i.e., crossing from the entrance to the exit of an alley and (ii) headland navigation, i.e., driving from the current alley to the next one. In order to mathematically define these two tasks, we rely on the spiral method presented in [START_REF] Boyadzhiev | Spirals and conchospirals in the flight of insects[END_REF]. As shown in Fig. 3, the position of a point O s in the robot frame F r can be expressed by its polar coordinates O s (t) = (α(t), d(t)), where α(t) is the oriented angle from the x r vector to the O r O s one, and d(t) the distance between the points O r and O s . It is shown in [START_REF] Boyadzhiev | Spirals and conchospirals in the flight of insects[END_REF] that if both υ(t) and α(t) are constant then O r describes a spiral, and the following equation holds:

A. Approach presentation

ḋ(t) = -υ cos(α) (2)
Thus, [START_REF] Reid | The impact of mechanization on agriculture[END_REF] shows that the type of spiral only depends on parameter α. Indeed, if α ∈ [-π; 0], then O r turns clockwise with respect to O s and counter This analysis allows to design the two previously mentioned tasks. For the first one, i.e., in-row navigation, it is proposed that the robot follows a path defined by a set of inner-row points (see Fig. 4). To do so, we define α * = 0 and d * = 0 where α * and d * are the desired values respectively for α(t) and d(t) for the considered inner-row point. Thus, regulating α(t) and d(t) towards the desired values α * and d * allows positioning the robot at a given position O s . To cross the alley, the robot must therefore perform a sequence of regulations towards a point O s defined for each regulation process by the inner-row point immediately facing the robot. For the second task, i.e., maneuvering in the headland, it is proposed to follow a circle centered on a pivot tree (see Fig. 4). To do so, following the analysis of (2), we impose α * = ±π/2 and d * = r c where r c is the radius of the circle.

-clockwise if α ∈ [0; π]. Moreover, if α ∈] -π; -π 2 [∪] π 2 ; π[, then the spiral is outward and inward if α ∈] -π 2 ; 0[∪]0; π 2 ; π[. It becomes a circle if α = ± π 2 ,
The adopted approach being stated, we now present the different techniques allowing to implement this solution. First, we detail the data processing methods used to compute the points of interest, i.e. the inner-row points and the position of the last trees in the row. Next, we introduce the sensor-based NMPC scheme intended to regulate the current α(t) and d(t) towards their desired values α * and d * . Finally, we present our strategy to deal with the sequence of tasks.

B. Data processing

The proposed navigation strategy relies on the position of the trees in the robot frame. These data are computed from the point-cloud provided by the on-board laser-range finder. The points are first clustered using a distance criterion in order to gather the ones related to the same tree. Next, the position of each tree is calculated by averaging each cluster and expressed using polar coordinates in F r . Now that the positions of the surrounding trees are computed, it is necessary to determine the lines ∆ L and ∆ R passing through the left and right tree rows. To do so, we rely on the normal parametrization proposed in [START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF]:

ρ = d i cos(α i -β) (3)
where (α i , d i ) are the polar coordinates of a tree T i in F r , ρ the algebraic shortest distance from O r to the line and β the angle from x r to the line normal (see Fig. 5).

First we consider that both tree rows are straight and parallel, i.e., their models share the same value for β while ρ has two distinct values. Thus, by varying β from 0 to π on the complete set of positions, the row orientation β row is found when only two values of ρ are obtained. However, the trees are never perfectly aligned and it is then impossible to obtain only two distinct values of ρ. Thus, in order to find the parameters of the two line models, it is necessary to look for two clusters of ρ in place of two distinct values (see Fig. 6(b) and 6(d)). Finally, we have to consider one special case for the line parameters computation: mismatches at the end of the rows when four or less trees are detected. In such cases, there are two values β for which there exists two distinct values of ρ after clustering. We select the value of β for which the inter-clusters distance is the closest one to the row width1 (see Fig. 6(c) and 6(d)).

It is now necessary to compute ∆ M , the central line to ∆ L and ∆ R . ρ m is obtained by averaging ρ L and ρ R and β M is equal to β row to be parallel to both lines. Finally, the innerrow points are obtained by orthogonally projecting the trees on ∆ M (see Fig. 6(a) and 6(c)). Points relative to the same tree rank are gathered when they are too close. (d) Related (ρ, β) plane at the end of the row Fig. 6. Examples of result for data processing. For the two presented cases, i.e. in the row and at its end, the color relies a tree position to its corresponding inner-row point and to its projection on the (ρ, β) plane.

At this step, the positions of the inner-row and of the pivot points have been computed from the laser range finder data.

We now present how those information are used by the NMPC scheme to allow the robot to navigate in the orchard.

C. Nonlinear Model Predictive Control

In order to achieve the two presented tasks, i.e., in-row and headland navigation, we rely on NMPC and we consider the following optimization problem:

U * (k) = min U(k) (J N p (O s (k), U(k))) (4)
with

J N p (O s (k), U(k)) = k+N p -1 ∑ p=k+1 [ Ôs (p) -O * s ] T [ Ôs (p) -O * s ] (5)
subject to

Ôs (p + 1) = f ( Ôs (p), U(p)) (6a) Ôs (k) = O s (k) (6b) C(U * (.)) ≤ 0 (6c)
Thus NMPC consists in computing an optimal control sequence U

* (k) of U(k), with U(k) = [U(k), ..., U(k + N P -1)],
that minimizes the cost function J N p over a prediction horizon of N p steps while taking account the set of user-defined constraints C(U * (k)). The optimal control sequence is of length N c , which represents the control horizon. Thus, if N p > N c , then the last N p -N c predictions are all obtained using the last control input of U * (k).

Cost function:

The cost function is the sum of the quadratic difference between Ôs (.), the predicted polar coordinates ( α, d) of a given point of interest in a set of N p predicted robot frames Fr (.), and O * s the desired coordinates of a point O s in the robot frame. These latter are equal to (0, 0) to make the robot cross the alley, and to (± π 2 , r c ) to drive the robot in the headland. Prediction model: In order to predict the polar coordinates Ôs (p + 1) as shown in (6a), we proceed as follows. First, for a command vector of the given sequence U(k) of N p commands, we compute the two dimensional homogeneous transformation matrices H p|p+1 such as

H p|p+1 = R p|p+1 T p|p+1 0 1×2 1 ( 7
)
where R p|p+1 and T p|p+1 are respectively the rotation matrix and translation vector between Fr (p) and Fr (p + 1), with p ∈ [1, ..., N p ]. They depend on ∆x p|p+1 , ∆y p|p+1 and ∆θ p|p+1 , respectively the displacement along the x r and y r axis, and the rotation around the z r axis between the frames Fr (p) and Fr (p + 1). ∆x p|p+1 , ∆y p|p+1 and ∆θ p|p+1 can be obtained by integrating the kinematic model given by (1) over a sampling time T s and considering a constant velocity vector U(p). Next, we deduce the homogeneous Cartesian coordinates Ôc s (p) from its polar ones Ôs (p). It is then possible to predict the homogeneous Cartesian coordinates of O s in the predicted frame F r (p + 1) as follows:

Ôc s (p + 1) = H -1 p|p+1 Ôc s (p) (8)
Finally, from Ôc s (p + 1) we deduce Ôs (p + 1). This process is initialized using (6b) and is repeated for the N p steps. Zero terminal equality constraint: The stability of the NMPC schema on the finite horizon N p is obtained by adding a zero terminal equality constraint. The respect of this constraint at each iteration guarantees the recursive feasibility, i.e., there exists a trajectory from the current state leading to the desired one. It is defined as the error, in F r , between the Cartesian coordinates ( x(.), ŷ(.)) of a predicted position of the robot φ(.) and the desired one φ ZEC . The zero equality constraint is thus expressed as:

∥ φ(.) -φ ZEC ∥= 0 (9)
In practice, the equality constraint is impossible to achieve and the following inequality constraint is used instead:

∥ φ(.) -φ ZEC ∥-δ ZEC ≤ 0 ( 10
)
where δ ZEC is a user-defined threshold sufficiently small to offer an efficient implementation of the constraint while impacting the optimisation process as the equality constraint. φ ZEC is different for each task. When regulating towards a point, φ ZEC is set to the Cartesian coordinates of O s to force the robot to reach the point. When maneuvering in the headland, φ ZEC is computed as the Cartesian coordinates of a point on the desired spiral. This allows the robot to keep moving forward even after reaching (α * , d * ). Input constraints: The constraints applied to the control inputs allow to take into account the physical limits of the actuators. Moreover, for the set of predictions, they bound the optimization problem. They are defined as:

U(i) -U u U l -U(i) ≤ 0 ( 11
)
where i ∈ [1, N c ], U l and U u the lower and upper boundaries.

D. Task switching

The presented control scheme allows to achieve two tasks: tracking a 'flat spiral' to reach a given position or tracking a 'circle-shape spiral' for the headland navigation. First, to cross an alley, it is necessary to perform position regulations one after the other with respect to a sequence of inner-row points. Second, the robot has also to switch from the sequence of position regulation tasks to the circle-shape tracking one, and vice-versa. Thus, there exists threes cases for switching: (i) from position regulation to position regulation, when the robot crosses the alley, (ii) from position regulation to circleshape tracking, when the robot leaves the alley and enters the headland, and (iii) from circle-shape tracking to position regulation, when the robot leaves the headland and enters the next alley (see Fig. 4). To make the robot switch from one task to another one, it is necessary to update O s , α * and d * . We now detail how the switch is performed for the three previously mentioned cases.

For the first switching case, the robot is initially regulating its position with respect to the first inner-row point in front of it, here referred to as O 1 , and the N p predictions are function of the predicted polar coordinates of O 1 , α * = 0 and d * = 0. While the robot drives towards O 1 , the predicted trajectory gets shorter and the n z last predicted positions, with n z ≥ 1, are equal, or extremely close, to the desired position. In such a case, these n z predictions, which are no more useful to reach the desired position O 1 , are now used to track O 2 , the second inner-row point in front of the robot. In other words, the n z last terms of the cost function are now function of the predicted polar coordinates of O 2 , α * = 0 and d * = 0. If necessary, the cost function is modified at each iteration, up to obtaining a cost function only depending of the coordinates of O 2 . In that case, O 2 is now the main reference and the process is repeated with O 3 , the next inner-point in front of the robot, until the end of the row is reached.

For the second switching case, the robot is first navigating with respect to the last inner-row point of the alley. Similarly to the first case, when the robot begins to be close enough to the point, the n z last predictions are equal, or extremely close, to the desired position. In this case, these predictions are used to track the circle-shape spiral centered on the pivot tree. In other words, the n z last terms of the cost function are now function of the predicted polar coordinates of the pivot tree, α * = ±π/2 and d * = r c . If necessary, the cost function is modified at each iteration, up to obtaining a cost function only depending of the coordinates of the pivot tree.

For the last case, the robot is initially tracking the circleshape spiral with a terminal constraint with respect to a moving point φ ZEC belonging to the spiral. When the first inner-row point of the next alley, here referred to as O 1 , becomes visible, two changes are made. First, r c is modified to make the circle pass through O 1 . Indeed, as the trees are not perfectly aligned and spaced, it is necessary to adjust the circle radius to make it pass through the inner-row point. O 1 is now candidate to become the terminal constraint point φ ZEC . As soon as O 1 is closer to the robot than the current φ ZEC , it becomes the terminal constraint point φ ZEC . This point defines the end of the circle-shape spiral following and it is now possible to rely on a switching strategy similar to the two other cases.

IV. SIMULATION A. Simulation environment and parameters

The simulation robot is a car-like mobile base inspired by the Hunter 2.0 robot (see Fig. 1(a)). Its wheelbase L is set to 0.65 meter, while the limits on the control inputs are defined as follows:

v(t) ∈ [0, v M ] with v M = 1 m/s and γ(t) ∈ [-γ M , γ M ]
with γ M = 0.69 rad. Finally, the embedded laser range finder is located at the robot front, providing the necessary data to detect the surrounding trees. Two examples of data acquisition are shown in Fig. 7 and Fig. 8.

The orchard is composed of 4 rows made of 10 trees. The rows are 20 meters long and 6 meters apart, whereas the space between two successive trees is fixed to 2 meters. A randomisation was applied on the tree positions (about 10 centimeters) and on their orientations (about 360 degrees). The robot has to drive through the first alley, turn on the left to reach the next alley, cross the second alley, turn to its right The presented simulation is implemented with ROS and coupled with Gazebo. The NMPC is implemented using the Python 3.8 language and the cost function is minimised with the SQP solver from the NLopt package [START_REF] Johnson | The nlopt nonlinear-optimization package[END_REF]. Now, concerning the simulation parameters themselves, the sampling period has been fixed to T s = 0.2 second, and both the prediction and control horizons have been taken equal to 12. Moreover, for the parameters of the NLopt, the stopping criteria are the absolute tolerance on the cost function and the maximum optimization time respectively set to 0.05 and 0.16 second, the zero equality constraint threshold δ ZEC = 0.006 and the boundaries of optimization problem are similar to the control inputs ones.

B. Simulation results

The obtained results2 are presented in Fig. 9 where the green vertical dashed lines represent the beginning of an inrow navigation and the red ones the beginning of the maneuver in the headland. In Fig. 9(a) we can see that the robot manages to achieve the navigation task by crossing the three alleys and maneuvering in the two headland sections, driving a total of 73.14 meters in 75.2 seconds. Regarding the control inputs, the linear velocity computed by the optimization process remains at its maximal value throughout the whole navigation (see Fig. 9(b)). This is mostly due to the inclusion of the zero terminal constraint that forces the robot to move forward to reach the desired position. Regarding the steering angle, it tends to 0 rad during the in-row navigation and to 0.2 rad during the headland maneuvers, which is the expected value to describe a circle of approximately 3 meters (width / 2) (Fig. 9(b)). In Fig. 9(c), the distribution of the interest points, i.e., the tasks, is shown. During the in-row navigation, the number of predictions assigned to the first inner-row point decreases from its maximal value to zero. The second innerrow points then becomes the first point of interest and the process is repeated up to the end of the alley. Before entering the headland (before the red dashed lines), we can see that the pivot tree is progressively included as the second point of interest. Once the robot starts the maneuver (red dashed lines), the whole prediction set remains exclusively dedicated to the circle tracking for a large period of time. It is only when the first inner-row point can be calculated (before the green dashed lines) that it is progressively included as the second point of interest. When the robot enters the next alley (green dashed lines) the transition on the prediction set is completed and it only contains the first inner-row point as point of interest. It should be noticed that despite the numerous switches between the tasks, the robot manages to achieve the navigation task. Finally, the processing times given in Fig. 9(d) show that the proposed approach is compatible with a robotic application. Indeed, the total processing time is lower than 100 ms, except for a couple of peaks that could be managed with a time out, offering a 10 hz command rate. Moreover, the solver does not seem to be greatly disturbed by the numerous changes on the optimization when switching tasks.

V. CONCLUSION

We have presented a navigation strategy relying on sensorbased model predictive control to make a robot autonomously navigate in an orchard. The proposed solution extends previous works and offers a sole framework based on spirals to deal the same way with both in-row and headland navigation tasks. It also makes possible the task switching by an adequate update of the spiral parameters and without risking to destabilize the robot. Finally, the design of a sensor-based NMPC instead of more classical non linear control law allows to explicitly deal with several constraints such as stability and actuator saturation. The obtained simulation results show the relevance and the efficiency of the approach.

For future works, we plan to experiment this solution on the Agilex Hunter 2.0 mobile base and test it in real orchards. We aim also at benefiting from the NMPC control structure to integrate new constraints allowing to avoid obstacles, thus performing more complex navigation tasks.

Fig. 1 .

 1 Fig. 1. The robotic system.

( a )

 a Orchard example. (b) General structure of an orchard.

Fig. 2 .

 2 Fig. 2. The orchard.

Fig. 3 .

 3 Fig. 3. Spiral model.

  with a radius equal to d. Finally, if α = 0, then O r follows a straight line towards O s or away from it if α = π.

Fig. 4 .

 4 Fig. 4. Example of navigation in the orchard.

Fig. 5 .

 5 Fig. 5. Normal parametrization of a line passing through the tree T i .

  (a) 3D situation in the row. (b) Related (ρ, β) plane in the row. (c) 3D situation at the end of the row.

Fig. 7 .

 7 Fig. 7. Data acquisition during the in-row navigation.

Fig. 8 .

 8 Fig. 8. Data acquisition during the headland navigation.

  Fig. 9. Simulation results.

Univ. de Toulouse, CNRS, UPS, Toulouse, France

CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France {avillemaze, cadenat}[at]laas.fr

Departamento de Engenharia de Mecânica, Universidade Federal de Pernambuco UFPE, Recife, PE, Brazil

The special case where the row width and tree distance are equal is not taken into account in this work.

A video is also available https://youtu.be/YMKyckTEpSA