domain: chim math

text:
Puneet Mishra
email: puneet.mishra@wur.nl

Jean Michel Roger

Douglas N Rutledge

Alessandra Biancolillo

Federico Marini

Alison Nordon

Delphine Jouan

Rimbaud Bouveresse

Delphine Jouan-Rimbaud-Bouveresse

Datasets for MBA-GUI demonstration Operating procedure and demo analysis Conclusion Validation

Keywords: Chemometrics and Intelligent Laboratory Systems data fusion, multi-sensor, chemometrics, graphical user interface

est destinée au dépôt

Introduction

Sensing technologies play a major role in chemical industries where they are implemented to monitor and optimise process and product properties [1]. Sensing technologies do this by rapid estimation of key critical quality attributes of the process and products. However, sometimes the products or processes are so complex that a single technique fails to obtain sufficient information about the samples. One such case, in the framework of process monitoring applications, is the use of Raman and mid-infrared (MIR) spectroscopy. Both techniques are complementary to one another as they are sensitive to different vibration modes corresponding to molecular groups. Furthermore, they complement one another's drawbacks as Raman signal may get influenced by fluorescence, but it will work with high moisture samples while MIR spectroscopy will be affected by the presence of moisture, but it will work well with fluorescent samples. In such a case, data from both techniques can be utilised in a complementary way: a combination of MIR and Raman spectroscopies could yield better results as each one will compensate for the drawbacks of the other.

Innovations in measurement technologies such as combining multi-spectral techniques for non-destructive estimation of process and product properties is now a major research domain requiring multi-block data analysis techniques. A possible application could be the integration of several different PATs such as MIRS, Raman, near-infrared spectroscopy (NIRS) and fluorescence spectroscopy (FS) into a single measurement probe (Fig. 1A). In the case of a process monitoring application, such a combined probe could be inserted into the process vessel with signals being obtained from multiple sensors, thus enabling the recording of continuous data from multiple techniques (Fig. 1B). However, this will require not only advances in hardware but also in chemometric procedures, such as variable selection methodologies for identification of key spectral regions, and data fusion algorithms for the combination of data coming from multiple techniques.

In some cases, it could be preferable for the data analysis to be performed using a sequential approach so as to highlight the added value coming from including new measurement modes, especially for industrial applications where financial justification is needed for new sensor purchase and integration. Data fusion is a feasible approach to combine all the information J o u r n a l P r e -p r o o f -7]. Data fusion can be dealt with in different ways, depending on the scientific domain. In the machine learning domain, data fusion can be performed at three levels, i.e. low, mid and high-level [8]. Low-level data fusion implies taking all the data together and performing an analysis in a similar way as for data from a single sensor. Midlevel fusion involves some preliminary data refining step where interesting features are extracted from each dataset and subsequently fused together. High-level data fusion is the fusion of the conclusions drawn from the output of the models created on each dataset, using decision rules such as majority voting or averaging. The main aim of data fusion in the machine learning domain is to improve the model accuracy with less attention being paid to the background process. However, in the chemometrics domain, data fusion, or multi-block data analysis, aims not only to improve the model accuracy but also to have a better understanding of the underlying characteristics involved. The aim is to identify the common and the specific hidden factors between and within multiple blocks of measurements, and to subsequently use them to build explainable and explanatory models [9-14].

There are two main tasks that need to be accomplished in multi-block data analysis, i.e., to enhance the data visualisation and improve the predictive performance of models. For data visualisation, a summary and comparison of methods can be found in [9,11,12] In the present work, the first version of a MATLAB-based graphical user interface for multiblock data analysis (MBA-GUI) is presented. The MBA-GUI can perform data visualisation, sequential regression and variable selection on up to 4 different data sources. However, the algorithms implemented in the GUI are not limited to 4 data sources and can be used for any number of data sources. The MBA-GUI can also be used to implement a technique called sequential pre-processing through orthogonalization (SPORT). The MBA-GUI can also be used for a single block scenario. In addition, cases are presented showing how to use the MBA-GUI for data visualisation, regression, classification, variable selection and SPORT in the multi-block scenario. This is the first version of the MBA-GUI and, as multi-block analysis methods progress, the toolbox will be updated to incorporate new algorithms. When looking at the available resources to perform multi-block data analysis, three main toolboxes can actually be found. Data pre-processing is an important step to clean and homogenise the data prior to analysis.

Proper data pre-processing can improve data modelling dramatically. There can be multiple steps in data pre-processing such as smoothing, scatter correction, normalisation and many others [23][24][25]. In the toolbox, we have provided a collection of common pre-processing methods.

Data smoothing reduces high-frequency noise from datasets prior to modelling. In the toolbox, several techniques are provided for performing smoothing in the variable domain. starts by normalising each block, X i , by its Frobenius' norm so that they all have the same total variance. Details regarding the algorithm can be found in [22,35]. Application of the original ComDim (which, for reasons which will become apparent later, we will call ComDim-PCA here) in the case of 2-blocks is presented in Fig. 3. The common components (CCs) are extracted sequentially in an iterative fashion, by extracting the eigenvector associated with the largest eigenvalue from a matrix W, which is the created by concatenating the weighted individual matrices, X i . The weightings (the saliences) are initially all set to 1, but they are recalculated during the iterations to reflect the contribution of each block to the dispersion of the individuals along that CC. After one CC is extracted, the X i matrices are deflated, and the procedure is repeated to obtain the next CC. where the sequential orthogonalization approach was used to fuse the information from different pre-processing techniques. Thus, instead of choosing between pre-treatments, SPORT allows us to make optimal use of the advantages of all pre-treatments. In this GUI, up to four different pre-processing techniques can be used for boosting. Boosting with sequential orthogonalization is a recent approach and has proven to be of high value to improve the prediction accuracies of the models. To perform SPORT, the same data can be loaded in four blocks and then different pre-processing can be applied to them. Further, based on the case of regression/classification, SO-PLS/SO-PLS-DA is used. More details on SPORT can be found in [41].

Apart from multi-block analysis, the MBA-GUI provides an option to perform standard chemometric analyses. For a single block, the MBA-GUI has options to do PCA, PLS as well as CovSel variable selection for regression and discriminant analysis. The one-block analysis will automatically start when the MBA-GUI detects that only one block of data is loaded.

Datasets for MBA-GUI demonstration

Use of the MBA-GUI is demonstrated with three datasets. All the three datasets can be accessed in the same GitHub repository. Operating procedure and demo analysis

The operating procedures are presented in separate sections to demonstrate the use of the MBA-GUI for the analysis of different datasets so that a person with minimal experience can repeat the analysis and use the GUI in their day-to-day multi-block data analysis tasks. The sections are data loading and pre-processing, data visualisation, regression, classification, variable selection and SPORT. All the figures presented in the analysis come directly from the MBA-GUI.

Fig. 5A shows the MBA-GUI interface for loading datasets. Currently, two data formats can be loaded, i.e., .xls and .csv. It is currently possible to load up to 4 blocks of data. The datasets can either come from several sensors or can be the same dataset imported multiple times to perform SPORT fusion. The datasets should be loaded in a logical order since sequential loading of the dataset may improve the performance of sequential methods where the natural order of blocks can provide meaningful insights into the datasets. Once the data are loaded, the figure in the MBA-GUI will be updated to show the data. Fig. 5B shows an example where two data blocks are loaded. Each data block can be separately pre-processed.

