Energy Efficiency Optimization in LoRa Networks -A Deep Learning Approach

Lam-Thanh Tu, Abbas Bradai, Olfa Ben Ahmed, Sahil Garg and Yannis Pousset

Abstract-The optimal transmit power that maximizes energy efficiency (EE) in LoRa networks is investigated by using deep learning (DL) approach. Particularly, the proposed artificial neural networks (ANNs) is trained two times; in the first phase, the ANNs is trained by the model-based data which are generated from the simplified system model while in the second phase, the pre-trained ANNs is re-trained by the practical data. Numerical results show that the proposed approach outperforms the conventional one which directly trains with the practical data. Moreover, the performance of the proposed ANNs under both partial and full optimum architecture are studied. The results depict that the gap between these architecture is negligible. Finally, our findings also illustrate that instead of fully retrained the ANNs in the second training phase, freezing some layers are also feasible since it does not significantly decrease the performance of the ANNs.

Index Terms-Energy Efficiency, LoRa Networks, Stochastic Geometry, Deep Learning, Poisson Cluster Process.

I. INTRODUCTION

Energy efficiency (EE) is one of the long lasting problems in wireless communications systems. However, in the past, the network operator/planning primarily focused on maximizing the spectral efficiency (SE) as well as enhancing the coverage area, as a result, the maximum transmit power was typically considered as one of the optimal solutions. Nonetheless, such the approach, of course, reduces dramatically the energy efficiency. The issue is even more serious since it is expected that by 2022 there will be 12.3 billion end-devices (EDs) to be connected to the wireless networks [START_REF] Index | Cisco visual networking index: Global mobile data traffic forecast update, 2017-2022 white paper[END_REF] and making the internet-of-things (IoTs) becomes feasible. As a consequence, minimizing the power consumption or maximizing the EE while guaranteeing the SE has been emerged as one of the most important issues in wireless networks. Despite the fact that the operation of the base stations (BSs) has been optimized to maximize the EE of the cellular networks [START_REF] Tu | On the Energy Efficiency of Heterogeneous Cellular Networks with Renewable Energy Sources -A Stochastic Geometry Framework[END_REF], superdense deployment makes its too bulky, thus, letting its less attractable for supporting such a massive low power networks. Fortunately, the low power wide area network (LPWAN) is regarded as one of the most suitable technologies thanks to its properties, i.e., low power consumption, low cost and wide coverage area [START_REF] Goursaud | Dedicated networks for IoT: PHY/MAC state of theart and challenges[END_REF]. Among all the available LPWAN techniques, i.e., SigFox, Weightless, NB-IoTs, etc., LoRa is gained lots of attraction in both academia and industry. By actively fine-tuning its parameters, i.e., the spreading factor (SF), the coding rate (CR) and the bandwidth (BW), LoRa is L.-T. Tu, A. Bradai, O. Ben Ahmed and Y. Pousset are with the Institute XLIM, University of Poitiers, France (email: lam.thanh.tu@univpoitiers.fr; abbas.bradai@univ-poitiers.fr; olfa.ben.ahmed@univ-poitiers.fr; yannis.pousset@univ-poitiers.fr).

S. Garg is with Ecole de technologie superieure, Universite du Quebec, Montreal (email: sahil.garg@ieee.org).

able to serve a wide range of IoTs applications/devices, e.g., e-Health, smart city, smart home, that, in general, have different quality-of service (QoS).

Deep learning (DL), on the other hand, is proven itself as one of the best ones among all machine learning (ML) techniques when a large number of data are available. Nevertheless, unlike other domains where the mathematical modelling is primarily difficult to be performed and the pure data-driven approach like deep learning seems to be the sole solution, wireless communications, on the contrary, have always depended on a robust mathematical modelling for system design, analyzing, optimization and can be considered as model-based approach. Nonetheless, as the dramatic evolution of the wireless networks, i.e., the exponential growth of enddevices with different applications are connected to wireless networks, making its more complex, hence, the mathematical modelling is steadily losing its accuracy as well as mathematical tractability. As a result, deep learning has recently been commenced applying into wireless communications. The application of DL in wireless networks, differently, does not mean that the prior mathematical modelling is ignored. In fact, the most feasible solution of applying DL in wireless communications is to combine the advantages of both datadriven and model-based approach [START_REF] Zappone | Model-Aided Wireless Artificial Intelligence: Embedding Expert Knowledge in Deep Neural Networks for Wireless System Optimization[END_REF]. The main target of this combination is to synergistically exploit the deep and expert knowledge from theoretical models even though it may be inaccuracy and cumbersomeness to facilitate the use of DL in wireless networks [START_REF] Zappone | Wireless Networks Design in the Era of Deep Learning: Model-Based, AI-Based, or Both?[END_REF]. This approach can also be considered as transfer learning [START_REF] Pratt | Discriminability-based transfer between neural networks[END_REF], where the artificial neural networks (ANNs) is first trained by faults data and then re-trained by the empirical data. Of course, such the approach requires a strong mathematical modelling which is as close as possible to the empirical data in order to attain the highest performance. As a result, in this paper, we maximize the energy efficiency in LoRa networks by combining the model-based and data-driven approach.

The performance of the energy efficiency in LoRa networks were studied in [START_REF] Beltramelli | LoRa beyond ALOHA: An Investigation of Alternative Random Access Protocols[END_REF], [START_REF] Su | Energy Efficient Uplink Transmissions in LoRa Networks[END_REF] In [START_REF] Beltramelli | LoRa beyond ALOHA: An Investigation of Alternative Random Access Protocols[END_REF], the energy efficiency was investigated by considering other medium accesses rather than pure ALOHA. By jointly exploiting user scheduling and SF assignment, the maximization of system EE was investigated in [START_REF] Su | Energy Efficient Uplink Transmissions in LoRa Networks[END_REF]. In spite of the use of DL has recently been received tremendous attention in cellular networks, its application in LoRa networks is still in infancy stage. There were a few works which applied ML/DL in LoRa networks, [START_REF] Sandoval | Optimizing and Updating LoRa Communication Parameters: A Machine Learning Approach[END_REF], [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF]. Particularly, in [START_REF] Sandoval | Optimizing and Updating LoRa Communication Parameters: A Machine Learning Approach[END_REF], the process of networks configuration was formulated as a rerinforcement learning (RL) problem. The time difference of arrival positioning method in LoRa networks was improved by applying deep learning in [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF].

In the present work, different to these above-mentioned papers, we maximize the energy efficiency in LoRa networks respect to the transmit power under practical scenarios where the distribution of the end-devices and the small-scale fading follow general distributions, i.e., the Poisson cluster processes (PCP) and the Nakagami-m distribution, combined with the imperfect power consumption at the end-devices. Of course, it is crystal clear that the optimal transmit power that maximizes the EE under such the general system model can not be attained based on the mathematical frameworks, thus, the deep learning approach is utilized instead. Our proposed DL approach (called double training approach), however, requires a relatively small practical data to attain the high accuracy output and is different to the conventional approach (called direct training approach) where a large number of practical data is imperatively needed in order to fully train the ANNs. To realize such the well-trained neural networks with a relative small practical data, the proposed artificial neural networks is trained two times. In the first training phase or pre-trained phase, the proposed ANNs is trained based on the model-based data and in the second training phase, the pre-trained ANNs is re-trained based on the practical data. Here, the model-based data in the first training phase are generated based on the mathematical frameworks under the simplified system model that may be inaccuracy but may provide some expert knowledge and can be helpful for the second training; the practical data in the second training phase, diversely, are uniquely generated by Monte Carlo simulation that is resources-consuming compared to the mathematical frameworks. To be more specific, the main contributions and novelties of this paper are summarized as follows: i) the distribution of the end-devices and the smallscale fading follow Poisson cluster processes and Nakagamim distribution, respectively; the power consumption at EDs is impaired by an additive noise that follows either Gaussian or uniform distribution; ii) the closed-formed expression of the energy efficiency under the simplified system model is provided; iii) the first training set is generated by deploying the mathematical frameworks under the simplified system model. Particularly, the optimal transmit power in this phase is attained by numerical solving a non-linear equation; iv) the proposed double training approach outperforms the direct training approach under all considered metrics, i.e., the mean square error (MSE) and R squared, R 2 , or the coefficient of determination; v) our findings show that increasing training data monotonically ameliorates the performance of the neural networks. Nonetheless, enhancing the ANNs performance by raising either the number of epochs or the number of neurons does not always benefit; vi) the performance of the ANNs under various optimal network architecture, i.e., partial and full optimum architecture, are investigated. The findings depict that the gap between these optimum architecture is negligible, thus, the partial optimum architecture is preferable provided that the resources consumption are taken into account; and vii) our findings also illustrate that it is feasible to freezing some layers during the second training phase which does not significantly affect the performance of the ANNs.

The rest of this paper is organized as follows. In Section II, the considered system model is presented. In Section III, we formulate the optimization problem and the design of the proposed neural networks is provided in Section IV. In Section V, the performance of the proposed ANNs is evaluated and discussed under various scenarios. Finally, Section VI concludes the paper. Notations: Pr (.) and E {.} are the probability and the expectation operators; max {.}, min {.} and exp (.) being the maximum, minimum and exponential function; uppercase boldface letters for vectors; |X| is the size of set X; 1 (x) is the indicator function which is equal to 1 if x > 0 and 0 otherwise; log (.) is the logarithm function; F X (x) and f X (x) being the cumulative distribution function (CDF) and the probability density function (PDF) of random variable (RV) X;

.

f (x) is the first-order derivative of function f over x, . f (x) = df (x) /dx.

II. SYSTEM MODEL A. LoRa Networks Modeling

Let us consider an uplink single gateway LoRa networks where the gateway is located at the center of the disc with area A = πR 2 , here, R is the network radius and a number of EDs which follow Poisson cluster processes (or doubly Poisson point processes), with density λ C > 0, inside A and λ C = 0, otherwise. In particular, in this work we consider two notable Poisson cluster processes, i.e., the Matérn and the Thomas cluster process with corresponding density λ Mat = λ P C Mat and λ Tho = λ P C Tho , where λ P is the density of the parent point process and C u , u ∈ {Mat, Tho}, is the average number of offspring per cluster in the u cluster process [START_REF] Haenggi | Stochastic Geometry for Wireless Networks[END_REF]. The offspring of the Matérn cluster process are independently and uniformly distributed in a disc of radius r Mat around the parent point while the offspring of the Thomas cluster process are scattered with variance σ 2 Tho around each parent point.

B. Channel Modelling