The pre-processing can be performed in three steps -smoothing, normalisation or scatter correction, and derivatives. Derivatives are of particular use for NIR data where they can reveal underlying peaks. An example of pre-processing is shown in Fig. 5B. The preprocessing choice in Fig. 5C involved SAVGOL smoothing, SNV normalisation and noderivative. Fig. 5C shows that once the pre-processing is selected the MBA-GUI will show the new pre-processed spectra. The user is free to explore multiple pre-processing methods by visualising how pre-processing affects the spectra. After pre-processing, the data are ready for multi-block analysis. There are five push-button options provided -visualisation, There is also an option to restart the complete analysis by clearing all the previous data and operation logs. A point to note is that if there is only one data block, the MBA-GUI options will take the user to standard chemometric analysis where analyses such as PCA, PLS and variable selection can be performed. Due to limited space, the one block chemometric analysis is not presented in this article.

J o u r n a l P r -p r o o f Once the pre-processing option is selected and the Apply button is pressed, the figures will contain the pre-processed spectra.

J o u r n a l P r e -p r o o f

Multi-block data visualisation is a challenging task where the interest is not only in the extraction of latent variables but also in how different blocks are linked with each other or contribute to the extraction of the global LVs. In the present case, the objective is also to limit the extraction of redundant information from the different blocks. In the MBA-GUI, multiple multi-block data visualisation methods are implemented. The presentation here will be limited to the use of ComDim-PLS, which is a supervised common dimension extraction method requiring a response variable to orient the decomposition of the blocks. It will be applied to the olive oil dataset. block that contributed to the CCs. In this case, the higher saliences for block 2 show that CC1 was dominated by the information from that block, whereas CC2 and CC3 were dominated by the information from block 1. Since they were normalised, each block has a total salience of 1 which means, in the present case, a total salience of two. In the figure, we can see the amount of salience extracted by each block for each CC as well as the total salience extracted by each CC. It is also possible to see that the total amount extracted from the 2 blocks by 3

CCs is about 1.75.

(B). (C). (D). (E). (F). (G). (H).

p r o o f

Multi-block regression is useful when multiple sensors are integrated to improve the predictive performance of the model. In the GUI, two types of multi-block regression are implemented, i.e., SO-PLS and MLR, both based on the ComDim-PLS scores. An example using SO-PLS is presented here. Fig. 7 presents the output of the SO-PLS analysis, where Fig. 8A shows the cross-validation error and Fig. 7B shows the number of LVs selected for each data block, i.e., 10 LVs for block 1 and 8 LVs for block 2. Fig. 7C and 7D show the final regression vector based on the LVs extracted from blocks 1 (Fig. 7C) and 2 (Fig. 7D).

The results from the SO-PLS modelling are presented as the calibration plots in Fig. 7E and 7F. The R 2 for calibration were 0.95 and 0.93, respectively. The calibration and prediction errors were 0.69 and 0.66 % dried matter, respectively. It should be noted that when only one block of data is used, the R 2 is lower and the error is higher (as shown in single block analysis in Fig. 10), showing the benefit of multi-block regression.

J o u r n a l P r e -p r o o f Variable selection is a key step to identify the predictive variables that are most responsible for explaining the response variables. Variable selection can give a better understanding of the important parameters and, in many cases, help in the development of cheap multi-spectral sensor systems. In this work, a demonstration of multi-block variable selection using SO-CovSel is given. The analysis was carried out on the olive fruits dry matter content dataset, where two data blocks in the NIR range are used to predict the dry matter in olives.

Currently, only CovSel variable selection is integrated into the MBA-GUI, however, it can be used for both continuous (regression) and discrete (classification) response variables. A cross-validation option is also provided which supports the selection of key variables. Once the response variables are loaded, the calibrate button can be used and the results will appear in new figures. Fig. 9 shows the outcome of the SO-CovSel analysis where Fig. 9A shows the cross-validation error, Fig. 9B shows the variables selected from block 1 and Fig. 9C shows the variables selected from block 2. A total of 11 variables were selected, 10 from block 1 and only 1 from block 2. The number of selected variables is almost 1/10 of the initial 146 in the two blocks. The calibration and prediction R 2 were 0.93 and 0.91, respectively (Fig. 9D and 9E), and the RMSEC and RMSEP values were 0.78 and 0.77 % dried matter. Although there was a slight decrease in R 2 and a slight increase in RMSEP with the models based on selected variables compared to SO-PLS regression, it should be noted that the model is now much simpler as it includes only 11 variables. Choosing the best pre-processing technique can sometimes be a challenging task such as in the case of spectral data. However, multiple pre-processing techniques can provide complementary information; for example, pre-processing with a derivative can help in revealing the underlying peaks and techniques such as SNV can help to reduce multiplicative effects. In the present MBA-GUI, a newly developed approach called SPORT is also integrated. Using SPORT is completely automated and just requires the user to load the same data in multiple blocks. The user can load the same data up to four times, and therefore, can apply a different combination of pre-processing which the user thinks are the best candidates for the type of data.

In the present case, three-blocks were assigned with three different types of pre-processing.

The first block was pre-processed using SNV, the second block has a combination of SAVGOL smoothing and VSN, and the third block had a combination of SAVGOL and MSC. After performing SPORT, the best pre-processing options were selected highlighting the LVs extracted from each block (Fig. 10A). The block having zero LVs is not useful, and therefore, the pre-processing associated with this block provides no improvement over the pre-processing previously used. In the present case, SNV and VSN pre-processing both had non-zeros LVs, in contrast to MSC pre-processing which had zero LVs, meaning that MSC is not useful in this case. Further, both SNV (5 LVs) and VSN (8 LVs) pre-processing has a significant number of LVs in the final calibration. The model obtained had a R 2 of calibration and prediction of 0.94 and 0.92, respectively. Further, the RMSEC and RMSEP were 0.73 % and 0.72 %, respectively. In summary, SPORT identified the best pre-processing method. In its ensemble, the developed toolbox constitutes a comprehensive software suite to address multi-block data analysis problems, which shows a great potential for attracting practitioners by making their life easier in scenarios that might exhibit particularly high complexities.

Disclaimer

The MBA-GUI is free to use for public as it also involves some algorithms from public sources. Great care has been taken while developing the MBA-GUI, however, the authors do not accept any responsibility or liability.

Fig. 1 :

 1 Fig. 1: Scheme of multispectral fibre system (figure courtesy of Art Photonics GmbH,

Fig. 2 :

 2 Fig. 2: Schematic of the tasks that can be performed with the MBA-GUI.

  Three window-based smoothing techniques, i.e., Savitzky-Golay (SAVGOL), moving average and moving polynomial are provided. Further, two data decomposition and reconstruction techniques, i.e., principal components reconstruction and independent ed with the MBA th the provided. In this toolbox, it is up to the user to choose a technique and decide, based on the model performance, which is best for their data. Preprocessing can also be explored automatically with the SPORT approach. All the spectral smoothing techniques are implemented using the codes explained in[23].Multivariate data, and especially spectral data, suffer from a range of physical and chemical factors leading to baseline, additive and multiplicative effects. Prior to data modelling, it is always recommended to perform correction and normalisation of the data. However, depending on the data, the correction or normalisation method may be different. In the toolbox, the user may select several scatter and spectral normalisation techniques, including detrending [26], offset correction, multiplicative scatter correction [27], spline correction (where spline fitting is used to approximate the baseline which is then subtracted from the signal), asymmetric least-squares (ALS) correction [28], standard normal variate (SNV) [26], variable sorting for normalisation (VSN) [29], probabilistic quotient normalisation (PQN) [30], robust normal variate (RNV) [31], logarithm transform, autoscaling, 1 st derivative (Savitzky-Golay), 2 nd derivative (Savitzky-Golay), min-max, norm, range and max correction. All the correction and normalisation methods were implemented using the codes presented in [23]. It is recommended to perform the smoothing step, if required, before baseline correction and normalisation as these techniques can be affected by high-frequency noise. Detailed understanding of pre-processing methods in chemometrics can be found in [23-25]. Many normalisation and baseline correction techniques can remove effects like baseline shift, additive and multiplicative scatter. Derivatives, however, are also able to reveal underlying peaks. Therefore, the user can choose to perform 1 st and 2 nd derivatives or define the order of J o u r n a l P r e -p r o o f a 3 rd pre-processing step. The algorithm for this operation is based on the codes provided in [23]. Data visualisation is the primary step prior to data analysis to gain an insight into the data [32]. Many factors, such as the presence of outliers, any groupings of samples or any strange patterns, can be detected by visualising the data. Based on the type of data, its visualisation can be very simple or very complex. If the data comprises only a few variables then major patterns can be observed via 1D, 2D and 3D plots. However, if the data consists of many variables, then dimension reduction techniques are used to transform the high-dimensional data into a space defined by interesting properties such as the variance. Once projected into a subspace (dimensionality reduction), the samples can be visualised via 1D, 2D and 3D plots together with plots of the transformed variables. Since spectral data usually contain many variables, ranging from hundreds to thousands depending on the spectral resolution, data transformation is almost always required [33]. In the toolbox, we implement a range of data transformation tools to enhance data visualisation task in such cases. The techniques are mainly based on capturing the major sources of variance in the sample domain. The multivariate data transformation techniques are specific to the multi-block scenario as they can highlight the contributions from the different blocks [9]. ComDim belongs to the family of multi-block methods, which aim to extract the global components that highlight the important dimensions as well as the local (block specific) components [22]. Originally, the ComDim method was developed with the name common components and specific weights analysis (CCSWA) [34]. ComDim extracts the global and local components from multiple blocks of data in a sequential way. Each data block has a specific contribution to each common component which is called its 'Salience'. ComDim J o u r n a l P r e -p r o o f

Fig. 3 :

 3 Fig. 3: The ComDim algorithm applied in the 2-block case ([36]): In the first step of the

  , variable selection and SPORT.

Fig. 5 :

 5 Fig. 5: MBA-GUI interface for data loading. (A) Up to four different data blocks can be

Fig. 6A ,

 6A Fig. 6A, 6B and 6C show the scores when a model with three common components (CCs) is

  Fig. 6G, 6H, 6I show the loadings for each CC presented in two figures each, corresponding to the two-blocks. Such loading plots can be used to understand which variables are of interest. J o u r n a l P r e -p r o o f

Fig. 6 :

 6 Fig. 6: Output from ComDim-PLS performed on the olive oil dataset. (A) Scores from CC1,

FigFig. 7 :Fig. 8 :

 78 Fig. Fig.

Fig. 9 :

 9 Fig. 9: Output results for SO-CovSel analysis performed on the olive fruits dry matter content

Fig. 10 :

 10 Fig.10: Selection of the best pre-processing options and their fusion with the SPORT

  menus are extremely easy to browse and their design is capable of guiding even non-expert users through the sequential steps of multi-block data analysis. They enable not only model training/calibration, but also (when it holds) model optimisation (by two different approaches for cross-validation) and external testing/validation. Furthermore, it is worth mentioning that the presented GUI offers as a valuable add-on an implementation of a recently proposed strategy named SPORT (Sequential Preprocessing through ORThogonalisation) for both the selection of the best pretreatment technique and the extraction of complementary information from the outcomes of multiple preprocessing operations.

  multi-blockThe first one is the multi-block toolbox by the

	from (https://github.com/puneetmishra2/Multi-block.git). In the GitHub repository, the
	standalone toolbox executable files can be downloaded as 'Multi-block_toolbox.zip' and the
	function for running the tools in command line as 'Toolbox.zip'. The dataset demonstrated in
	this article can be downloaded as 'Dataset.zip' from the same GitHub repository. All three
	files are available in the link (10/June/2020). To run the toolbox from command line user
	should use the toolbox folder as the current folder and type T1 on the command line which
	University will start the main GUI interface. The user should put the password: 'welovedata' without of Copenhagen, Denmark
	(http://www.models.life.ku.dk/~courses/MBtoolbox/mbtmain.htm), which, having been last parallel and sequential partial least-squares regression' toolbox by NOFIMA [20]. Both these toolboxes provide command line functionalities (which may be difficult to implement for people unfamiliar with the Matlab environment) and, anyway, consist of a limited number of tools. There is also a basic GUI available for performing multi-block component analysis in the domain of Behaviour research [21]. However, the GUI can only perform multi-block component analysis for data visualisation. Therefore, there exists a need for a GUI which has updated J o u r n a l P r e -p r o o f htm), which, ), which o roaches which w aches whic t analysis and nalysis a ulti-block regre i-block regre x by NOFIMA x by NOFI hich may be di h may be t) and, anyway, nd, anyw e for performin or perform rch [21]. How ch [21]. H isualisat sual comma and click run. Then user can load data and run the analysis. See also supplementary file to have a visual understating on how to download and setup the GUI The GUI supports data format of .csv, .xlsx and .mat. A summary of the functionalities is presented in Fig. 2. In summary, the toolbox has options for loading data, three levels of pre-processing, i.e., smoothing, scatter correction and normalisation, and derivative estimation, as well as multi-block data visualisation, regression, classification, variable selection and SPORT. Multi-block variable selection methods are available for both regression and classification cases. Two main types of regression and classification are available, i.e., sequential orthogonalization [15], and common components and specific weights analysis (CCSWA), aka common dimensions (ComDim) [22]. alysis. See al sis. See nd setup the GU setup the he functionalitie unctiona ding data, three g data, three alisation, and d lisation, an on, classificatio classific thods are availa ds are av f regression regression n [15], and com [15], and c J o u r n a l P r e -p r o o f nsions (C sion
	tools and a complete set of functionalities to perform multi-block data analysis.
	The MBA-GUI was built utilising the application builder in MATLAB version 2018b
	(Natick, MA, USA). The application can be downloaded and installed in MATLAB
	(preferred version 2018b or higher) or can be used as a stand-alone executable or can be run
	through the '.mlapp' files in MATLAB command line. If user does not have MATLAB
	version 2018b or greater then it is recommended to install free MATLAB runtime tool and
	run the app as standalone. All the executable and MATLAB function can be downloaded

revised in 2001, focuses only on the two data fusion approaches which were most popular in that period, i.e, multi-block principal component analysis and multi-block partial least squares regression. The second one is the 'multi-block regression by

  be coded to account for class membership and that, as a further step, linear discriminant analysis is applied to either the predicted response or the concatenated scores block. The function integrated into the GUI is the one freely downloadable at https://www.chem.uniroma1.it/romechemometrics/research/algorithms/.The aim of ComDim based classification methods is to sequentially extract global components form multi-block data and then to use the scores to perform a linear discriminant

	MBA-GUI	utilising	in-house	codes	freely	downloadable	at
	https://www.chem.uniroma1.it/romechemometrics/research/algorithms/.
	Novel use of multi-block methods can also be understood as boosting of different pre-
	processing techniques. Such a methodology called SPORT was recently developed by [41],
					the glob e gl perform a form
					mDim methods Dim meth DA are include are inclu
				that were extra at were ex the decomposit decompo
				ComDim-PLSR omDim-PL ses PLS2 to s PLS2 to
			way by the PLS way by the P trix by ICA, wh rix by ICA,
			m. To use new m. To use s can be direct can be d
		to calculate th alculate ming the new g the new
		d MLR model MLR mod en inputting into inputting
		formed med	usin u
		m.uniroma1. niroma J

analysis (LDA). In the GUI, two variants of ComDim-LDA are included. ComDim-PLS2-LDA simply uses a PLS2 inside ComDim to orient the decomposition of the W matrix. On the other hand, ComDim-PLS2-ICA-LDA uses PLS2 to orient the extraction of the independent components from the W matrix by ICA, which has replaced PCA nested inside ComDim. In both cases, the scores can be directly used for LDA. Test set prediction is performed by first transforming the new data to the same space using the ComDim calibration model and then inputting into the trained LDA model. SO-CovSel is a multi-block variable selection technique recently developed by [16]. The SO-CovSel is an extension of the CovSel technique to the multi-block scenario [40]. CovSel extracts the 'k' variables from a matrix X that are most correlated to Y and independent to each other. SO-CovSel performs CovSel in a sequential orthogonalized way. It works as SO-PLS, replacing the PLS-scores by the COVSEL selected variables. CovSel in such an approach performs variable selection so that extraction of the variables from the consecutive block improves the model. SO-CovSel is designed for both regression as well as classification cases. SO-CovSel for regression and SO-CovSel-LDA are integrated into the J o u r n a l P r e -p r o o f

Table 1 :

 1 The reference dry matter was measured by the hot air oven method by noting the fresh and the dry weight of the samples. More details on the olive fruits data can be found in[44]. The third dataset to demonstrate the classification task consists of the NIR spectra of mayonnaise samples made from oils of 6 different origins. The classification task can be understood as a 6-class problem. The mayonnaise dataset was obtained from the official website of ChemHouse (www. chemproject.org). To make the data fit for multi-block analysis, each dataset was split into two blocks in the spectral domain. A further description of the datasets is provided in Table1. Description of the datasets used in the demonstration.

	pre-treatmen -treatm
	r boosting. Boo oosting. B
	oven to be of h n to be
	m SPORT, the SPORT, the
	can be applied can be app
	SO-PLS-DA is u -PLS-DA
	lti-block anal block an
	s F

The first dataset for data visualisation consists of the FTIR spectra of olive oils of 4 different origins (4 classes) measured in transmission mode

[42, 43]

. The second dataset to demonstrate the regression and variable selection task relates J o u r n a l P r e -p r o o f

  ConclusionA MATLAB based GUI for multi-block data analysis (MBA-GUI) is presented. The toolbox can perform a range of common pre-processing methods on blocks of multivariate data.Multi-block data analysis for regression, classification, visualisation, variable selection and SPORT are proposed. The performance of the MBA-GUI for each of the data analysis tasks was demonstrated with several data sets. The results showed that the MBA-GUI performed well, and all the options are fully functional. The main advantage of the toolbox is that it can be easily understood and used by non-experts. The first version of the GUI can be downloaded at (https://github.com/puneetmishra2/Multi-block.git). Other features will be added to the GUI with the development of new methods. All the data analysis presented in this work can be replicated with the supplied data. The GUI supports data format of .csv, .xlsx and .mat. The users are welcome to notify the authors if they find any bug or problem related to the use of the toolbox, so that the toolbox can be continuously improved along time. The app can be directly installed in MATLAB or can be used as stand by installing the //nl.mathworks.com/products/compiler/matlab-runtime.html). The password to the start the toolbox is "welovedata" without double colon. GUI Toolbox v. 1.0 was successfully executed on three different versions of MATLAB (R2016b, R2017b and R2019a) without any issue to report. The Graphical User Interface (GUI) for data loading and visualisation allows a maximum number of 4 distinct numerical arrays to be imported, represented and pretreated by means of a significant amount of computational tools for smoothing, derivation and other types of corrections (e.g., Standard Normal Variate -SNV -Multiplicative Scatter Correction -MSC -etc.). It is intuitive and very well-conceived.The imported datasets can afterwards processed by a remarkable collection of multi-block latent variable-based dimensionality reduction methodologies permitting tasks of various nature to be carried out:1. Variable selection (Sequential Orthogonalised Covariance Selection and Sequential

	MATLAB The MBA-Orthogonalised Covariance Selection-Linear Discriminant Analysis); 2018b run time compiler of the f th t of rsion o lti-block.git). O block.git) p r o methods. All th ods. A ied data. The G data. Th e notify the a o notify th x, so that the t so that tly installed in installed 018b b J o u r n a l P r e -p r o o f works.com/prod ks.com/p J o veda 2. data exploration/visualisation (Common Dimension -ComDim -Principal tool at Component Analysis, ComDim Canonical Correlation Analysis, ComDim Independent Component Analysis and ComDim Partial Least Squares regression); (https:Validation 3. multivariate regression (Sequential Orthogonalised Partial Least Squares regression
	Dr. Raffaele Vitale and ComDim Regression); and
	1) U. Lille, CNRS, LASIRE, Laboratoire de Spectroscopie pour les Interactions, la Réactivité 4. multivariate classification (Sequential Orthogonalised Covariance Selection-Linear
	et	l'Environnement, Discriminant Analysis, ComDim Independent Component Analysis-Linear Cité Scientifique, F-59000, Lille, France
	2) Molecular Imaging and Photonics Unit, Department of Chemistry, KU Leuven, Discriminant Analysis, ComDim Principal Component Analysis-Linear Discriminant
	Celestijnenlaan 200F, B-3001, Leuven, Belgium
	Email: rvitale86@gmail.com

Analysis and Sequential Orthogonalised Partial Least Squares Discriminant Analysis).