Let us consider a generic signal from an arbitrary ED to the gateway, it is impaired by both the small-scale fading and large-scale path-loss. It should be emphasized that the impact of the shadowing is implicitly studied by modifying the density of the EDs [START_REF] Tu | On the Energy Efficiency of Heterogeneous Cellular Networks with Renewable Energy Sources -A Stochastic Geometry Framework[END_REF].

1) Small-scale fading: Let us denote h q is the fading from an arbitrary node q to the gateway which follows Nakagami-m distribution with corresponding shape and spread parameters, i.e., m ≥ 1/2 and Ω. There is no doubt that Nakagami-m fading is one of the most general fading distributions that can represent other well-known distributions by properly adjusting its shape parameter, e.g., m = 1, is the Rayleigh fading and m = 1/2, is the single-sided Gaussian distribution. In addition, assuming that time is slotted (slotted ALOHA medium access [START_REF] Beltramelli | LoRa beyond ALOHA: An Investigation of Alternative Random Access Protocols[END_REF]) and the fading remains constant during one time-slot and changes between time-slot.

2) Large-scale path-loss: Consider a transmission link from a generic node q to the gateway, the large-scale pathloss is formulated as

ρ q = K 0 r β q , (1)
where β > 2 and K 0 = (4πf c /c) 2 are the path-loss exponent and the path-loss constant, respectively; f c is the carrier frequency and c = 3 × 10 8 (in meter per second) is the speed of light; and r q is the distance from node q to the gateway.

C. Power Consumption Modelling

Assuming that the ED is operated either in transmission mode or sleep mode. To be more precise, the ED is considered in transmission mode providing that it transmits packets to the gateway and in sleep mode otherwise. The power consumption under transmission mode comprises of two parts, i.e., the transmit power, P tx , and the static (circuit) power, P sta ; while in the idle mode, the ED consumes P idl power. Here, the static power consumption comprises of all other power consumption excluding the transmit power in the transmission mode, i.e., wake up, radio preparation, wait and receive the 1st and 2nd windows, radio off and so on. Furthermore, we also consider the impact of the hardware impairment on the static and idle power; indeed, the static and idle power are impaired by an additive noise which follows either Gaussian or uniform distribution. As a result, the practical static and idle power denoted as P sta and P idl are written as follows:

P sta = P sta + ω sta , P idl = P idl + ω idl , (2)
where ω x , x ∈ {sta, idl}, is the random variable (RV) and P sta and P idl are the ideal static and idle power.

III. PROBLEM FORMULATION

In this work, the principal objective is to identify the optimal transmit power denoted as P * tx that maximizes the energy efficiency (in bits/Joule) of whole networks given a set of input parameters, i.e., the fading parameters included both shape and spread factor, the transmission bandwidth, the spreading factor, the coding rate, the inter-arrival time between two packets, T in , the packet length, L pac , the path-loss exponent, the density of the PCP, λ C , and the imperfection static and idle power as shown in Fig. 1. Particularly, the problem in Fig. 1 can be formulated as follows:

max

Ptx∈[P min tx ,P max tx ] EE = PSE Pcon = λ A BWlog 2 (1+γD)Pcov(γD) λ A (Ptx+P sta)+(λC -λ A )P idl . (3)
It is noted that the range of the transmit power in (3), without loss of generality, can go from zero to infinity, i.e., P min tx → 0 and P max tx → ∞; where PSE refers to the potential spectral efficiency (in bits/s/m 2 ) that measures number of bit successfully transmitted per unit area and is formulated as

λ A BWlog 2 (1 + γ D ) Pcov (γ D ); Pcov (γ D )
is the coverage probability of an arbitrary link and formulated as

Pcov(γ D ) =Pr SIR = S I S ≥ γ I , SNR = P tx S σ 2 ≥ γ D , ( 4
)
where SIR is signal-to-interference ratio and is computed as SIR = PtxS0 PtxIS = Ptxh 2 0 /ρ0 Ptx i∈Φ A \(0) h 2 i /ρi ; S 0 is the signals from the ED of interest to the gateway; I S is the aggregate interference from all active EDs except for the desired one; assuming that packets with different spreading factor are perfectly orthogonal, hence, there is no inter spreading factor interference at the gateway [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF]; h 2 0 and ρ 0 being the channel gain and large-scale path-loss of the ED of interest while h 2 i and ρ i are the channel gain and the large-scale path-loss from interferer i to the gateway; Φ A \ (0) is the set of active EDs except for the desired ED which density λ A = p A λ C under the considered area; where p A = 1 Tin Lpac Rbit being active probability. In this work, assuming that the length of the packet, L pac , is identical among all EDs and R bit is the bit rate which is computed as R bit = SF BW 2 SF CR [START_REF] Goursaud | Dedicated networks for IoT: PHY/MAC state of theart and challenges[END_REF]. SNR is the signal-tonoise ratio and is formulated as SNR = P tx S 0 /σ 2 where σ 2 = 10 (-174+NF+10 log 10 BW)/10 [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF] is the noise variance of the AWGN noise; NF is the noise figure (in dBm) at receiver. In (3), Pcon is the average power consumption of the whole networks measured in Watt/m 2 ; the first term of Pcon, λ A P tx + P sta , accounts for the power consumption under transmission mode and the remain term, (λ C -λ A ) P idl is under idle mode; where P x = P x + E {ω x }, x ∈ {sta,idl}, is the average dissipation of the static and idle power.

It is apparent that the most intuitive approach to solve the optimization problem in (3) is to compute the EE in closedform expression followed by solving a non-linear equation to obtain P * tx . However, as shown in the sequel, it is impossible to obtain the closed-form expression of the EE under the considered system model. Moreover, even the distribution of the distance from an arbitrary node to the gateway as well as the aggregate interference are also unfeasible to represent in the closed-form expressions. Especially, by direct inspection (3), it is evident that the framework of the Pcov is essential in order to compute the EE, let us formulate the Pcov as follows:

Pcov (γ D ) = ∞ s=σ 2 γD/Ptx F I S (s/γ I ) f S (s) ds, (5)
where f S (s) and F I S (x) are the PDF of the intended signal and the CDF of the aggregate interference. Nonetheless, the CDF of I S has never existed even with the simplest scenario, i.e., the EDs follow Poisson point process (PPP); the PDF of the intended signal, on the contrary, can be computed in the closed-form expression under some special cases, for example, the EDs follow PPP combined with Rayleigh fading and special value of the path-loss exponent. It, as a result, is unworkable to obtain the closed-form expression of both Pcov and EE even with the simplest case [START_REF] Tu | A New Closed-Form Expression of the Coverage Probability for Different QoS in LoRa Networks[END_REF].

As a result, in this paper, we are going to find the P * tx by using deep learning approach. Nevertheless, different to the conventional DL approach that requires enormous training data to fully train the ANNs or simply called the data-driven approach. In this manuscript, we synergistically combine the model-based with the data-driven approach to derive the optimal transmit power in (3). To do so, we split the training process into two phases instead of one as the conventional approach. In the first training phase, the ANNs is trained by

P (γ D ) = M -1 (γ I ) -α (1 -exp (-M min {(Aγ I ) α , 1})) + exp (-M) min {A α , 1} -(γ I ) -α 1 (Aγ I -1) (6)
.

P (x) = . A (x)α(A (x)) α-1 (exp (-M(A (x) γ I ) α ) 1 (1 -A (x) γ I ) + exp (-M) 1 (A (x) γ I -1) 1 (1 -A (x))) (8)
the model-based data which may be simplified and inaccuracy but can bring some expert knowledge and useful for the second phase. In the second training phase, the pre-trained ANNs is re-trained by utilizing a small amount of practical data which are generated via Monte Carlo simulation. In next section, the design of the proposed neural networks as well as the detail description of the training process are provided. IV. DESIGN OF THE NEURAL NETWORKS In this section, the proposed ANNs networks is designed from creating the data to selecting the optimal architecture of the ANNs as well as identifying the performance metrics which are used to evaluate the performance of the ANNs. A. Generate data set 1) First training phase (pre-train phase) data set: It is obvious that data is the most important element of any neural networks. Under the considered system model, it is unreasonable to obtain a large number of data from either Monte Carlo simulations due to the resources constraint or the frameworks owing to the mathematical intractability as discussed in Section III. Consequently, in this section, we are going to simplify the system model in Section II so that the closed-form expression of the EE and the optimal transmit power, P * tx , can be numerically obtained. Particularly, the simplified system model is described as follows: i) assuming that the ED is no longer followed PCP but PPP with density λ P ; ii) the ideal static and idle power at EDs; and iii) the aggregate interference is approximated by the dominant interferer [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF] and the instantaneous fading by its average [START_REF] Tu | A New Closed-Form Expression of the Coverage Probability for Different QoS in LoRa Networks[END_REF]. Other assumptions are remained the same as Section II. Based on these assumptions, the approximated coverage probability denoted as P (γ D ) is computed in [START_REF] Pratt | Discriminability-based transfer between neural networks[END_REF] at the top of this page. Here, α = 2/β, M = p A λ P πR 2 , A = mθPtx R β K0σ 2 γD and 1 (x) is the indicator function. Having the approximated Pcov in hand, the optimal transmit power under the simplified system model denoted as P * ,1 tx , is obtained by solving following non linear equation:

.

P (P tx ) Q (P tx ) -P (P tx ) . Q (P tx ) = 0, (7)
where .

P (x = P tx ) is the first-order derivative of the approximated Pcov respect to the transmit power and is computed in [START_REF] Su | Energy Efficient Uplink Transmissions in LoRa Networks[END_REF] at the top of this page and . Q (x) = p A λ P is the firstorder derivative of the average power consumption respect to P tx ; Q (x) = λ A P tx + P sta + (λ C -λ A ) P idl and .

A (x = P tx ) = dA/dx = mθ R β K0σ 2 γD , respectively.
2) Second training phase (re-train phase) data set: The data in this phase are created imperatively by employing Monte Carlo simulation. However, owing to the resources constraint, a limited amount of data are created compared to the previous phase.

B. Data normalization

The data set is normalized before putting into the ANNs. In this work, the simple max-min normalization is applied as follows:

x = x ori -min (x ori ) max (x ori ) -min (x ori ) , (9)
where x ori and x ∈ [0, 1] are the input and output of the normalization process. In the sequel, without explicit explanation, we assume that the data set has already been normalized.

C. Networks Architecture

To maximize the performance of the ANNs, optimizing the networks architecture, i.e., the optimal number of hidden layers and/or number of neurons of each layer, is essential along with other hyper-parameters optimization. As the considered ANNs will be trained two times; it, in theory, exists two optimal typologies, the first one which merely optimizes from the 1st training set and another which optimizes from both the 1st and 2nd training set. In mathematics, the optimal architecture based on the 1st training set (partial optimum) and from both 1st and 2nd training set (full optimum) denoted as N *