Acknowledgments

PM and AN acknowledge the funding received from the European Union's Horizon 2020 research and innovation program named MODLIFE (Advancing Modelling for Process-Product Innovation, Optimization, Monitoring and Control in Life Science Industries) under the Marie Sklodowska-Curie grant agreement number 675251.


Le Flécher

Simultaneous Control of Two Robotics Arms Sharing Workspace via Visual Predictive Control

Keywords: Agricultural robotics, multi-arms, Non-linear Model Predictive Control, Image Based Visual Servoing

This work addresses the problem of simultaneously controlling two robotic arms to automatically pick up fruits in an orchard. When considering such a scenario, the design of the controller has to face several challenges: (i) the arms share the same workspace, (ii) the presence of obstacles such as branches and other fruits, and (iii) the system evolves in a dynamic environment, i.e., the positions of the fruits and branches change over time. In this paper, it is proposed to control the two robotic arms relying on a Visual Predictive Control (VPC) scheme. It allows to control the system in a reactive fashion while taking into account a large number of constraints. The global and local models of the VPC scheme as well as constraints dealing with collisions or joint boundaries are presented. Finally, the efficiency of the proposed approach is highlighted with numerous simulation results using the PR2 arms model.

Introduction

To meet the increasing food demands from a growing world population, agriculture will need to double its production by 2050 [START_REF] Foley | Solutions for a cultivated planet[END_REF][START_REF] Grift | A review of automation and robotics for the bio-industry[END_REF]. This cannot be achieved by simply doubling the inputs (land, water, seeds, labor, etc.) because of constrained resources and environmental concerns. Among the solutions presented in [START_REF] Foley | Solutions for a cultivated planet[END_REF][START_REF] Grift | A review of automation and robotics for the bio-industry[END_REF] to increase the global agricultural total factor productivity, robotics is one with the potential to effect agriculture in a broad, systemic way, and contribute significantly to meeting our future needs. Up to today, the three most significant impacts robotics and automation have had on agriculture are: (i) Precision agriculture, i.e., the use of sensors to precisely control when and where to apply fertilizers or water; (ii) Auto-guidance on field crop machinery, which today can drive down a field with an accuracy unattainable by human drivers; (iii) Machines that harvest fruits and vegetables for processing (e.g., tomato paste and orange juice) [START_REF] Bergerman | Robotics in agriculture and forestry[END_REF]. Research is now focused on the next wave of sensing, mobility, and manipulation technologies that promise to increase agriculture output and productivity.