1 and N * 1+2 as a function of the number of hidden layers, L ∈ {L min , . . . , L max }, the number of neurons of all hidden layers are formulated as follows:

N * 1 (L * 1 , N * 1 ) = min L L, N|Ψ (1)
N * 1+2 (L * 1+2 , N * 1+2 ) = min L L, N|Ψ (1)
+ Ψ (2) , [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF] where N = [N 1 , . . . , N L ] being a vector which contains number of neurons of all hidden layers; N l ∈ {N min , . . . , N max }, l ∈ {1, . . . , L}, is the number of neurons in l-th hidden layer; Ψ (u) , u ∈ {1, 2}, is the training set of u phase; L (.) is the objective function which needs to be minimized; (L * a , N * a ), a ∈ {1, 1+2}, is the optimal solution based on u training set; {L min , . . . , L max } and {N min , . . . , N max } are the minimum and maximum of number of hidden layers and number of neurons of each layer, respectively.

It is evident that optimizing DL topology in [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF] is always a cumbersome task as it requires not only the skills and experience of the user, but also the mastery of the data features [START_REF] Talbi | Optimization of deep neural networks: a survey and unified taxonomy[END_REF]. As a result, the grid search method is applied to find out the optimal architecture with the help of [START_REF] Subramanian | Deep Learning with PyTorch: A practical approach to building neural network models using PyTorch[END_REF]. In addition, to simplify the search space, we further assume that the number of neurons of all hidden layers is identical, i.e., N = N . Thus, the optimization problem in [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF] can be re-written as follows:

N * 1 (L * 1 , N * 1 ) = min L L, N |Ψ (1) N * 1+2 (L * 1+2 , N * 1+2 ) = min L L, N |Ψ (1) + Ψ (2) , (11)
In fact, the search space has massively reduced from T = (Lmax-Lmin) i=1

(N max -N min ) i in [START_REF] Cho | Improving TDoA Based Positioning Accuracy Using Machine Learning in a LoRaWan Environment[END_REF] to T = (L max -L min ) (N max -N min ) in [START_REF] Haenggi | Stochastic Geometry for Wireless Networks[END_REF].

Remark 1: It is apparent that the full optimum will theoretically provide better performance compared to the partial optimum. Thus, the main purpose is to clarify the gap between the partial optimum versus the full optimum. In case the gap is relatively small, it is more beneficial to employ the partial optimum than its counterpart provided that the consumed resources are taken into consideration.

D. Loss function

As the considered problem belongs to the regression one, the typical mean square error is deployed as the loss function of the ANNs and is formulated as follows:

min MSE = (|Ψ|) -1 j∈Ψ (y j -y j ) 2 , ( 12
)
where Ψ and |Ψ| are the training set and its size; y is the observed output and y is the predicted output.

E. Performance Metric

The mean square error and the R squared, R 2 , or the coefficient of determination are utilized as the main metrics to evaluate the performance of the ANNs. As the MSE measures the average squared difference between the estimated and true values and is provided in [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF]; the R 2 , differently, measures how close of two sets of data in terms of the distribution. Mathematical speaking, R 2 can be formulated as follows:

R 2 = 1 - SS res SS tot = 1 - j∈Ψ (y j -y j ) 2 j∈Ψ (y j -y) 2 , (13)
where SS res and SS tot are the residual sum of squares and the total sum of squares. Here, SS res measures the amount of variability that is left unexplained after performing the regression while SS tot measures the total variance. As a consequence, if R squared towards 1, it explicitly means that most of the variability are explained by the model; y = (|Ψ|) -1 j∈Ψ y j is the mean of the observed output.

F. Training process

Considering a fully-connected feedforward neural networks with L hidden layers, N neurons per hidden layer, one input layer with I input and single output as shown in Fig. 2. The input of the ANNs are the path-loss exponent, the fading parameters, the bandwidth, the coding rate, the spreading factor, the packet length, the inter-arrival time, the density of EDs, the practical static and idle power. The single output is the optimal transmit power, P * tx , which is obtained by [START_REF] Beltramelli | LoRa beyond ALOHA: An Investigation of Alternative Random Access Protocols[END_REF] in the pre-trained phase or Monte Carlo simulation in the re-trained phase.

Given a training set, Ψ (u) , u ∈ {1, 2}, the training process is commenced with normalizing the training set by using (9), followed by forward propagation. Particularly, let us denote x l n as the output signal at neuron n, n ∈ {1, . . . , N }, in hidden layer l, l ∈ {1, . . . , L}, and is computed as

x l n = ζ V i=1 w l i,n x l-1 i +b l n , V = V = I, l = 1 V = N , l = 1 , ( 14
)
where w l i,n is the weight from node i in layer l -1 to node n in layer l; x l-1 i is output of neuron i in layer l-1; when l = 1, we have x 0 i , i ∈ {1, . . . , I}, is the i-th input of the ANNs; tx is obtained via Monte Carlo simulation which absolutely requires huge resources and efforts, thus, a limited search space is considered. P * ,1 tx , contrarily, is obtained by maths which is effortless compared to Monte Carlo method thus the range can be arbitrarily meaningful numbers.

In addition, in the second training set, Ψ (2) , we have the average number of offspring per cluster, C u = {2, 6}, u ∈ {Tho, Max}, the radius of each cluster in Matérn cluster process, r Mat = {2, 6} m, and the standard derivation in Thomas cluster process, σ 2 Tho = {0.5, 3}. Moreover, also in the second training set, Ψ (2) , both the static and idle power are impaired by an additive noise which follows either Gaussian distribution with zero mean and unit variance (in dBm) or 1) , and retrained with practical data, Ψ (2) ; the red curves are the conventional approach which directly trains with Ψ (2) . The ANNs architecture is (L = 4, N = 55); |Ξ| = 500; Matérn cluster process with C Mat = 2, r Mat = 2 and ωx, x ∈ {sta,idl}, follows Gaussian distribution.

uniform distribution from [-0.5, 0.5] (dBm). The size of the 1st training set is equal to 30000, i.e., |Ψ (1) | = 30000, while the size of the 2nd training set is less than 2000, i.e., |Ψ (2) | ≤ 2000, and the size of the test set denoted as |Ξ| is always fixed at 500, i.e., |Ξ| = 500. Furthermore, the test set is solely available for the second training phase. It means that we do not evaluate the performance of the ANNs after the first training. It should be noted that the training set in both phases are able to cover most of the practical environment, i.e., urban or rural area, as well as the different applications/end-devices which are applied into various domains, i.e., smart home, smart city and smart health, etc [START_REF]Machine to Machine (M2M) Evaluation Methodology Document (EMD)[END_REF]. In this section, without explicit stated, the performance from the best epoch is chosen and the size of the training data indicates the second training phase, |Ψ

|.

A. Direct training versus double training approach

In this section, the performance between the conventional approach (direct training with the second data set, Ψ (2) ) and the proposed approach are examined. Fig. 3 illustrates the mean square error of both the training and test set in second training phase versus number of epochs of the proposed approach and the conventional approach. To be more specific, the curves from the double training approach are firstly trained by 30000 model-based data, i.e., Ψ (1) = 30000, in 150 epochs then re-training by either Ψ (2) = 500 or Ψ (2) = 2000, in another 150 epochs. The conventional approach, on the contrary, is direct training with the same amount of data of Ψ (2) in 300 epochs. We observe that our approach outperforms the conventional one in both training set and test set. Particularly, our approach is better almost ten times compared to its counterpart. In addition, it is expected that increasing training set slightly ameliorates the MSE performance under the proposed method as its weights and bias have already been configured to the near optimal values after the 1st training phase. The conventional one, differently, improves dramatically by soaring number of 2) ; each point are chosen from the best epoch of that set up. Solid curves are conventional approach and dash dot curves are the proposed approach. Thomas cluster process with C Tho = 2, σ Tho = 0.5 and the ideal static and idle power, 4 hidden layers, i.e., L = 4, and |Ξ| = 500.

training data. Fig. 3 also reveals that increasing the number of epochs is not a wise choice to enhance the performance of the MSE especially in direct training approach where the MSE start constant from around 50 epochs (case Ψ (2) = 500). Moreover, we experience that 300 and 150 epochs are sufficient to obtain the stable results for both the conventional and proposed approach. As a result, in the sequel, this number of epochs are yielded if no explicit stated. Fig. 4 shows the performance of the MSE and R 2 of both the conventional and the proposed methods versus the number of neurons, N , with different values of |Ψ (2) |. Each point of all curves are chosen from the best epoch that provides the minimum MSE or the maximum R 2 among all epochs. The results are measured on the test set and it should be clarified that the best epoch over the test set is not necessary the same over the training set. We see that the performance of the conventional approach is constant with N and is regardless of the metrics, i.e., MSE and R 2 . It can be interpreted that the higher number of neurons, N , is not necessary the better performance under the direct training approach. The double training approach, on the contrary, is fluctuated and the largest number or neurons, N = 70, is not necessary the best performance, the best one, however, is case N = 55 where the MSE = 1.236×10 -4 and R 2 = 0.9953 for case Ψ (2) = 2000. It is interesting to point out that although increasing number of neurons does not always benefit, one is always true that the larger the training set the better the performance of the ANNs. Indeed, by rising Ψ (2) from 500 to 2000, the R 2 of the conventional approach improves from 0.9252 to 0.9926 and is better than the R 2 of proposed approach case Ψ (2) = 500. Nevertheless, the MSE of the double training approach is never worse than the direct training one. Even taking the worst case from the proposed method into comparison, i.e., Ψ (2) = 500, its performance is still better than the best one of the conventional approach, i.e., Ψ (2) = 2000. In addition, in Fig. 4(b), we also experience that the R 2 of the double training one is never lower than 0.985. It means that only (2) ; solid lines and dot lines are corresponding to Ψ (2) = 500 and Ψ (2) = 1000, respectively; |Ξ| = 500; the curves 'conventional' are direct training approach and selected from the best architecture based on Ψ (2) ; the curves 'partial optimum' and 'full optimum' are proposed approach and selected from the best architecture based on the Ψ (1) and Ψ (1) + Ψ (2) , respectively. Matérn cluster process with C Mat = 6, r Mat = 2 and ωx, x ∈ {sta,idl}, follows uniform distribution.

around 1.5% the variability of the test set can not be explained by the ANNs. B. Partial vs. Full optimum architecture

In the following, we are going to study the performance of the proposed approach by measuring the performance of the partial optimum, N * 1 , versus full optimum architecture, N * 1+2 . The partial optimum architecture means that we select the best topology based on Ψ (1) , and solely utilize this architecture for Ψ (2) . It is re-emphasized that the selection is based only on the training set. The full optimum architecture, N * 1+2 , selects the architecture that provides the best performance from both Ψ (1) and Ψ (2) . Fig. 5 unveils the performance of the MSE in both training set and test set of the second phase versus the number of epochs with various optimum architecture, i.e., partial and full optimum. In Fig. 5, in order to have a fair comparison, the curves from the conventional approach is also chosen from the best architecture in the same search space as the proposed approach, i.e., T = (L max -L min ) (N max -N min ). The conventional approach still underperforms both optimum architecture from the double training approach though the best architecture has been selected. In Fig. 5(a), it is crystal clear that the full optimum architecture is superior compared to the partial one. Nevertheless, the gap between these curves, in general, are negligible, for example, 3.99 × 10 -5 vs. 6.79 × 10 -5 for Ψ (2) = 2000, and 1.34 × 10 -4 vs. 1.19 × 10 -4 for |Ξ| = 500. In addition, the expenses for this superior performance is quite expensive in terms of the consumed resources, namely, in order to attain the best architecture under full optimum approach, each topology requires to train exactly 300 epochs like the conventional approach while the partial optimum approach requires only half of each, i.e., 150 epochs. Mathematical speaking, assuming that the search space has T = (L max -L min ) (N max -N min ) architecture, then both the conventional and full optimum approach need 300T epochs to attain the best architecture while the partial optimum (2) ; the curves 'conventional' are direct training approach and selected from the best architecture based on the Ψ (2) ; the curves 'partial optimum' and 'full optimum' are proposed approach and selected from the best architecture based on the Ψ (1) and Ψ (1) + Ψ (2) , respectively. Thomas cluster process with average C Tho = 6 and σ Tho = 0.5 and ωx, x ∈ {sta,idl}, follows uniform distribution. approach, the epochs' requirements are only 150 (T + 1). Again, we experience that surging the training data will monotonically increase the performance of the ANNs. As a consequence, in next figure, we investigate in detail the impact of training data on the performance of MSE and R 2 .

Fig. 6 illustrates the performance of R 2 (a) and MSE (b) versus |Ψ (2) |. We observe the same pattern as Fig. 4 that, in general, increasing training data will monotonically improve R 2 and decline the MSE of all approaches, i.e., conventional, partial and full optimum. However, the pace of the improvement among these methods are different, namely, the direct training one enhances almost 10 times in terms of the MSE when training data rises from 100 to 2000 while for both partial and full optimum approach, it merely increases around 3 times. As for the R 2 , the performance of the partial and optimum approach are nearly stable and close to one while the remain one boosts around 5% in case the training data goes from 100 to 2000. From Fig. 6, it is apparent that in order to obtain the MSE equal to 10 -4 , 2000 training data is sufficient for both optimum architecture. Furthermore, this figure also confirms the above statement that the gap between the partial and optimum approach is negligible. Fig. 7 shows the optimal transmit power, P * tx , of three different ways versus the test set where the samples are sorted in ascending order. Looking at this figure, there is no surprising that the direct training approach experiences the biggest error between the observed output and the predicted one while the error under the partial and full optimum are relatively small. To better measure the error between the observation and the predicted one, let us utilize the mean absolute error (MAE) which refers to the arithmetic average of the absolute errors between the observed and predicted values. Then the MAE under the conventional, the partial and the full optimum approaches in Fig. 7 the 'no' curves denote without freezing any layers; the '1st4th' curves denote freezing from the 1st to the 4th hidden layers (all hidden layers); the 'in3rd' curves mean freezing from the input layer to the 3th hidden layer; the '2ndout' curves mean freezing from the 2nd hidden layer to the output layer. Thomas cluster process with C Tho = 6 and σ Tho = 0.5 and ωx, x ∈ {sta,idl}, follows Gaussian distribution.

i.e., 0.046, 0.009 and 0.008 for the conventional, partial and full optimum architecture.

C. Re-train whole vs. parts of the ANNs Fig. 8 illustrates the MSE in both training set and test set versus |Ψ (2) | by freezing some layers of the considered ANNs during the second training phase. In particular, we consider four distinct scenarios; i) re-train whole ANNs or no freezing denoted as 'no'; ii) freezing all hidden layers denoted as '1st4th'; iii) freezing from the input layer to the 3rd hidden layer denoted as 'in3rd'; and iv) freezing from the 2nd hidden layer to the output layer denoted as '2ndout'. In the last 3 cases, the ANNs only re-trains 2 layers, e.g., the input and output layers for case '1st4th', the input and the 1st hidden layer for case '2ndout' and the last hidden layer (the 4th hidden layer) and the output layer for case 'in3rd'. From the figure, we observe a fairly large gap between case 'in3rd' and others, in fact, its performance is worse than almost 10 times in both training and test set. The second worst case is '1st4th' where the ANNs only updates the weights and bias of the input and output layers. The best scenario is, of course, the case without freezing any layers. However, the difference between the best and the second best case, i.e., case '2ndout', is significantly small, especially when the training size is adequately large, i.e., around 2000 data. From Fig. 8, we conclude that the impact of freezing the input layer and/or hidden layers close to input on the performance of the ANNs is more serious than ones close to the output layer. Moreover, by carefully freezing some layers, the performance of the ANNs slightly reduces compared with fully re-trained the networks, thus, it is feasible to freeze some layers in order to save time and resources when re-training the neural networks.

VI. CONCLUSION

In this paper, the maximization of the EE respect to the transmit power is studied by combining the advantages of both model-based in wireless networks and the data-driven of the deep learning technique. Our findings show that the proposed training approach outperforms the conventional one almost 10 times under some scenarios. Moreover, the re-training process can be done either in whole or parts of the ANNs. Finally, the application of the proposed approach is general to cover most of EDs distributions as well as to overcome the hardware impairment at EDs.

Fig. 1 .

 1 Fig. 1. Schematic of the considered artificial neural networks

Fig. 2 .

 2 Fig. 2. The fully-connected feedforward neural networks

  ζ () is the activation function; b l n denotes the bias of node n in layer l. In this work, we utilize the sigmoid activation function as the training set is normalized to [0, 1]. The predicted output of the ANNs denoted as y is linearly combined the output of the last hidden layer with the weight from the last hidden layer to the output as y = N +1 i=1 w L+1 i x L i + b y , b y is the bias of the output. The forward pass ends by computing the loss function given in (12). The next step is backward propagation where we compute the gradient of the loss function respect to all the weights and bias of the networks. This gradient is then used to update the weights and bias to minimize the loss function by deploying Adam algorithm with adaptive learning rate [16] and we finish one training epoch. The training process is repeated by re-computing the forward and backward propagation; updating the loss function, weights and bias until reaching the stopping conditions, i.e., the maximum epoch or the MSE is smaller than a predefined threshold or the MSE starts increasing, etc. Remark 2: In the step back propagation during the second training phase, it is unnecessary to update all the weights and bias. In fact, the ANNs can re-train partially instead of fully re-trained or the weight and bias of some layers can be kept constant during the second training phase. V. NUMERICAL RESULTS AND DISCUSSION In this section, numerical results are provided to confirm our findings. Unless otherwise stated, following values are used to generate the input training set in the first and second phases: β ∈ [2, 6], m ∈ [0.5, 30], Ω ∈ [1, 50], SF ∈ {7, . . . , 12}, BW ∈ {1, . . . , 500} KHz, CR = 4/(4 + o), o ∈ {1, . . . , 4}, L pac ∈ {1, . . . , 120} (bytes), T in ∈ {1, . . . , 90} (minutes), P cir ∈ [0.5, 5] (dBm), P idle ∈ [-10, -5] (dBm), λ P = N /πR 2 where N ∈ {500, . . . , 100000} and R = 2000 m; γ I = 1 dBm, NF = 6 dBm, f c = 868 MHz; γ D is chosen from one of the following values: γ D ∈ {-6, -9, -12, -15, -17.5, -20} (dBm), which absolutely depends on the utilized SF, e.g., if SF = 7 then γ D = -6 dBm and SF = 12 then γ D = -20 dBm, respectively. The range of the output, P * tx , in the first and second phases, on the other hand, are different. Particularly, the range of the optimal transmit power in the first training set denoted as P * ,1 tx is from -100 dBm to 60 dBm while the range of P * tx in the second training set denoted as P * ,2 tx , is only from -40 dBm to 50 dBm. The chef reason behinds this difference is that P * ,2

Fig. 3 .

 3 Fig. 3. Mean square error (MSE) versus number of epochs of training set (a) and test set (b) in the second training phase; the blue curves are proposed approach where the ANNs is pre-trained by model-based data, Ψ(1) , and retrained with practical data, Ψ(2) ; the red curves are the conventional approach which directly trains with Ψ(2) . The ANNs architecture is (L = 4, N = 55); |Ξ| = 500; Matérn cluster process with C Mat = 2, r Mat = 2 and ωx, x ∈ {sta,idl}, follows Gaussian distribution.

4 Blue: |Ψ ( 2 )Fig. 4 .

 424 Fig.[START_REF] Zappone | Model-Aided Wireless Artificial Intelligence: Embedding Expert Knowledge in Deep Neural Networks for Wireless System Optimization[END_REF]. MSE (a) and R 2 (b) vs. number of neurons of each hidden layer, N , of the test set with different training data size, i.e., Ψ(2) ; each point are chosen from the best epoch of that set up. Solid curves are conventional approach and dash dot curves are the proposed approach. Thomas cluster process with C Tho = 2, σ Tho = 0.5 and the ideal static and idle power, 4 hidden layers, i.e., L = 4, and |Ξ| = 500.

Fig. 5 .

 5 Fig. 5. MSE vs. number of epochs in training set (a) and test set (b) of Ψ(2) ; solid lines and dot lines are corresponding to Ψ(2) = 500 and Ψ (2) = 1000, respectively; |Ξ| = 500; the curves 'conventional' are direct training approach and selected from the best architecture based on Ψ(2) ; the curves 'partial optimum' and 'full optimum' are proposed approach and selected from the best architecture based on the Ψ(1) and Ψ(1) + Ψ(2) , respectively. Matérn cluster process with C Mat = 6, r Mat = 2 and ωx, x ∈ {sta,idl}, follows uniform distribution.