Over the last years, one has focused on robotic systems performing weeding, spraying, trimming, or harvesting in orchards. In order to achieve any of these tasks, the following capacities are expected for the robotic system: (i) autonomously navigate in an orchard; (ii) detect and localize a set of natural landmarks of interest; (iii) approach the landmarks with some robotic manipulators; (iv) manipulate/interact with the landmarks of interest without damaging them. After addressing the autonomous navigation problem in [START_REF] Durand-Petiteville | Design of a sensor-based controller performing u-turn to navigate in orchards[END_REF][START_REF] Flecher | Implementation on a harvesting robot of a sensor-based controller performing a uturn[END_REF], the focus is now on the landmark approach problem when considering numerous robotic manipulators embedded on a mobile platform. This problem has already been addressed to design harvesting systems for cucumbers [START_REF] Van Henten | Collision-free Motion Planning for a Cucumber Picking Robot[END_REF] or tomatoes [START_REF] Zhao | Dual-arm Robot Design and Testing for Harvesting Tomato in Greenhouse[END_REF]. However, despite the promising results, the obtained harvesting throughput does not match with the farmers needs and must be increased [START_REF] Zhao | A review of key techniques of visionbased control for harvesting robot[END_REF]. To achieve this aim, one of the most adopted approach consists of embedding several robotic arms on the mobile platform, as it is done in [START_REF] Zhao | A review of key techniques of visionbased control for harvesting robot[END_REF] for example. However, the throughput increase is not proportional to the number of manipulators, it also has to do with their distribution on the machine and their workspace. Indeed, it is pointed in [START_REF] Vougioukas | A Study of Fruit Reachability in Orchard Trees by Linear-Only Motion[END_REF] that fruits are not uniformly distributed on commercial orchards trees. To maximize the fruits accessibility, growers use specific trees with high fruits density, thus facilitating the manual harvesting. As a consequence, when designing a harvesting embedding multiple arms, it is also important to consider the workspace of each arm and the possible overlapping. Thus, even if the arms do not collaborate to manipulate one single object, they share a common workspace. It is then possible to use several arms to pick fruits in the same high density area, which increases the overall throughput.

In this work, it is proposed to design a control scheme for two robotic arms aiming at positioning their end effector to achieve an agricultural task in an orchard. When considering performing a task in orchard, it seems relevant to focus on reactive controllers. Indeed, a tree is a highly dynamic environment because of the flexibility of the branches and the weight of the fruits: each time a fruit is harvested, the local distribution of fruits and branches is modified. Among the available control scheme, it is proposed to use an image based visual servoing (IBVS) scheme [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF]. It requires the use of a camera, and nowadays it is simple to obtain a camera small enough to be install close to the end-effector without creating any trouble. Moreover, it seems to be the appropriate sensor to track the fruit of interest. IBVS consists of computing the error between the current coordinates of points of interest and a set of desired values. This error is then minimized via a controller based on the image Jacobian matrix, allowing to control the displacements of the end effector. Computing the control law directly in the image space makes this approach robust to modeling errors and reactive objective modification, however it does not provide any guarantee regarding the trajectory of the camera in the workspace. This represents a major issue when considering arms sharing a common workspace.

Numerous works have been developed in order to modify the camera trajectory in the workspace or to handle issues such as the loss of the visual features, collisions with obstacles,... For example, the work presented in [START_REF] Chaumette | Asservissement visuel. In: La commande des robots manipulateurs[END_REF] proposes to use advanced visual features to obtain more appropriate camera trajectories. However, they do not guarantee the non-collision with obstacles. Another approach consists of coupling IBVS with other elementary controllers, and to switch from one controller to the other based on the current situation (obstacles, occultations, joint boundaries, ...) [START_REF] Cadenat | Comparison of two sequencing techniques to perform a vision-based navigation task in a cluttered environment[END_REF] [9] [START_REF] Durand-Petiteville | A complete sensor-based system to navigate through a cluttered environment[END_REF]. This approach requires to design numerous ad-hoc controllers, and defining the switching rules can happen to be tenuous. Other works use the robot redundancy to modify the trajectory in order to satisfy the given constraints [START_REF] Marchand | Using the task function approach to avoid robot joint limits and kinematic singularities in visual servoing[END_REF] [START_REF] Chaumette | A redundancy-based iterative approach for avoiding joint limits: application to visual servoing[END_REF]. However, this method can only be used if the main task does not use all the degrees of freedom [START_REF] Mansard | Task sequencing for high-level sensorbased control[END_REF], which might be the case when using systems with large number of degrees of freedom, such as humanoid robots. Finally, it has been proposed to couple IBVS with path planning [START_REF] Mezouar | Path planning for robust imagebased control[END_REF], [START_REF] Kazemi | Path-planning for visual servoing: A review and issues[END_REF], in order to generate a path free of occultations and collisions. However, to do so, it is mandatory to have a model of the environment prior to the task execution, which might be challenging in a dynamic environment.

In this work, it is proposed to rely on a well adopted solution allowing to take into account constraints on position, velocity, or workspace: the model predictive control scheme [START_REF] Grüne | Nonlinear model predictive control[END_REF]. This approach relies on the minimization of a cost function made of the predicted states of the system and defining the task to achieve. Moreover, when minimizing the cost function with a numerical solver, this latter can handle numerous constraints [START_REF] Grüne | Nonlinear model predictive control[END_REF]. Thus, to obtain a reactive controller dealing with constraints, it is proposed to design an IBVS-based NMPC (Nonlinear Model Predictive Control) controller, also named Visual Predictive Control (VPC). For this particular case, the cost function is defined as the error between the image coordinates of the visual features representing the landmark of interest at the current pose and the desired one. To develop a VPC scheme, it is then mandatory to (i) obtain the prediction models of the system state, (ii) design the cost function defining the task to achieve, (iii) design the mandatory constraints to achieve the task, and (iv) select the numerical solver allowing to minimize the cost function under the given constraints. Thus, this paper addresses the aforementioned points for the specific case of two simultaneously controlled robotic arms. Moreover, it is an extended version of the work presented in [START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF] which has been complemented with an analysis of the related works, a discussion regarding the terminal constraint and a more comprehensive set of simulation results.

The paper is organized as follows. First, section II is dedicated to an overview of the existing VPC scheme. Next, in section III, one presents the different used predictive models, before addressing in section IV the design of the cost function and of the different constraints. In section V, a set of simulations is presented to support the proposed approach and highlight its efficiency. Finally, in section VI, one summarizes and analyses the obtained results before developing the expectations for the future.