Fig. 6 .

 6 Fig.6. R 2 (a) and MSE (b) versus Ψ(2) ; the curves 'conventional' are direct training approach and selected from the best architecture based on the Ψ(2) ; the curves 'partial optimum' and 'full optimum' are proposed approach and selected from the best architecture based on the Ψ(1) and Ψ(1) + Ψ(2) , respectively. Thomas cluster process with average C Tho = 6 and σ Tho = 0.5 and ωx, x ∈ {sta,idl}, follows uniform distribution.

Fig. 7 .Fig. 8 .

 78 Fig.7. P * tx versus the test set (sorted in ascending order) under various case studies, e.g., conventional (a), partial optimum (b) and fully optimum (c); the 'Observed' curves mean the observation output and the 'predicted' curves mean the predicted output from the ANNs. Matérn cluster process with C Mat = 6, r Mat = 2 and ωx, x ∈ {sta,idl}, follows uniform distribution; Ψ(2) = 500, and |Ξ| = 500.

  are 4.14, 1.47 and 1.35 dBm, respectively. It is noted that this MAE is computed based on the original range of the P

* tx . Hence, if we compute the MAE after normalizing by (

9

), the MAE of these approaches are smaller,


Mikhail Kamalov
email: mikhail.kamalov@inria.fr

Konstantin Avrachenkov
email: k.avrachenkov@inria.fr

GenPR: Generative PageRank framework for Semi-Supervised Learning on citation graphs

Keywords: semi-supervised learning, generative model, PageRank, citation graphs, neural networks

Nowadays, Semi-Supervised Learning (SSL) on citation graph data sets is a rapidly growing area of research. However, the recently proposed graph-based SSL algorithms use a default adjacency matrix with binary weights on edges (citations), that causes a loss of the nodes (papers) similarity information. In this work, therefore, we propose a framework focused on embedding PageRank SSL in a generative model. This framework allows one to do joint training of nodes latent space representation and label spreading through the reweighted adjacency matrix by node similarities in the latent space. We explain that a generative model can improve accuracy and reduce the number of iteration steps for PageRank SSL. Moreover, we show that our framework outperforms the best graph-based SSL algorithms on four public citation graph data sets and improves the interpretability of classification results.

Introduction

The main idea of SSL is to solve a classification task with an extremely low number n l of labeled data points in comparison with the number n u of unlabeled data points (n l n u ). Therefore, with regard to citation graphs with a huge amount of nodes (e.g. Pubmed, MS Academic) SSL is a good technique to avoid preparing data points for supervised learning. The area of SSL focusing on the classification of nodes in citation graphs, in particular, citation graphs is called a graph-based SSL. The standard input for graph-based SSL algorithms is a graph G = (V, E) with n = n l + n u = |V| nodes (papers), e = |E| edges (citations), A ∈ R n×n is the adjacency matrix and X = (x i,j ) n,d i,j=1 is a matrix of nodes where each node x i = (x i,1 , . . . , x i,d ) ∈ R d has a feature representation in d-space. In the context of citation graphs X is a bag-of-words representation for the nodes (papers). Moreover, each node belongs to one of c classes {C 1 , . . . , C c }. Also we have the labels matrix Y = (y i,j ) n,c i,j=1 ∈ R n×c such that y i,j = 1 if x i ∈ C j and y i,j = 0 otherwise. Nowadays, the area of the graph-based SSL consists of two main research directions:

Supported by MyDataModels company.

the classical diffusion-based linear algorithms that spread the class information through the adjacency matrix A: Label Propagation (LP) [START_REF] Zhu | Learning from labeled and unlabeled data with label propagation[END_REF], PageRank SSL (PRSSL) [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]; the graph convolution-based neural network (NN) algorithms that apply the dot product of adjacency matrix A with NN nonlinear transformation of features X for the classification. The recently proposed: approximated Personalized graph NN (APPNP) [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF], Graph Attention Network (GAT) [START_REF] Veličković | Graph attention networks[END_REF], Graph Convolution Network (GCN) [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF].

Regarding graph-based SSL algorithms, one can notice that they use A with binary weights on edges (citations), which can cause a loss of information about node similarities and further may negatively affect label diffusion through A. We address this issue through the following contributions:

-We propose a novel graph-based SSL inductive (I)/ transductive (T) framework, created by embedding PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF] in generative model (GenPR); -We show that the generative model can be used to reweight A to further improve PRSSL label spreading; -We show that GenPR improves the interpretability of NN classification results based on the information about nodes similarity in the latent space; -We show that GenPR outperforms the recently proposed algorithms for graph-based SSL and reduces the number of steps of PageRank [START_REF] Page | The pagerank citation ranking: Bringing order to the web[END_REF] to obtain more accurate classification results.

Related work

Our framework is based on the combination of the following two ideas:

1. PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF] gives a PowerIteration based explicit solution for the graph classification:

F t = αD -σ AD σ-1 F t-1 + (1 -α)Y ; t ≥ 0
where F t is a result of the t-th iteration and α is a regularization parameter in the range [0, 1] and σ is a power of D i,i = n j=1 A i,j ; 2. generative semi-supervised model (M2) [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]

: p(x, ŷ, z) = p(x|z, ŷ)p(ŷ)p(z)
where p(z) = N (z|0, I), p(ŷ) is a categorical distribution of latent class variable and p(x|z, ŷ) is a nonlinear transformation of the latent variables z and ŷ.

Since our framework is the graph convolution-based NN algorithm with PageRank, we need to define the main difference with APPNP [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF]. The difference is that GenPR jointly trains a redefined generative model [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF] and PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF] with a linear combination of A and similarity matrix in latent space, while APPNP applies PageRank with default A as a preprocessing step for output of multilayer perceptron (MLP).

3 Generative PageRank (GenPR)

Intuition of GenPR

Before we go into details of our framework let us define the motivation and the intuition behind GenPR. The main idea of GenPR is to resolve the following issues:

1. A i,j = 1 does not provide the information about impact of cited paper j on the citing one i; 2. A i,j = 0 may show that author i did not cite the paper j, but he could have used some information from it.

Let us define some useful notation for the GenPR intuition:

x i ∈ X is a i.i.d. samples of some continuous random variable x, then an output of MLP Y * = (y * i ) n i=1 ∈ R n×c is a sample from random variable y * given x as an input; Z = (z i ) n i=1 ∈ R n×d where z i is a latent representations of each node x i sampled from latent random variable z in d -space; W = (w i,j ) n i,j=1 ∈ R n×n is a similarity matrix where each element w i,j = h(z i , z j ), ∀z i , z j ∈ Z is an output of some positively defined kernel h; A = A + γW is a reweighted adjacency matrix A with a parameter γ ∈ [0, 1] for W ; D i,i = n j=1 A i,j is a diagonal matrix.
Let us redefine the recurrent formula of PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF] using Y * at each training epoch as a replacement of real labels Y = Y * :

F t = αD (-σ) A D (σ-1) F t-1 + (1 -α)Y * ; (1)
where F t = (y pr i,j ) n,c i,j=1 ∈ R n×c . Here y pr i = (y pr i,1 , . . . , y pr i,c ) is a sample from random variable y pr since (1) is a transformation of the random variable y * and

F 0 = Y * .
Then assume that F t will improve the accuracy of Y * by using the information of nodes similarity in latent space during the t-th iterations. We named it the PageRank spreading assumption. Moreover, we propose to use Y * as a new labels. Let us notice that and y * ∼ y pr due to the PowerIteration PageRank property

||F t -Y * || 1 ≤ 1 1-α ||F 1 -Y * || 1 [2](Property 12
). This allows us consistently use the aforementioned PageRank spreading assumption in training process of the generative model:

p(x, y * , z) ≈ p(x|z, y pr )p(y pr )p(z) (2)
where p(•) is a PDF of a random variable.

Objective function of GenPR

In this subsection we consider the inductive regime of GenPR which allows us to train jointly the generative model (2) and PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]. Now let us define GenPR objective function. It is obtained by maximizing the variational lower bound of the data log-likelihood of (2) with variance φ and generative θ parameters [START_REF] Kingma | Auto-encoding variational bayes[END_REF]:

log p(x, y * ) ≥ E q φ (z|x,y * ) log p θ (x|z, y pr ) + E q φ (z|x,y * ) log p θ (y pr ) -D KL (p(z)||q φ (z|x, y * )) (3)
where q φ (z|y * , x) = N (z|µ(y * , x), σ 2 (x)) is a multivariate Gaussian distribution parameterized by µ(y * , x) and σ(x) that are inferred from NN layers for expectation and variance respectively; p θ (x|z, y pr ) = f θ (z, y pr ) is a nonlinear transformation of z and y pr by NN layer; p θ (y pr ) = P R(y * , µ(y * , x), A) is a linear transformation of y * by (1) (the NN layer version will be defined in the next subsection 3.3), p(z) = N (z|0, I) is a multivariate Gaussian distribution and D KL (•||•) is the Kullback-Leibler divergence.

Since we can trade the quality generation of x for the quality of y pr i and estimate y pr i using the information from n l , we can use β ∈ [0, 1] as a weight parameter for p θ (x|z, y pr ) and the categorical crossentropy U(F t , Y ) = n l i=1 c j=1 (y i,j • log(y pr i,j )) for y pr i estimation. Thus, we obtain from (3) the final inductive (I) GenPR objective function:

L(θ, φ, x, Y ) = βE q φ (z|x,y * ) log p θ (x|z, y pr ) + log p θ (y pr ) -D KL (p(z)||q φ (z|x, y * )) -U(F t , Y ) (4)
The difference between inductive (I) and transductive (T) regimes of GenPR is that transductive GenPR does not use the proposition that y * is a new labels and an objective function looks as follows:

L T (θ, φ, x, Y ) = βE q φ (z|x) log p θ (x|z) -D KL (p(z)||q φ (z|x)) -U(F t , Y ) (5)
.

Fig. 1. Sample nodes from Citeseer data set: a -A before GenPR, where colored nodes are labeled and grey are unlabeled, straight black edges are citations between nodes (papers); b -A after GenPR, where all colored nodes are result from F t , and color of an edges by weights from A (cyan is a lower weights, maroon is a higher weights); cthe result of filtering lower weight edges for the node 545.

Architecture of GenPR

Since we have defined the objective function of GenPR (4) we can explain the GenPR layers architecture. The part of z inference contains the following layers:

Y * = π θ (X); π θ (X) = h 1 (XW 1 + B 1 ) (6) µ(X, Y * ) = h µ (concat(X, Y * )W µ + B µ ) (7) σ(X) = h σ (XW σ + B σ ) (8)
where h . and B . are activation functions and biases for NN layers respectively; W 1 ∈ R d×c , W µ ∈ R (d+c)×d and W σ ∈ R d×d are trainable weight matrices of MLP [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF], expectation [START_REF] Page | The pagerank citation ranking: Bringing order to the web[END_REF] and variance (8) for a NN layer respectively; (m i ) n i=1 = µ(X, Y * ) is an output of ( 7) layer with m i ∈ R d ; concat(•, •) is a matrix concatenation column-wise.

To avoid the issues with high variance of the gradient estimation of E q φ (z|x,y * ) log p θ (x|z, y pr ) by Monte Carlo method, we follow [START_REF] Kingma | Auto-encoding variational bayes[END_REF] in using the reparameterization trick to compute a low-variance gradient estimator for q φ (z|x, y * ):

q φ (z|x, y * ) ∼ Z, Z = µ(X, Y * ) + σ(X) , ∼ N (0, I)
where is an element-wise product and is a random variable. Now we can define (1) as a sequential sublayers in P R(Y * , µ(Y * , X), A):

1. the reweighting of A:

w i,j = h m i , m j ; ∀ w i,j ∈ W ; (9)
A = A + γW ; ( 10
)
where γ is a parameter of involvement W in reweighting of A within the range [0, 1]. Here we compute the similarities between the outputs of ( 7) because we assume that the expectation of the latent variable z more correctly defines the differences between nodes in latent space.

the regularization of

A : Â = D (-σ) A D (σ-1) ; D i,i = n j=1 A i,j (11)
where σ is a parameter for selection of regularization type: σ = 1 is a Standard Laplacian; σ = 0 is a PageRank; σ = 1/2 is a Normalized Laplacian; 3. the redefined PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]:

F t = α Â F t-1 + (1 -α)Y * ; t ≥ 0; ( 12
)
where F 0 = Y * (6) and F t is a result of the t-th iterations, smoothly changing the node labels Y * during iterations.

The final layer is the reconstruction of nodes (papers

) X = f θ (Z, F t ) where X ∈ R n×d : f θ (Z, F t ) = h 2 (concat(Z, F t )W 2 + B 2 ) ( 13
)
where W 2 ∈ R (d +c)×d , B 2 are weight and bias for x generation p θ (x|z, y pr ) = f θ (z, y pr ). We can turn to transductive regime of the aforementioned GenPR layers architecture by using modified loss as in [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]. The Figure 2 presents the difference between inductive (I) and transductive (T) GenPR architectures.

Experimental setup

For conducting an experiment in the graph-based SSL area we have taken following citation-graph data sets: Citeseer, Cora-ML, Pubmed, MS-Academic (MSA).

The description and statistics of data sets available in work APPNP [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF]. These data sets are available here1 . As a baseline algorithms we have considered: the classical diffusion -based linear algorithms (LP [START_REF] Zhu | Learning from labeled and unlabeled data with label propagation[END_REF], PRSSL [START_REF] Avrachenkov | Generalized optimization framework for graph-based semi-supervised learning[END_REF]); the recent graph convolution -based NN (GCN [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF], GAT [START_REF] Veličković | Graph attention networks[END_REF], APPNP [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF]); the deep generative model (M2 [START_REF] Kingma | Semi-supervised learning with deep generative models[END_REF]). To avoid overfitting issue we applied L 2 regularization with parameter λ = 0.05 for weights W . , dropout for  with rate dr = 0.5 at each PowerIteration step and learning rate l = 0.001 for Adam optimizer. Moreover, we have used the random train-test-validation splitting strategy described in [START_REF] Klicpera | Predict then propagate: Graph neural networks meet personalized pagerank[END_REF] and repeated experiments on each data set 500 times. For a fair model comparison we have made an architecture and parameters of GenPR that are very close to APPNP and GCN. In particular, for all data sets use the intermediate embedding layer f 0 (X) = relu(XW 0 + B 0 ) with W 0 ∈ R d× d as the input for (6) with d = 64, W 1 ∈ R d×c and h 1 (•) = sof tmax(•), d = 64 in ( 7) and ( 8), σ = 0.5 and t = 4 in (11), B . = 0. In (12) for MSA α = 0.8, for Cora-ML, Pubmed and Citeseer α = 0.9.