Related Works

As it has been explained, NMPC schemes require to use a model to predict the future states of the system. In the case of VPC, it consists of predicting the coordinates of the visual features knowing a sequence of control inputs. Several works have addressed the design of VPC schemes, and according to [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF], they can be classified into two categories: the local and global models. The first ones relies on the interaction matrix (or image Jacobian) [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF], whereas the second ones use the camera pose. For both cases, the models are built considering a flying camera, i.e., the system embedding the camera is not taken into account in the model. Finally, none of the models proposes an analytic form for the prediction of the visual features. Thus, the predicted values are obtained using a linear approximation, or with advanced numerical integration schemes, such as Runge-Kutta.

Most of the existing works belong to the first category, i.e., relying on local models. For example, in [START_REF] Assa | Robust model predictive control for visual servoing[END_REF] the authors use a first order approximation to control a camera embedded on a robotic arm. Similar approaches are presented in [START_REF] Heshmati-Alamdari | Robustness analysis of model predictive control for constrained image-based visual servoing[END_REF] for a flying camera and in [START_REF] Ke | Visual servoing of constrained mobile robots based on model predictive control[END_REF] for a mobile robot. In [START_REF] Wang | Quasi-min-max model predictive control for image-based visual servoing[END_REF] and [START_REF] Hajiloo | Robust On-Line Model Predictive Control for a Constrained Image Based Visual Servoing[END_REF], the visual servoing scheme is represented by a linear parameter varying (LPV) system. The objectives of these works is to demonstrate the feasibility of VPC subject to constraints on the camera field of view or command vector. The results show the relevance of the approach, but are obtained for small displacements between the initial and desired pose, minimizing the impact of the predictions. On the opposite, the work presented in [START_REF] Mcfadyen | Aircraft collision avoidance using spherical visual predictive control and single point features[END_REF] proposes to control an UAV while avoiding other aircraft. In this case, the distance between the vehicles is large enough to use an approximate value of the depth in place of the real one.

In this work, it is proposed to first perform an evaluation of the accuracy of the prediction models using the real or approximated depth values. This evaluation is done in the context of a camera embedded on a robotic arm. Following this evaluation, we then present the design of our VPC scheme to simultaneously control two robotic arms.

Modeling

In this section, one first focuses on the modeling of the system, i.e., the robotic arms, the cameras and the targets [START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]. Then, two models to predict the visual feature coordinates of a point are presented [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF] [START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]. One ends with a discussion regarding the viability of these models for the visual predictive control problem.

System modeling

The robotic system considered in this work is composed of two identical arms, each of them containing m q revolute joints. A camera is attached to the end effector of each arm. To model the system, let us first define a global frame

F b = (O b , x b , y b , z b ) attached to the system base. Moreover, two frames F ci = (O ci , x ci , y ci , z ci ), with i ∈ [1, 2]
, are used to represent the pose of the camera attached to the i th arm.

Let define q ij and qij , with j ∈ [1, ..., m q ], respectively the angular position and velocity of the j th joint of the i th arm. Thus, for the i th arm, one obtains the configuration vector Q i = [q i1 , q i2 , ..., q imq ] T and the control input vector Qi = [ qi1 , qi2 , ..., qimq ] T . Finally, for the whole system one obtains

Q = [Q T 1 , Q T 2 ] T and Q = [ QT 1 , QT 2 ]
T . The cameras mounted on the end effectors are modeled using the pinhole model. Thus, the projection matrices H Ii/ci mapping the 3D coordinates (x, y, z) of a point in the camera frame F ci to its 2D projection (X, Y ) on the image plan F Ii is defined as follows:

    X Y z 1     =     f /z 0 0 0 0 f /z 0 0 0 0 1 0 0 0 0 1         x y z 1     (1)
In this work, the arms are controlled to make each camera reach a pose defined by the coordinates of point visual features in the image space. To do so, one uses two landmarks made of four points. When considering the i th camera/target couple, the coordinates of the projection of each point on the image is denoted [START_REF] Assa | Robust model predictive control for visual servoing[END_REF][START_REF] Bergerman | Robotics in agriculture and forestry[END_REF][START_REF] Cadenat | Comparison of two sequencing techniques to perform a vision-based navigation task in a cluttered environment[END_REF]. Thus, the visual features vector for each camera is defined as S i = [S i1 , S i2 , S i3 , S i4 ] T , whereas the system one is given by S = [S T 1 , ..., S T na ] T . In the same way, a vector for the visual features desired coordinates of each camera is defined as

S il = [X il , Y il ] T , with i ∈ [1, 2] and l ∈ [1,
S * i = [S * i1 , S * i2 , S * i3 , S * i4 ] T .

Prediction of the visual features

A visual predictive control scheme requires to be able to predict the values of the visual features for a given sequence of control inputs. In [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF], the authors present two ways to perform the prediction step for a flying camera: the local and global approaches. This section is devoted to the presentation of both methods for the case of a 6 degrees of freedom camera embedded on a robotic arm.

The global model Global models rely on the pose of the camera. It can be decomposed as follows: from the image coordinates, the location of the visual features in the current camera frame is calculated (2D to 3D); the visual features coordinates are then computed in the predicted camera frame (3D to 3D); and finally the visual features are projected on the predicted image plane (3D to 2D). To do so, the camera pose is generally predicted by integrating the camera kinematic screw [START_REF] Allibert | Predictive Control for Constrained Image-Based Visual Servoing[END_REF]. In this work, one integrates the robotic arm, which allows to predict the camera pose based one the robot command vector, minimizing the approximation. To do so, let us consider a discrete time system with a sampling period T s , where t k+1 = t k + T s , and the notation .(k) = .(t k ). The configuration

Q i (k + 1) of the i th arm, with i ∈ [1, 2]
, at the next iteration obtained after applying the control inputs Qi (k) is given by:

Q i (k + 1) = Q i (k) + Qi (k)T s (2)
Knowing the arm configuration at instants t k and t k+1 , it is then possible to compute the homogeneous transformation matrices between the camera frame F ci and the base one

F b . ōF b (.) = H b/ci (.)ō Fc i (.) (3)
where ō are the homogeneous coordinates of a point feature expressed in the base or camera frame. Finally, coupling (3) with ( 1), one obtains a prediction of the coordinates of a point feature in the image space.

S il (k + 1) = H Ii/ci (k + 1)H -1 b/ci (k + 1)H b/ci (k)H -1 Ii/ci (k)S il (k) (4)
The local model Local models rely on the differential equation mapping the derivative of the visual feature coordinates Ṡil to the robot control input vector Qi . To obtain this equation, one uses the robot's Jacobian J i , mapping the camera kinematic screw T ci and Qi [START_REF] Chaumette | A redundancy-based iterative approach for avoiding joint limits: application to visual servoing[END_REF], and the interaction matrix L il , mapping Ṡil to T ci [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF].

T ci = J i Qi (5) Ṡil = L il T ci (6)
where

L il = -1 z l 0 X l z l X l Y l -(1 + X 2 l ) Y l 0 -1 z l Y l z l 1 + Y 2 l -X l Y l X l (7)
By combining equations ( 5) and ( 6), one obtains the following local model for each arm/camera couple:

Ṡil = L il J i Qi (8)
Usability of the models When considering using one of the models on an experimental system, several aspects have to be taken into account. First, both models require a measure/estimate of the visual feature depth z (see equations ( 1) and ( 7)). However, when using a monocular camera, this value can only be retrieved by estimation. Another widely adopted option consists of using a constant value, such as the depth value at the initial or desired pose [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF].

The second issue is related to the integration of the local model to obtain the prediction equation. A first option consists of using an advanced numeric scheme to accurately compute the integration, whereas a second one relies on a linear approximation of the model, such as:

S il (k + 1) = S il (k) + L il J i Qi (k)T s (9)
Thus, depending on the user choices regarding these issues, the prediction feature coordinates will be more or less accurate. In section 5, different configurations will be evaluated for our specific application in order to select the most appropriate and realistic solution.

Controller Design

In this section the control of the arms motion to reach their own targets and compounding necessary constraints are presented. To achieve these goals, NMPC control has been implemented to find the optimal inputs vector which minimizes a cost function subject to constraints. To do so, we first present the cost function to be minimized in order to reach the desired targets. Then, the expression of each constraint will be defined. A similar approach was presented in [START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF] for a multi-arm robotic system.

Visual Predictive Control

The proposed VPC scheme couples NMPC (Nonlinear Model Predictive Control) with IBVS (Image-Based Visual Servoing). Thus, similarly to NMPC, it consists of computing an optimal control sequence Q(.) that minimizes a cost function J N P over a prediction horizon of N P steps while taking into account a set of userdefined constraints C( Q(.)). Moreover, similarly to IBVS, the task to achieve is defined as an error in the image space. Thus, the cost function to minimize is the sum of the quadratic error between the visual feature coordinates Ŝ(.) predicted over the horizon N p and the desired ones S * . It then possible to write the VPC problem as:

Q(.) = min Q(.) J N P (S(k), Q(.)) (10)
with

J N P (S(k), Q(.)) = k+Np-1 p=k+1 [ Ŝ(p) -S * ] T [ Ŝ(p) -S * ] (11)
subject to

Ŝ(.) = (4) or (8) (12a) Ŝ(k) = S(k) (12b) C( Q(.)) ≤ 0 (12c)
As mentioned in the previous section, the visual feature coordinates can be predicted either using the global or local models (equation (12a)). To compute the prediction, one uses the last measured visual features (equation (12b)). Finally, equation (12c) defines the set of constraints to be taken into account when minimizing the cost function [START_REF] Foley | Solutions for a cultivated planet[END_REF]. Thus, solving equation ( 10) leads to the optimal sequence of control inputs Q(.). In this work, and as it is usually done, only the first element Q(1) of the sequence is applied to the system. At the next iteration, the minimization problem is restarted, and a new sequence of optimal control inputs is computed. This loop is repeated until the task is achieved, e.i., the cost function ( 11) is vanished.

Constraints

As seen in the previous section, it is possible to add a set of constraints in the minimization problem. In this section, one presents the constraints related to the geometry of the robot, the velocity of the joints, and the field of view of the camera. One finishes with defining a set of constraints that guarantee the non-collision of the arms despite sharing a common workspace.

System constraints

The joint velocity lower and upper bounds are respectively denoted Qmin and Qmax . They are two N Q long vectors, with N Q = 2m q N p . They allow to define the lower and upper limits of each joint of each arm over the prediction horizon. Consequently, the constraint on Q(.) is described in the domain κ = [ Qmin , Qmax ], and the velocity constraint is then written as:

Q -Qmax Qmin -Q ≤ 0 ( 13
)
Mechanical limits follow the same principle. The lower and upper bounds of the angular values are denoted Q min and Q max . They are two N Q long vectors, with N Q = 2m q N p . These constraints are then given by:

Q -Q max Q min -Q ≤ 0 (14)
Finally, the visual constraints are defined by the image dimension. One denotes the bounds by S min and S max . They are two N S long vectors, with N S = 8n a N p . Thus, all the predicted visual features coordinates need to be kept into this space, and the constraints can then be written as:

S -S max S min -S ≤ 0 ( 15
)
Overlapping workspace constraint Due to the overlapping workspace, it is mandatory to guarantee the non-collision between the arms. To do so, one computes the shortest distance d i,j|i ,j between the j th joint of the i th arm and the j th joint of the i th arm. The closest distance is evaluated as the norm of the segment perpendicular to both assessed links as it is illustrated in figure 1.

Finally, in order to avoid collisions, one defines the following constraint:

d min -d ≤ 0 ( 16
)
where d = [d 1,1|2,1 , ..., d 1,mq|2,mq ] and d min is a pre-defined safe distance. Terminal constraint When using a NMPC scheme, it is necessary to guarantee that the set of control inputs makes the camera converge towards the desired pose. Due to the fact that a NMPC scheme minimizes the distance between a predicted trajectory and a desired one, there is no guarantee that the ultimate predicted state has converged towards the desired one. Moreover, the prediction horizon might be too short or the constraints on the control inputs might be too restrictive to reach the goal. [START_REF] Grüne | Nonlinear model predictive control[END_REF]. Thus, to check the feasibility, one adds a final constraint in the optimization problem. It is defined as the error between the prediction of the visual feature coordinates Ŝ(k + N p -1) obtained at the end of the prediction horizon, and the desired ones S * . If the solver cannot compute a sequence of control inputs Q that respects this constraint, then the problem is not feasible, and there is no guarantee regarding the system convergence.

| Ŝ(k + N p -1) -S * | -δ tc ≤ 0 ( 17
)
where δ tc is a user defined threshold. Two remarks can be made concerning the use of the terminal constraint. First, none of the works studied during the literature review proposes to include a terminal constraint, despite its ease of implementation and utility. Moreover, authors usually weight the last predicted value based on the distance to the desired one. This method helps the solver to converge towards an optimal solution but it does not guarantee the convergence as the terminal constraint does. Second, the terminal constraint guarantees the convergence of the predicted states towards the desired ones. In the case the predictions are strongly erroneous, then the real system cannot converge towards the desired sate. Thus, using a terminal constraint does not provide an absolute guarantee of the task realization.

All the constraints are now defined and can be integrated into the non-linear function C( Q(.)):

C( Q) ≤ 0 ⇒             Q -Qmax Qmin - Q Q -Q max Q min -Q S -S max S min -S Dmin -D | Ŝ(k + N p -1) -S * | -δ tc             ≤ 0 (18)
Thus, the optimal input must belong to the domain defined by C( Q(.)).

Simulation

In this section, one presents the results obtained by simulating the control of two arms with MATLAB software. A 7 DOF model has been selected, based on the arms embedded on the PR2 robot. In these simulations, the task consists of positioning each camera with respect to a given target. Each final pose is defined in the image space by S * i = [x 1 * , y 1 * , x 2 * , y 2 * , x 3 * , y 3 * , x 4 * , y 4 * ] T , which are the coordinates of four points representing the corresponding target. Moreover, one setups the prediction horizon to N p = 5. Regarding the constraints, the angular velocity limits for the first prediction step are given by κ = [-0.2 0.2] for each joint, whereas it is extended to κ = [- 5 5] for the rest of the prediction step. This allows to decrease the computation time while guaranteeing a feasible control input for the first step. The size of the image is S min = [-4 -4] and S max = [ 4 4]. The angular constraints on each joint are setup as described in the PR2 documentation [START_REF] Garage | Pr2 user manual[END_REF]. Moreover, the minimum collision avoidance distance is fixed to D min = 0.1. The initial coordinates of the visual features are computed based on the initial configuration of the robot given by q init = [0 π 2 0 0 0 0 0 0 π 2 0 0 0 0 0] T for every simulation. Finally, one fixes T s = 50ms.

Comparison between the different prediction models

The first set of simulation assesses the accuracy of different predictive models presented in section 3.2. To do so, the two robotics arms are controlled via VPC from q init to q * = [-π 8 π 2 0 0 0 0 0 0 π 2 0 0 0 0 0] T , using the global model and the real depth values. When simulating a VPC scheme, this configuration can be considered as the ground truth, as there is no approximation in the calculations. Next, at each iteration, the optimal sequence of control inputs calculated by the VPC scheme, is used to compare the different prediction models. To do so, the predicted values at the end of the prediction horizon using a particular model are compared with the reference ones.

In figure 2, one can see the evolution of the quadratic error between the predicted coordinates of one visual feature and the ground truth. Six models are presented by combining the global model (solid lines), local Euler approximation (dashed lines) and local Runge Kutta approximation (dotted lines) with the depth values real z (black) and desired z * (blue). First, it can be seen that the error decreases over the servoing. Indeed, when the camera is close to the desired pose, the sequence of optimal control inputs is quite small, as there is no much distance left to perform. Having small control inputs then reduces the error in the predicted values. Moreover, at the end of the simulation, z(t) and z * are very similar, which explains the convergence of the blue and black curves for each model. Next, it can be noticed that the global models and the local Runge Kutta are significantly more accurate than the local Euler one. Indeed, using this latter prediction scheme introduces a large error, and might not allow the system to converge toward the desired configuration. In figure 3, the errors are summed over the simulation to offer a better view of model performances. When considering the use of the real depth, the global model is the most accurate. However, this result is obtained in a simulation, and its use in a real context would require an estimator of the depth. Thus, one now considers one more realistic scenario where the depth is approximated with the constant value z * . In this case, the local Runge Kutta is the most accurate model. Thus, for the rest of the simulations the local Runge Kutta model with z * is used to evaluate the performances. dealing with the previously presented constraints. The VPC scheme relies on the local model coupled with a Runge Kutta algorithm and the desired depth z * . Using this configuration, the VPC scheme achieves to control the robotic system to converge towards the desired poses, as it can be seen in the different results presented hereafter. First, the evolution of visual features is displayed in figure 4 and highlights the achievement of the task for both cameras. Indeed, the proposed scheme was able to control both arms to make the current visual features converge towards the desired ones from their initial locations (respectively represented by the red and green crosses in figures 4(a) and 4(b)). To achieve the task, the arms had to reach the configurations presented in figure 5.

Regarding the constraints, the image limits were never reached and there was no risk of collision during the servoing. Moreover, it can be seen in figures 6(b) and 6(a) that the joint angular values remain within the defined domain. Similarly, in figures 6(d) and 6(c), one can see that the the joint velocities were bounded by the constraints. Indeed, the command inputs stay within the given range κ = [-0.2 0.2] all along the simulation.

Auto-collision scenario

This new simulation is setup in order to create a collision between the two robotic arms. Thus, the system has to reach the desired configuration

q * = [ π 16 3π 8 0 π 4 0 π 8 0 -π 16 π 2
π 4 π 8 0 0 0] T . Moreover, as previously done, one uses the Runge Kutta local model coupled with the desired depth values. As it can be seen in figure 7, one more time the desired visual features are reached from the initial poses. However, as shown in figure 8, the arms had to cross each other to reach the desired locations. Thanks to the shared environment constraint included in the control problem, the collision was avoided. Indeed, figure 9 shows that the distance between the different links was never smaller than the given value D min = 0.1m. In addition to the auto-collision avoidance, the other constraints on the joints angular values and velocities were respected (see figures 10 and 11).

External collision scenario

In this scenario, the goal is to achieve the task while dealing with an obstacle on the right arm trajectory. First, one defines the desired pose such as q * = [-π scene at the coordinates ŌObs| F b = (0.79, -0.27, -0.08) (see figure 12). Finally, the minimal distance is setup to 0.1m.

As it can be seen in figure 13(a), the positioning task is achieved, and the visual features have reached the desired coordinates. Moreover, the corresponding trajectory is plotted in figure 12. It should be noticed that the obstacle on the right arm trajectory has been avoided thanks to the distance constraint. Indeed, as it can be seen in figure 13(b), the distance between any segment of the arm and the obstacle is always higher than the defined threshold. Thus, it is possible to safely achieve the positioning task despite the presence of unknown obstacle.

* = [ π 20 π 2 -π 2 π 6 π π 4 0 -π 20 π 2 π 2 π 6 0 0 0].
The obstacle is positioned at the coordinates ŌObs F b = (0.79, -0.27, -0.08) (see figure 14) and the minimal distance for both arms and obstacle is set up to 0.1m.

One more time, the task is successfully achieved as it can be seen in figure 15. Indeed, the visual features manage to converge towards the desired values despite the presence of an obstacle and the possible auto-collision. The corresponding trajectory is plotted in figure 14. This trajectory is obtained by taking into

Conclusion

In this paper, the problem of fruit picking using a dual-arm robot in a commercial orchard has been considered. A strategy based on VPC has been developed to make its end-effectors reach a desired pose close to the corresponding fruits, while respecting constraints due to the environment dynamics, the visual manipulation and the shared workspace. This strategy thus allows to benefit from the advantages of a reactive controller to perform a highly complex task in a strongly evolutive environment where various constraints must be taken into account. To validate the approach, a simulation campaign has been conducted. A first step was to evaluate the different prediction models to select the best one for the servoing. Then, the VPC strategy itself was tested, considering different scenarii and using a PR2 model as the considered dual-arm robot. The obtained results have shown its interest and its efficiency to safely perform various positioning tasks in a shared workspace. For future works, two main leads can be followed. First, these works only consider a fixed dual-arm system. To go further, it is necessary to mount it on a mobile base and to coordinate the motions of the complete robotic setup. To do so, the proposed control strategy will have to be coupled to our previous works about autonomous navigation in orchards. A second interesting research axis concerns the realization of experimental tests on a real robotic system. To achieve this goal, it will be required to speed up the processing time. The prediction models will have to be selected depending not only on their accuracy but also on their computational time. In this context, designing new models specifically adapted to Graphical Processing Units could offer a nice opportunity. Such aspects are now at the core of our current work in order to perform a complete fruits picking task in a real environment.

Fig. 1 .

 1 Fig.1. Computation of the shortest distance d i,j|i ,j between links Lij and L i j . Source[START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]

Fig. 2 .Fig. 3 .

 23 Fig. 2. Quadratic errors of prediction over the simulation -Solid: Global model -Dashed: Local Euler -Dotted: Local Runge Kutta -Black: z -Blue: z * . Source[START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]

Fig. 4 .Fig. 5 .

 45 Fig. 4. Evolution of the visual features using local Runge-Kutta model and z -Blue dotted: Trajectories -Red circles: Final locations -Green circles: Initial locations.Source[START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 6789 Fig.6. Joint configurations and velocities. Source[START_REF] Flécher | Visual predictive control of robotic arms with overlapping workspace[END_REF]

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Configuration of the right arm over the simulation -Blue solid: Angular joint value q1j -Red dashed: Angle boundaries. Source [10]

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Final configurations of the robotic arms

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Final configurations of the robotic arms

Fig. 16 .

 16 Fig. 16. Collision distance between the right arm and the obstacle, and between the arms