We have selected the specific parameters of GenPR by the 5 fold crossvalidation grid search2 . For all data sets use h(m i , m j ) = (m T i m j ) 3 in ( 9) and β = 0.001 in ( 4) and [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]. In particular, we have used: for Citeseer: γ = 1 in (10), 7) and ( 8

h µ (•) = h σ (•) = relu(•) in (

Experimental results

Table 1 presents performance of the classification based on the default adjacency matrix A or on the node features X leads to loss of classification quality because we do not use all available information. In the case of the combination of X and A, GenPR significantly and consistently outperforms the others due to the intuition that default A contains incomplete information about nodes similarity. Since we have reached the best results with GenPR (I) and γ = 1 for Citeseer, it means that latent information is helpful for reweighting default adjacency matrix A (citation graph). In particular, Figure 1 (c) shows that GenPR can be used for the explanation of classification results, by filter the edges by weight and observe nodes with more influence on considered one (e.g. node 545). Figure 3 shows the GenPR outperforms the APPNP not only in terms of accuracy, but also in number of PowerIteration steps, because GenPR takes less steps to converge for better accuracy than APPNP. Moreover, the GenPR is less complex than APPNP because it uses just one layer for MLP rather than 2 in APPNP.

Conclusion

In this work, we propose a graph-based SSL (I)/(T) framework created by embedding PRSSL in generative model. Based on the experimental results, we show that the generative model application for PRSSL can be used not only for the label spreading improvement, but also for interpretation of the classification results. We also show that GenPR significantly and consistently outperforms all other algorithms on every data set and requires less number of PageRank PowerIteration steps. Since GenPR produces complete weighted graph defined by A we can use PageRank properties to split A into batches that are complete subgraphs, which opens an opportunity to explore a distributed version of GenPR. The other interesting direction to investigate is an application of GenPR on data sets without default graph structure (e.g. images).

Fig. 2 .

 2 Fig. 2. The I-inductive (a) and T-transductive (b) architectures of GenPR.

  ), h 2 (•) = sigmoid(•) in (13); for Cora-ML, Pubmed and MSA: γ = 0.001 in (10), and h . (•) = linear(•).

Fig. 3 .

 3 Fig. 3. Average accuracy of GenPR (I) inductive, GenPR (T) transductive and APPNP over the t-iteration steps.

Table 1 .

 1 Average accuracy (%) on citation graphs. and denote the statistical significance (t-test) of GenPR for p < 0.05 and p < 0.01, respectively, compared to the APPNP.

	Input Data set	Citeseer Cora-ML Pubmed MSA
	A	PRSSL LP	71.21 45.32	78.12 68.31	72.51 63.12	76.12 65.32
	X	M2	70.81	79.22	77.6	86.12
		APPNP	75.74	85.09	79.71	93.28
		GAT	75.43	84.41	77.73	91.18
	X,A	GCN	75.31	83.52	78.65	92.09
		GenPR (I) 77.18	85.52	80.09	94.08
		GenPR (T) 76.91	86.19	81.13 93.81

https://github.com/klicperajo/ppnp/tree/master/ppnp/data

https://scikit-learn.org/stable/modules/grid search.html


Marco Faroni
email: m.faroni003@unibs.it

Manuel Beschi
email: manuel.beschi@itia.cnr.it

Antonio Visioli
email: antonio.visioli@unibs.it

Lorenzo Molinari
email: lorenzo.molinaritosatti@itia.cnr.it

Lorenzo Molinari Tosatti

A Global Approach to Manipulability Optimisation for a Dual-Arm Manipulator

come

I. INTRODUCTION

Cooperative manipulators have been deeply investigated due to the advantages they guarantee respect to a single-arm system. In particular, their use allows larger payloads and assembly operations that would be infeasible for a single robot [START_REF] Smith | Dual arm manipulation -a survey[END_REF]. Furthermore, multi-arm systems usually present intrinsic redundancy, which may be exploited to achieve secondary goals through optimisation [START_REF] Lee | Relative task prioritization for dual-arm with multiple, conflicting tasks: Derivation and experiments[END_REF]. The drawback is a greater complexity in control and planning, which have to deal with extra degrees of freedom and the choice of appropriate cost functions.

Many approaches have been proposed that define a performance index to be optimised in the redundancy resolution.

For cooperative multi-robot systems, resolution criteria have been proposed that properly redistribute the torque over each manipulator [START_REF] Owen | Acceleration and torque redistribution for a dual-manipulator system[END_REF] or that maximise compliance during the task [START_REF]On-line trajectory resolution for two-armed systems with conflicting performance criteria[END_REF]. Other examples exploit the redundancy to ensure joint-range availability, obstacle avoidance [START_REF] Lewis | Trajectory generation for cooperating robots[END_REF], or to generate human-inspired motion [START_REF] Lee | A cost function inspired by human arms movement for a bimanual robotic machining[END_REF] [START_REF] Lamperti | A redundancy resolution method for an anthropomorphic dual-arm manipulator based on a musculoskeletal criterion[END_REF].

In particular, a great effort has been put in extending the concept of manipulability from single to multi-arm systems. In [START_REF] Lee | Dual redundant arm configuration optimization with taskoriented dual arm manipulability[END_REF], the author approximates the manipulability ellipsoid of a dual-arm robot with the intersection of the singlearm ellipsoids. The approach proposed in [START_REF] Chiacchio | Global task space manipulability ellipsoids for multiple-arm systems[END_REF] extends the definition of manipulability ellipsoid to the multi-arm case for cooperative manipulation tasks by means of a closedchain kinematic analysis. Further contributions have been presented in [START_REF] Bicchi | On the mobility and manipulability of general multiple limb robots[END_REF] and [START_REF] Bicchi | Manipulability of cooperating robots with unactuated joints and closed-chain mechanisms[END_REF], which analyse manipulability in actuated and non-actuated closed-chain mechanism.

A particular case occurs when the assigned task of a dualarm manipulator may be expressed in terms of relative poses between the two end-effectors. As shown in Section II-B, this allows the reduction of the dual-arm system to a single-armlike manipulator. In this way, the manipulability ellipsoids may be directly derived from the so-called relative Jacobian (see [START_REF] Jamisola | Modular relative jacobian for dual-arms and the wrench transformation matrix[END_REF]). This strategy has been applied, for instance, in [START_REF] Lee | Relative task prioritization for dual-arm with multiple, conflicting tasks: Derivation and experiments[END_REF] and [START_REF] Zanchettin | Dual-arm redundancy resolution based on null-space dynamically-scaled posture optimization[END_REF].

A peculiarity of the relative-task approach is the expression of the motion with respect to the end-effectors only. As a consequence, since the system is usually redundant, the motion with respect to the world frame is not univocally defined, and the end-effectors assume different poses depending on the redundancy-resolution criteria.

In any case, the typical approach to redundancy resolution consists in a punctual optimisation, in the sense that the resolution tends to minimise the performance index relative to the actual configuration of the system, as in [START_REF] Zanchettin | Dual-arm redundancy resolution based on null-space dynamically-scaled posture optimization[END_REF] and [START_REF] Chiacchio | Dexterous reconfiguration of a two-arm robot system[END_REF]. This is the case, for instance, of the pseudo-inverse method presented in Section II. The great advantages of these methods are low computational burden and good effectiveness in typical applications. On the other hand, they do not ensure that the mean value of the index is globally optimum. In other words, the solution may be the best for the actual configuration, but the average value of the performance index along the task might not be the optimal one.

A better approach should perform a one-step optimisation to find the configurations that minimise the chosen index throughout the whole task.

This paper, indeed, proposes an algorithm that solves the redundancy issue for a dual-arm system, maximising the manipulability of the whole task, when the desired path between the end-effectors is known a priori. In other words, this method takes into account not only the current configuration of the system, but also the previous and next ones, based on the known path.

A straightforward solution to this issue would be the optimisation of a cost function composed by the sum of the performance indexes relative to every point of the path. Since the resulting computational burden of such approach would be too heavy, an approximated method is proposed.

In particular, the optimisation splits apart into two steps: the first one optimises the sum of the indexes of a certain number of points along the path, whereas the second one optimises the inner points, between the formers.

In this way, the method does not take into account just the punctual configuration of the system, but it tends to maximise the whole task, considering previous and successive configurations as well.

The paper is organised as follows. In Section II, the theoretical background of manipulability and its extension to dual-arm manipulators is presented. In particular, the case of relative tasks and the use of the relative Jacobian is illustrated. Section III discusses the proposed method, illustrating the performance index and the optimisation problem to solve. A simulated implementation of the method is proposed in Section IV, whereas Section V compares the results with the ones given by the typical pseudo-inverse-based method. Finally, conclusions are exposed in Section VI.

II. BACKGROUND

A. Manipulability measure

As mentioned in the previous section, manipulability has been extended to dual arm manipulators in different ways (e.g. [START_REF] Lee | Dual redundant arm configuration optimization with taskoriented dual arm manipulability[END_REF], [START_REF] Chiacchio | Global task space manipulability ellipsoids for multiple-arm systems[END_REF]). Many contributions recall the concept of manipulability ellipsoid, developed for a single arm in [START_REF] Yoshikawa | Manipulability of robotics mechanisms[END_REF].

Briefly, given an n-degree-of-freedom manipulator and an m-dimensional task, the differential kinematics gives:

v = J q (1)
where v ∈ ℝ m is the task velocity vector, q ∈ ℝ n is the joint velocity vector, and J is the m × n Jacobian matrix. Thus, the set of joint velocities with norm equal to one is represented by the unit sphere:

qT q = 1 ( 2 )
which is mapped into an ellipsoid in the task space:

v T (JJ T ) -1 v = 1 ( 3 )
Exploiting the kineto-static duality, the set of joint torques with unitary norm is given by:

τ T τ = 1 ( 4 ) from which: γ T (JJ T )γ = 1 ( 5
)
where τ and γ are the joint and task forces respectively. The former ellipsoid is usually referred to as velocity ellipsoid, whereas the latter as force ellipsoid.

In the case of a 6-dof manipulator, it is easy to notice that the first three rows of the Jacobian refer to the linear velocity transformation, whereas the last three rows refer to the angular one. This means that a manipulability ellipsoid may be decomposed into two further ellipsoids, referring to the linear and the angular velocity respectively. Moreover, given a unit vector u representing a Cartesian direction, the intersection of the force/velocity ellipsoids with u gives respectively:

α = ( u T ( JJ T ) u ) -1/2 (6) β = ( u T ( JJ T ) -1 u ) -1/2 (7)
A meaningful interpretation of ( 6) and ( 7) is the force/velocity transmission ratio of the manipulator. This means that, given a Cartesian task, ( 6) and ( 7) may be exploited to maximise the manipulability of the robot along a desired direction.

Unfortunately, this definition of α and β is not physically consistent, since it involves both linear and angular velocities with different units of measurement (see [START_REF] Schwartz | Non-commensurate systems in robotics[END_REF]). On the other hand, its use as cost function may be considered as a multiobjective optimisation problem composed of the linear and angular terms, weighted by a suitable diagonal matrix W:

ᾱ = ( u T ( JWJ T ) u ) -1/2 (8) β = ( u T ( JWJ T ) -1 u ) -1/2 (9)

B. Relative Jacobian and manipulability

When the desired task of a dual-arm system may be described in terms of relative poses between the two endeffectors, the extension of manipulability ellipsoids comes straightforwardly.

As a matter of fact, in this case a transformation matrix T e 1 e 2 ∈ ℝ 4×4 may be defined, which expresses the pose of the second end-effector frame with respect to the first one. Let T b 1 e 1 , T b 2 e 2 ∈ ℝ 4×4 be the transformation matrices of the first and second end-effectors with respect to their base and T b 1 b 2 ∈ ℝ 4×4 the transformation matrix between the bases (see Figure 1), it results:

T e 1 e 2 = T -1 b 1 e 1 T b 1 b 2 T b 2 e 2 (10)
Moreover, a relative Jacobian J r may be defined that represents a linear transformation between joint and task velocities, such that:

e 1 v r = J r q (11)
where

e 1 v r = [ ẋ ẏ ż ω x ω y ω z ]
T is the velocity of the second end-effector with respect to the first end-effector frame and q = [ qT 1 qT 2 ] T is the vector composed by the joint velocities of the two arms. For the sake of clarity, superscript e 1 will be omitted in the remain of the paper when referring to the relative velocity e 1 v r .

Different ways to derive the relative Jacobian has been illustrated in [START_REF] Lewis | Trajectory generation for cooperating robots[END_REF], [START_REF] Jamisola | Modular relative jacobian for dual-arms and the wrench transformation matrix[END_REF] and [START_REF] Mohri | Cooperative path planning for two manipulators[END_REF]. The present work applies the procedure presented in [START_REF] Jamisola | Modular relative jacobian for dual-arms and the wrench transformation matrix[END_REF].

Once J r has been obtained, the manipulability ellipsoids may be derived as shown previously for the single-arm case.

Moreover, due to the intrinsic redundancy of the system, the inverse kinematic solution becomes analogous to a single redundant manipulator. In particular, deriving the desired relative velocity v d r from a desired task-space trajectory T d e 1 e 2 (t), the inverse kinematic may be written as:

q = J † r v d r + (I -J † r J r ) q0 (12)
where J r is the relative Jacobian of the dual-arm system, J † r = J T r (J r J T r ) -1 is its pseudo-inverse, I is the identity matrix and q0 is an arbitrary vector, which is projected into the null space of J r to satisfy a secondary goal. The projection of q0 into the null space results in self-motions of the joints and, since the task is relative, it generates the motion of the endeffectors with respect to the world frame.

The vector q0 is typically chosen at every step by implementing a gradient-type law, such that:

q0 = -k 0 ( ∂ c(q) ∂ q ) T ( 13
)
where c(q) is a cost function to be minimised and k 0 is a positive scalar. This method implements an effective online redundancy resolution, allowing a local optimisation of the chosen performance index. However, it presents some drawbacks:

• it converges to a local optimum, but it does not ensure the convergence to a global one; • it needs an initial configuration q 0 and the found optimum is usually the nearest to this configuration; • even if the desired trajectory was known a priori, the optimisation would take into account only the actual configuration of the system. Thus, the action q0 might not be optimal for the overall task. Note that the resolution of the kinematics for a relative task does not express any constraints with respect to the base frame. This means the manipulators are allowed to move with respect to the base frame as long as they respond to the relationships of the relative task. As a matter of fact, the joint self-motion makes the end-effectors move with respect to the base frame and the manipulators assume optimal postures depending on the optimisation criterion.

III. PROPOSED METHOD

According to the previous section, it is possible to maximise the manipulability of the system along a direction u by choosing a proper performance index.

Using the relative Jacobian J r in (8) (where W has been omitted for the sake of simplicity), a performance index c(q) may be defined that is inversely proportional to the force transmission ratio along u at each configuration q of the system (see also [START_REF] Chiu | Task compatibility of manipulator postures[END_REF]):

c(q) = 1 α 2 = u T J r J T r u (14)
Minimising c(q), then, corresponds to maximising the force transmission ratio along u.

From a practical point of view, a lower value of c(q) is preferred when the relative task presents high forces on the tool along a given direction u. For instance, a milling process implies elevated forces along the direction tangent to the desired path. Thus, a good value of c(q) in each point of the task ensures lower torques in the joints and energy saving as well.

However, considering a desired path composed by h points and represented by a tensor T d e 1 e 2 ∈ ℝ 4×4×h , a global approach to the optimisation of c(q) over the whole task would consist in finding the h configurations that minimise the total cost ∑ h i=1 c i .

Remark 1:

The choice of a different performance index does not affect the basic concept of the method. In any case, a total cost function ∑ i c i could be defined, in order to optimise the selected performance index c(q) along the overall task.

It is clear that the computational burden of this approach would be too heavy. Thus, let Γ = {γ 1 . . .γ j } be a subset of indices relative to j points on the path (with j < h) and let w(q Γ ) be a cost function, such that:

w(q Γ ) = ∑ i∈Γ c(q i ) (15)
where q Γ = [q T γ 1 , . . . , q T γ j ] T is a column vector composed by the j configurations relative to each element of Γ.

If the chosen points are the most relevant in the definition of the overall task manipulability, then, (15) represents a satisfactory approximation of the complete cost function.

Remark 2:

The choice of the relevant points depends on how the path curvature and the selected direction u change along the task. In other words, in each section of the path, the larger the variation of curvature and u, the greater the number of points to consider.

The minimisation of (15) represents a nonlinear optimisation problem, which leads to the search of j configurations q i∈Γ in the joint space. As a consequence, the optimisation problem is constrained, since:

• the solutions must lay between the mechanical boundaries of the joints:

lb ≤ q i ≤ ub ∀i ∈ Γ (16)
• the solutions must obviously satisfy the nonlinear constraints given by the geometrical path:

T e 1 e 2 (q i ) = T d e 1 e 2 (i) ∀i ∈ Γ (17)

A. Smoothness in joint motion

In order to avoid abrupt changes in the configurations obtained from the optimisation resolution, it is necessary to limit the excursion between the configurations q i∈Γ . For this reason, consider a term w ′ (q Γ ), such that:

w ′ (q Γ ) = ∑ i∈Γ∖{γ 1 } ∥q i -q i-1 ∥ 4 ( 18
)
Note that w ′ (q Γ ) depends on the distance between the j configurations derived from the optimisation and may be utilised to ensure that the optimisation solutions are not abruptly different.

Considering ( 15) and ( 18), the cost function becomes:

w(q Γ ) = ∑ i∈Γ u T i J r (q i )J T r (q i )u i + λ ∑ i∈Γ∖{γ 1 } ∥q i -q i-1 ∥ 4 ( 19
)
where λ ∈ ℝ ≥0 is a weighting factor for the smoothness term.

It is important to mention that the smoothness term in the cost function avoids large differences between consecutive configurations, but does not directly take into account the velocity and acceleration constraints. On the other hand, the proposed method aims at obtaining a series of configurations, each one relative to a geometrical point on the path, and without defining the motion over time. These boundaries have to be taken into account when defining the motion law.

B. Interpolation between the obtained configurations

The optimisation problem presented in [START_REF] Corke | Robotics, Vision & Control: Fundamental Algorithms in Matlab[END_REF] returns the optimal configurations q Γ , but does not calculate the configurations relative to the remaining ones. In order to solve this issue, it is necessary to find a criteria to interpolate the given q i∈Γ . The chosen method is the iterative optimisation of the cost function [START_REF] Corke | Robotics, Vision & Control: Fundamental Algorithms in Matlab[END_REF] for a certain number k of inner points between each couple of q i∈Γ . However, in this case, the cost function must take into account the distance from the initial and final configurations of the interpolation.

In other words, consider each couple of consecutive configurations qγ , qγ +1 , with γ ∈ Γ and let Δ = {δ 1 . . .δ k } be the set of k internal points between qγ and qγ +1 . The cost function to optimise becomes:

ŵ(q Δ ) = ∑ i∈Δ u T i J r (q i )J T r (q i )u i + λ ∑ i∈Δ∖{δ 1 } ∥q i -q i-1 ∥ 4 + λ ( ∥q δ 1 -qγ ∥ 4 + ∥ qγ +1 -q δ k ∥ 4 ) (20)
In this way the algorithm tends to maximise the manipulability between each configuration. As mentioned, this optimisation has to be performed for all the intervals between each couple of consecutive configuration q i∈Γ . This method may be applied recursively to interpolate between the new k+1 intervals until the total amount of points reaches the points of the path (i.e. h).

IV. ILLUSTRATIVE EXAMPLE

A. System setup

The method has been simulated on a dual-arm system composed by two 6-dof Comau NS16 manipulators. The robots lay one opposite to the other, at a distance of 2.5 m, as shown in Figure 2, where the Matlab Robotics Toolbox [START_REF] Corke | Robotics, Vision & Control: Fundamental Algorithms in Matlab[END_REF] has been used as simulation environment.

Thus, the transformation matrix between the bases is:

T b 1 b 2 = ⎡ ⎢ ⎢ ⎣ -1 0 0 2.5 0 -1 0 0 0 0 1 0 0 0 0 1 ⎤ ⎥ ⎥ ⎦
Considering the first robot base as the world frame, [START_REF] Bicchi | On the mobility and manipulability of general multiple limb robots[END_REF] is still equal to:

T e 1 e 2 = T -1 b 1 e 1 T b 1 b 2 T b 2 e 2 (21)
The mechanical bounds of the joints are the same for both arms and are reported in Table I.

B. Desired path

The desired path is defined with respect to the first endeffector frame and consists of a corner with side equal to 150 mm in the horizontal plane z = 0.3 m, with a fillet of radius r = 10 mm, as shown in Figure 3. Fillets and changes of direction represent typical issues in technological application, and this motivates the choice of this path.

Therefore, considering the Cartesian path as composed by h = h 1 + h 2 + h 3 points, where h 1 , h 2 , h 3 are the number

Arm

Bound q 1 q 2 q 3 q 4 q 5 q 6

1 Lower of points of the vertical, rounded and horizontal segments respectively, the relative position x d e 1 e 2 = [x i y i z i ] T is given by:

x i = ⎧   ⎨   ⎩ l ∀i = 1 . . .h 1 (l -r) + r cos ( π 2 i-h 1 h 2 +1 ) ∀i = h 1 +1. . . h 1 +h 2 ( l -r ) ( 1 -i-(h 1 +h 2 +1) h 3 -1 ) ∀i = h 1 +h 2 +1. . . h y i = ⎧  ⎨  ⎩ (l -r) i-1 h 1 -1 ∀i = 1 . . .h 1 (l -r) + r sin ( π 2 i-h 1 h 2 +1 ) ∀i = h 1 +1. . . h 1 +h 2 l ∀i = h 1 +h 2 +1. . . h z i = 0.3 ∀i = 1 . . .h
The relative orientation between the end-effectors is constant and is given by the following rotation matrix:

R d e 1 e 2 = ⎡ ⎣ 1 0 0 0 -1 0 0 0 -1 ⎤ ⎦
Thereby, the desired path may be written as:

T d e 1 e 2 = [ R d e 1 e 2 x d e 1 e 2 0 1 ] ∀i = 1 . . . h
Moreover, the unit vector tangent to each point of the path is:

t i = ⎧  ⎨  ⎩ [0 1 0] T ∀i = 1 . . . h 1 [-sin ( π 2 i-h 1 h 2 +1 ) cos ( π 2 i-h 1 h 2 +1 ) 0] T ∀i = h 1 +1. . . h 1 +h 2 [-1 0 0] T ∀i = h 1 +h 2 +1. . . h
In the simulation, h 1 , h 2 , h 3 were set to 13, 5, 13 respectively, so that h = 31 results.

C. Algorithm setup

First of all, the algorithm has to optimise the cost function [START_REF] Corke | Robotics, Vision & Control: Fundamental Algorithms in Matlab[END_REF], which considers the manipulability indexes of a subset of points Γ. As mentioned previously, the manipulability ellipsoid may be decomposed in a linear and an angular one. In this simulation, the cost function takes into account only the linear ellipsoid, which means the weighting matrix introduced in (8) has to be W = diag(1, 1, 1, 0, 0, 0).

Moreover, the manipulability index has to be optimised along the direction tangent to the path. Thereby, for all points, it must be u i = [t T i , 0] T . Following the guidelines given in Remark 2, the number of elements of Γ is chosen equal to six, which is a good trade-off between computational burden and effectiveness. These points are shown in Figure 3 and correspond to Γ = {1; 7; 13; 19; 25; 31}.

For this first step, the weight for joint smoothness λ was selected equal to 1.5 in order to ensure a good trade-off between the smoothness and the manipulability terms. In the second step, the algorithm minimises the cost function [START_REF] Corke | Robotics, Vision & Control: Fundamental Algorithms in Matlab[END_REF] in order to optimise the manipulability index of the five inner points between each of the previous configurations. In this case, λ was set to 3.

The optimisation has been performed by means of a multiple start gradient-based algorithm, starting from 100 random initial configurations for each problem (for details see [START_REF] Messac | Optimization in Practice with MATLAB ®[END_REF]).

V. RESULTS

The purpose of this section is the comparison between the results obtained with the proposed method and the typical approach presented in Section II. In order to do this, the proposed algorithm is executed as described in the previous section. The optimisation determines h configurations, relative to each point of the path.

In the meanwhile, the pseudo-inverse resolution presented in Section II is performed with k 0 = 1 and using an initial configuration q 0 equal to the first configuration obtained by the proposed algorithm. The cost function c(q) is the same as that defined in [START_REF] Chiacchio | Dexterous reconfiguration of a two-arm robot system[END_REF].

Figure 4 shows the value of the manipulability index (14) along the path, obtained by means of the proposed method. Similarly, Figure 4 shows the same index obtained through the pseudo-inverse-based method.

Comparing the figures, the proposed method shows a flatter trend of the performance index, although the preferred direction u changes through the path. On the contrary, the pseudo-inverse method shows an abrupt increment of the performance index, which corresponds to the change of direction due to the fillet. As a matter of fact, the average performance indexes are equal to 1.17 and 1.52 for the first and second method respectively. This shows that the proposed method ensures a better global manipulability, since it takes into account the whole task during the optimisation.

Moreover, the pseudo-inverse method optimises the cost function from a given initial configuration. This means it does not give any information about the optimal initial posture of the manipulators for a given task. The first configuration found by the proposed method, instead, represents an optimal solution in this sense. This means the task will be performed starting from an optimal initial posture of the manipulators. Thus, a critical issue for the pseudo-inverse Comparison of the pseudo-inverse method for different initial configurations (Blue: solution using the first configuration of the proposed method as initial configuration. Red, green, magenta: solution using random initial configurations; note that the initial configurations must satisfy the constraints ( 16)-( 17)). method is the choice of a suitable initial configuration of the system. As a matter of fact, the use of a non-optimal q 0 leads to a remarkable decay in the trend of the performance index. This is depicted in Figure 5, which shows the trend of the manipulability index, using different initial configurations. In particular, it is easy to note that the trend worsens with respect to the case in which the first configuration of the proposed method is used as initial one (blue line in Figure 5).

Finally, the proposed method allows an easy implementation of joint and Cartesian constraints, which are not easily handled by the pseudo-inverse method.

VI. CONCLUSIONS

A method for the global optimisation of the task manipulability of a dual-arm system has been presented. This approach is based on the optimisation of a cost function, which depends on several points along the desired path.

Since the task is expressed in terms of relative poses between the end-effectors, the dual-arm kinematics is treated by using the relative Jacobian, which makes the resolution of the redundancy similar to a single-arm manipulator.

As shown by simulation, the proposed method allows a better overall manipulability along the task. Moreover, this approach is able to find an optimal initial configuration for the system, which represents a great advantage respect to other methods, like the pseudo-inverse-based one.

As a drawback, the optimisation problem presents a heavy computational burden, and needs to be performed offline.

Fig. 1 .

 1 Fig. 1. Kinematic chain of the dual-arm system.

Fig. 2 .

 2 Fig. 2. Simulation of the dual-arm manipulator setup by means of the Matlab Robotics Toolbox.

Fig. 3 .

 3 Fig. 3. Relative path to be performed by the manipulators.

Fig. 4 .

 4 Fig.4. Comparison of the proposed method and the pseudo-inverse method using the same initial configuration. (Blue: cost function trend using the proposed method. Red: cost function trend using the pseudo-inverse-based method).

TABLE I MECHANICAL

 I BOUNDS OF THE JOINTS FOR THE TWO ARMS.