Henrique Nascimento
email: henrique.nascimento@tecnico.ulisboa.pt

Paulo Rogério Pereira

V-GRADIENT: A Density-Aware Geocast Routing Protocol for Vehicular Delay-Tolerant Networks

Keywords: wireless communications, geocast routing, vehicular delay-tolerant networks

Vehicular Delay-Tolerant Networks (VDTNs) are networks of vehicles that communicate wirelessly, where there are no permanent end-to-end connections. VDTNs have a highly variable topology, with frequent partitions, and possibly low node density. Thus, delay-tolerant routing adopts a store-carryand-forward message transfer paradigm, where messages have a useful Time To Live (TTL) and are stored until a good contact opportunity arises. Multiple message replicas can be generated to improve delivery probability at the cost of increasing network congestion. In this paper, we propose the V-GRADIENT geocast routing protocol that monitors node density and buffer occupancy, to adapt dynamically the forwarding techniques used to disseminate messages within the geographic region of interest. Simulation results show that V-GRADIENT is capable of controlling network congestion and efficiently deliver messages resulting in better delivery ratios (13-99%) and lower latencies when compared with existing protocols.

Introduction

Delay Tolerant Networks (DTNs) are characterized by long or variable delays, intermittent connectivity, asymmetric data rates and high error rates. In this type of networks, routing strategies rely on the store-carry-and-forward paradigm, i.e, on the dissemination of messages to intermediate nodes that can retain them, transport them and deliver them either to the final destination or to other intermediate nodes. The complexity of the routing decisions is, for that reason, closely related with the selection of a proper message replication mechanism.

The geocast concept refers to a delivery scheme that filters the eligibility of the network's target nodes using their location as a criterion. With this in mind, the diversity of geocast environments encompassed by DTNs suggests the demarcation of the strategies according to generic properties such as the definition of regions of interest or the replication mechanisms applied.

Moreover, the coupling between the advances in wireless communications technologies and the car industry makes Vehicle-to-Vehicle connectivity an enabler for a broad range of applications that can be grouped into two categories: safety-related and commercial-related. Therefore, in light of the advances mentioned, this communication paradigm has received a lot of research attention and has motivated the design of routing strategies in Vehicular Delay Tolerant Networks (VDTNs) that can properly handle the following challenges of such distinguishable environment.

In this sense, the design of strategies in VDTNs that can efficiently deliver messages to a group of recipients that would particularly benefit from the reception of such messages has become a widespread topic of discussion. To that extent, the implementation of algorithms that exploit information from the vehicles' navigation system to route messages within specific Regions of Interest (ROIs) motivates the development of the V-GRADIENT protocol.

Contribution to Industrial and Service Systems

The proposed V-GRADIENT geocast routing protocol is the basis for a message dissemination service within a ROI in vehicular highly dynamic or sparse networks. Applications include Vehicle-to-Vehicle (V2V) commercial of safety-related message dissemination. Examples are disseminating fuel prices, current weather or incident warnings. The protocol provides the following contributions: mechanisms to monitor node density and node's buffer occupancy; mechanisms to control message replication avoiding network congestion and buffer exhaustion; efficient message delivery within the ROI.

State of the Art

Unicast Routing

As described in [START_REF] Spyropoulos | Single-copy routing in intermittently connected mobile networks[END_REF], the Direct Delivery protocol is the simplest way to conceive a unicast delivery scheme. The node that carries the message forwards it if and only if it establishes a direct contact opportunity with the target node.

The Spray and Wait protocol [START_REF] Spyropoulos | Spray and wait: An efficient routing scheme for intermittently connected mobile networks[END_REF], in its turn, is based on the dissemination of a limited set of replicas of a given message to distinct contacted nodes during a first phase called spray phase. Afterwards, during the wait phase, relay nodes assume the responsibility of delivering the message to the final recipient directly. However, with regard to the spray phase, several heuristics can be envisioned to distribute properly the messages.

Finally, the Epidemic [START_REF] Vahdat | Epidemic Routing for Partially-Connected Ad Hoc Networks[END_REF] protocol is classified as a flooding-based scheme, since the idea behind it is based on the provision of unlimited message replication throughout the network by exchanging all possible messages whenever contact opportunities arise.

Geocast Routing

On the one hand, protocols such as the GSAF [START_REF] Rajaei | GSAF: Efficient and flexible geocasting for opportunistic networks[END_REF] and the Geoopp [START_REF] Lu | Geoopp: Geocasting for opportunistic networks[END_REF] focus on distributing messages by all nodes located inside a region set remotely, relatively to the position of the nodes that generate such messages. This implies, however, that these algorithms break down the routing process into the two following phases: (1) forwarding (and carrying) the messages to the destination region; (2) delivering the message to all nodes inside the ROI. While, regarding the latter phase, both protocols mentioned follow an epidemic-based behaviour, their strategy differs with respect to the first phase. For instance, when using the GSAF protocol a spraying procedure is applied in a first phase. Messages are assigned a predetermined number of tickets T that denote the number of times a given message can be forwarded to encountered nodes outside the destination region. Thus, each time a message is replicated, the T field is sequentially decreased on both messages until T equals zero.

On the other hand, geocast protocols can be applied to distribute content in ROIs whose center is defined by the position of the nodes responsible for creating the respective messages at the generation instant. In this sense, the authors of the Floating Content algorithm [START_REF] Ott | Floating content: Information sharing in urban areas[END_REF] include in the messages' headers not only the center coordinates and radius of the circularly shaped ROI but also an additional radius parameter to determine the extension of a buffer zone. The buffer zone is defined as an extension of the ROI where carrier nodes may keep the messages storage beyond the boundaries of the ROI.

Geocast Strategies Comparison

This work focuses on distributing messages, during a specified Time To Live (TTL), by all nodes that come across an invariant circularly shaped ROI defined around the node which creates the messages. Moreover, the development of the intended message dissemination algorithm is performed by grasping which features of the considered schemes yield better results under a broad variety of simulations.

In this study, in order to drive such conclusions, the following routing strategies are taken into account, which result from the extension of the definitions regarding the aforementioned unicast routing protocols to a geocast context: Multiple Copy GeoDirectDelivery (MC-GeoDirectDelivery), GeoSpray-and-Wait (GeoS&W), GeoEpidemic and GeoEpidemic Constrained.

Firstly, when using the MC-GeoDirectDelivery protocol, on the one hand, messages are delivered to and kept by nodes residing inside the respective ROI at the delivery instant. On the other hand, messages are retained when nodes leave the ROI, as they may return to the ROI in the future. Because messages are only dropped when their TTL expires or the nodes' buffers are filled up, all nodes that receive a message also become responsible to deliver it to the remaining recipients, enabling the presence of multiple copies of the same message in the network.

Secondly, the GeoS&W scheme couples the principles of the Spray and Wait and the GSAF protocols. In this sense, during the spray phase, messages are replicated, even when nodes are located outside the ROI, using the spraying heuristic adopted by the GSAF protocol. Regarding the wait phase, the GeoS&W behaves as the MC-GeoDirectDelivery protocol.

The Epidemic Protocol nature may assist message distribution in a geocast context as it already distributes an unbounded number of copies of the same messages by multiple receivers. However, to distinguish the protocol's operation in unicast scenarios from its operation in a geocast routing paradigm, this protocol will be referred as the GeoEpidemic protocol throughout the article. Likewise, the GeoEpidemic Constrained scheme follows the Epidemic protocol's principles as well. However, in this case, the routing scheme under consideration constrains the dissemination of messages within the boundaries of the ROI, since messages are dropped as soon as the nodes carrying them leave this region.

Additionally, in this line of reasoning, it is also crucial to make explicit the following performance assessment metrics utilized: delivery ratio and delivery latency. The delivery ratio is computed using the formulation developed in [START_REF] Rajaei | GSAF: Efficient and flexible geocasting for opportunistic networks[END_REF]. By monitoring the recipients that reside inside a given message's ROI during its TTL, it is possible to compute a ratio between the number of successful receivers and total number of nodes belonging to the list of recipients. Obviously, this ratio only translates itself into a per message result. For that matter, to obtain an overall delivery ratio estimation, the per message metric is extended to all created messages through an average operation. Latency is also monitored, in a first instance, taking into account a single message. However, in this case, the per-message metric is computed by summing the amount of time it takes for the message under consideration to reach each one of its successful receivers and, then, by dividing this outcome by the same number of recipients.

V-GRADIENT Design

The major points of improvement attained by benchmarking the different geocast strategies relate not only with a moderated replication process applied outside the messages' regions of interest but also with a dropping policy that can dynamically estimate if it is beneficial to keep messages being carried by nodes located outside the ROI. With this in mind, on the one hand, with regard to the replication process, the V-GRADIENT algorithm incorporates the principles of the Spray-and-Wait protocol [START_REF] Spyropoulos | Spray and wait: An efficient routing scheme for intermittently connected mobile networks[END_REF] in a geocast context and, during the spray phase, the methodology adopted by the GSAF protocol [START_REF] Rajaei | GSAF: Efficient and flexible geocasting for opportunistic networks[END_REF]. On the other hand, the deployed strategy follows the idea behind the Floating Content Algorithm [START_REF] Ott | Floating content: Information sharing in urban areas[END_REF], i.e., defines a buffer zone whose range varies according to the network conditions.

For such purposes, to select dynamically the buffer zone range, the V-GRADIENT algorithm takes advantage of estimations of the network density level and of the buffer occupation level, updated periodically. Because in this study it is assumed that nodes do not have access to any knowledge oracles, vehicles rely only on information exchanged during contact opportunities to perceive the surrounding network conditions and, for that reason, both awareness metrics are kept as dictionaries in the nodes' buffers in a decentralized manner. Moreover, the mathematical formulation used to compute both metrics follow an Exponentially Weighted Moving Average (EWMA). For the density level metric ρt, the EWMA is performed every 30 seconds taking into account the number of distinct contacts nc established, as shown in Eq. 1. α was set to 0.25.

ρt = α•nc + (1-α)•ρt-1 . (1)
The buffer occupation metric βt is determined by averaging the fractions of the buffer occupancy βlevel, varying from zero to one, retrieved from neighbouring nodes in 30 seconds intervals with an EWMA, as shown in Eq. 2, where γ was set to 0.7.

𝛽 𝑡 = 𝛾 ∑ 𝛽 𝑙𝑒𝑣𝑒𝑙 𝑛 𝑐 𝑛 𝑐 + (1 -𝛾)𝛽 𝑡-1 .
(

) 2
In order to carry out a reliable implementation of the mechanism responsible for adjusting the extent of the buffer zone, a wide variety of simulations were executed using the ONE simulator [START_REF] Keränen | The ONE Simulator for DTN Protocol Evaluation[END_REF] under different density scenarios, considering ROIs with different radius values, and in environments that led to serious buffer congestion issues, as shown in Table 1. Accordingly, in each simulation an additional dropping policy was included that periodically removed messages from storage if the nodes were located outside of the buffer zone range, defined as the product between the ROI radius, RROI, value and a threshold R+ratio. Vector of 100 evenly spaced points between 1 and 10 Furthermore, to exhibit the impact of the network density conditions and the ROIs' radius values on the additional dropping policy formulation, Table 2 includes the optimal R+ratio values, bearing in mind the delivery ratio results obtained. It can also be seen in the same table that, for a fixed number of nodes deployed in the network, the optimal R+ratio expression can be modeled by an equation obtained through a power regression. The profile of such equations, since it also follows a trend depending on the number of nodes and consequently on the density level, also suggests a mathematical formulation of a variable threshold value associated with the extension of the buffer zone as shown in Eq. 3.

𝑅 +𝑟𝑎𝑡𝑖𝑜 (𝑅 𝑅𝑂𝐼 , 𝜌 𝑡 ) = 𝑔(𝜌 𝑡 ) • 𝑅 𝑅𝑂𝐼 ℎ(𝜌 𝑡 ) . (3)
With respect to the formulation of the 𝑔(𝜌 𝑡 ) and the ℎ(𝜌 𝑡 ) expressions, by selecting their corresponding target parameters presented in Table 2, it can be observed, as depicted in Fig. 1, that both describe a well modeled decay relatively to the average density level ρavg measured in simulations deploying a different number of nodes. In this sense, in scenarios characterized by extreme buffer congestion issues, the V-GRADIENT algorithm allows a dynamic adjustment of the range of the buffer zone for all messages carried by a particular node according to the estimated threshold expression presented in Eq. 4, as a result from the best fit from Fig. 1. Additionally, the buffer occupation level metric plays an important role on detecting the occurrence of congestion, providing nimbleness to the dropping policy operation. If we let ourselves picture a scenario where only a few messages are generated per hour, it can be concluded that, because no congestion issues arise, the drop policy should be less restrictive. For this reason, a logistic function f(βt) responsible for regulating the stringency of the policy, with βt = 0.5 defining the sigmoid's midpoint, is included in the final definition of R+ratio.

𝑅 +𝑟𝑎𝑡𝑖𝑜 (𝑅 𝑅𝑂𝐼 , 𝜌 𝑡 , 𝛽 𝑡 ) = 1 𝑓(𝛽 𝑡 ) 10413 • 𝜌 𝑡 -0.243 • 𝑅 𝑅𝑂𝐼 1.3731-0.016•𝜌 𝑡 . ( 4
)

Evaluation and Results

In this section, the results regarding the different strategies described in Section 3 are compared against the ones concerning the V-GRADIENT algorithm, considering a simulation scenario following most of the parameters' configuration presented in Table 1. However, in this case, different message generators are deployed in each simulation simultaneously. Additionally, the nodes' buffer size is doubled to 10 MB and four different seed values are used for the movement models' pseudo random number generator to reinforce the statistical significance of the results obtained. The corresponding 95% confidence intervals are presented in the following graphs. Moreover, to assess the performance of the V-GRADIENT strategy under different density conditions, results are presented in parallel according to the number of nodes deployed in the network. Regarding the delivery ratio assessment, as depicted in Fig. 2, the V-GRADIENT algorithm yields better results than the remaining protocols. In fact, for the simulations executed with 150 nodes deployed, the V-GRADIENT algorithm shows a 13% increase on the delivery rate relatively to the GeoS&W protocol, 15% relative to the MC-GeoDirectDelivery, 39% relative to the GeoEpidemic Constrained and 99% relative to the GeoEpidemic. Furthermore, the remaining protocols' outcome also reinforce the choices related with the V-GRADIENT design. For instance, the pronounced upward trend shown by the GeoEpidemic Constrained protocol's results, as the density level increases, validates the exercise of a dropping policy that renders itself stricter in high density conditions. When analyzing latency results, one should also bear in mind the delivery ratio results to draw respectable conclusions, since only the delivered messages are taken into account when it comes to measuring latency. Accordingly, in spite of the GeoEpidemic Constrained achieving better results, it can be assumed that, because messages are deleted from the nodes' buffers as soon the carrier nodes leave the ROI, they only live a fraction of their total assigned TTL, which compromises the delivery ratio efficiency. In this sense, the V-GRADIENT provides a moderate and adequate solution, putting a limit on latency, by means of a dynamic dropping policy, while still achieving good performance regarding the delivery ratio: 20% better than GeoS&W.

Conclusions and Future Work

We proposed the new V-GRADIENT geocast routing protocol for VDTNs that monitors node density and buffer occupation to disseminate efficiently messages within a ROI centered in the source node. The protocol achieves a better delivery ratio (13-99% with 150 nodes) and lower latency than the existing protocols by controlling network congestion. As future work, we intend to extend the operation of the V-GRADIENT protocol to an environment encompassing distinct geocast groups associated with several types of messages' applications, taking into account the penetration rate of such groups in the network. We also intend to test other scenarios and evaluate the protocol overhead.

Fig. 1 .

 1 Fig. 1. Regression-based model for the g(ρt) and the h(ρt) expressions.

Fig. 2 .

 2 Fig. 2. Delivery ratio and delivery latency as a function of node density and protocol.

Table 1 .

 1 Assigned parameters during the simulation procedure.

	Map	Downtown part of the city of Helsinki [4500m x 3400m]
	Simulation time	12 hours
	Movement model	Shortest Path Map-Based Movement
	Type of nodes	Vehicles (cars)
	Nodes' speed interval	10-50 km/h
	Nodes' buffer size	5 MB
	Nodes' wait time	5-30 seconds
	Message size interval	500 kB-1 MB
	Message generation interval	30-90 seconds
	Interfaces' data rate	250 kbps
	Interfaces' transmission range	50 meters
	Initial Number of Copy Tickets 3
	Message's TTL	150 minutes
	Number of nodes	[50;100;200;400] nodes
	ROI Radius	[250;500;750;1000] meters
	R+ratio

Table 2 .

 2 Optimal R+ratio values obtained for simulations considering different density conditions and distinct radius values defining the ROI.

Acknowledgments. This work was supported by national funds through Fundação para a Ciência e a Tecnologia (FCT) with reference UID/CEC/50021/2019.


Tiago A Antunes
email: tiago.antunes@siemens.com

Rui Castro
email: rcastro@tecnico.ulisboa.pt

P J Santos
email: paulo.santos@estsetubal.ips.pt

A J Pires

Matthias Foehr
email: matthias.foehr@siemens.com

Optimized Electrification of Subsea Oil & Gas Infrastructures Based in Genetic Algorithm

Keywords: subsea, substations, oil & gas, optimization, genetic algorithm

Offshore field development relies on multiple optimization techniques targeting a feasible and cost-effective production solution yet are focused on the field itself. While so, advancements in offshore engineering bring increasingly complex subsea infrastructures to depths in the excess of 3,500 meters. Many offshore production topsides which currently rely on costly and harmful onboard thermal-based power generation are turning to high voltage power-from-shore electrification solutions to cope with the challenges being brought by subsea infrastructures. An optimal electrification of these subsea templates is a challenge on its own as the seafloor morphology and well distribution is far from consistent. This paper presents a combined k-means and genetic-algorithm optimization to assess how the combined deployment of high voltage umbilical, wellheads and subsea substations can be optimized for the lowest cost possible. Results show a significant improvement in optimization of the total umbilical length as well as the substation positioning on the seabed.

Introduction

The International Energy Agency (IEA) stresses that projects in the main areas such as the Gulf of Mexico and the North Sea, were once characterized by a break-even ratio of 60-80USD/barrel. Equivalent projects are made viable and as robust as breakeven oil prices of 25-40 USD/barrel [START_REF]Offshore Energy Outlook[END_REF]. This can only be achieved via proper simplification, modularization and downsizing; Electrification is at the core [START_REF] Rajashekara | Power electronics for subsea systems: Challenges and opportunities[END_REF][3].

Offshore power production is based in thermal combustion systems whether this is accomplished via diesel or gas and that can reach several dozen megawatts [START_REF] Shi | Grid integration of offshore wind farms and offshore oil/gas platforms[END_REF] [START_REF] Hamdan | Using a VSC Based HVDC Application to Energize Offshore Platforms from Onshore -A Life-cycle Economic Appraisal[END_REF]. Modularization and optimization mean that owners can increase the complexity of the templates, reduce the size of the topside and increase the overall efficiency and resilience. Furthermore, unmanned operations, expected to dramatically reduced Operational Expenditure (OpEx), will require a reliable power supply with a large maintenance window, possible by means of power-from-shore (PFS) solutions [START_REF] Yamamoto | Integrated subsea production system: An overview on energy distribution and remote control[END_REF].

In subsea facilities commonly ranging up to 100MW [START_REF] Hazel | Taking power distribution under the sea: Design, manufacture, and assembly of a Subsea electrical distribution system[END_REF] the matter is addressed on a point-to-point or direct in-place analysis of the most suitable electrical infrastructure and not on the optimization of the actual electrical grid on the seafloor. The authors have reviewed field development cases [START_REF] Wang | A new partition model for the optimization of subsea cluster manifolds based on the new definition of layout cost[END_REF] [START_REF] Wang | Modeling for the optimization of layout scenarios of cluster manifolds with pipeline end manifolds[END_REF], yet those references are focused on the subsea production setup and not on the actual electrification.

The authors were able to review several case studies where satisfactory results were obtained by means of GA [START_REF] Taroco | A Faster Genetic Algorithm for Substation Location and Network Design of Power Distribution Systems[END_REF][11] [START_REF] Gen | Hybrid genetic algorithm for multi-time period production/distribution planning[END_REF] -which sustains the suitability of the approach. The same algorithm was also used in other examples, such as in the design of S/S's earthing grids [13], switchgear internal design [START_REF] Hinow | Genetic algorithm based methodology for optimisation of innovative switchgear design[END_REF] or in life-cycle assessment of electrical S/S [15]. References show that simulated annealing can also be applied for this sort of optimization, often times combined with GA [16][17].

The rest of the paper is organized as follows. On the second chapter, a connection is established between this research and the resilient systems topic. The third chapter provides an overview of the main industry challenges, which shall serve, on the forth chapter, to support the model design. The results and conclusions are presented on the fifth and sixth chapter, respectively.

Relationship to Industrial and Service Systems

Although PFS on itself is not an innovation, and field development is a consolidated environment, no references to support the use of Genetic Algorithms (GA) or others for optimized subsea electrification were found. However, the application of GA on the optimization of power systems, more specifically, in the substation (S/S) expansion planning (SEP) is widespread and quite common. Companies like SIEMENS want to drive large-scale deployment of PFS solutions as an alternative to onboard thermal allowing utilities to take full advantage of the massive renewable capacity available.

Having so, and whilst the authors assess new ways to power Oil & Gas (O&G) infrastructures from shore, this paper presents their novelty approach towards bringing a much-needed solution for the cost optimization of the subsea electrical setup. It is designed by means of a Python-based mixed Kernel-Means (KM) and GA, considering simultaneously the layout of the field and the electrification concept to get the most cost-effective solution.

Industry Challenges

• An analysis of over one hundred offshore sites showed electrification costs ranging between 0.3-0.5 MUSD/kboed 1 despite the variations on the resource or water depth. In the case of subsea infrastructures, the ratio was 0.5; • Subsea production facilities are characterized by largely dispersed and smallsized equipment. An analysis of all the possibilities can become timeconsuming and the GA does not require the evaluation of all possibilities [START_REF] Rexhepi | Analysis of the impact of parameters values on the Genetic Algorithm for TSP[END_REF];

• Industry leaders are working towards subsea electrical power standardization, which supports the claim that robust subsea electrical systems are both new and indeed needed at a fast pace [19]; On one hand, unsupervised methods such as KM clustering bring a straightforward method to analyze the positioning of the wells without initial candidate locations [17]. As it does not require initial mapping it serves as a solid method to suggest individual candidates on the GA application [START_REF] Vahedi | Pseudo-dynamic substation expansion planning using hybrid heuristic and genetic algorithm[END_REF]. On the other hand, GA are known for providing high quality and performance results, especially for a high number of possibilities [START_REF] Rothlauf | Representations for Genetic and Evolutionary Algorithms, 2nd Editio[END_REF]. In addition, they show enough flexibility to accommodate several restrictive measures -such as the number of power outlets allowed on a S/S and the maximum cable distance.

Proposed Method

Initial Assumptions

An extensive analysis of the offshore portfolio was made. Wellheads are distributed along a squared area, representing field section. The number of wells is based on current projects. Focus was given to the evaluation of the best possible electrification for any given field development considering the following portfolio restrictions:

• Electrical distribution is made typically using voltages ranging between 3.3 and 36kV [19] [START_REF] Hazel | Taking power distribution under the sea: Design, manufacture, and assembly of a Subsea electrical distribution system[END_REF] and oil-filled, pressure-balanced and coated electrical distribution units (EDU or SDU) are used to power the subsea systems. These units rely on templates pre-assembled on surface [START_REF] Bai | Representations for evolutionary algorithms[END_REF]; • Electrical distribution manifolds (EDM) -a subsea substation, can also be used to provide power to multiple locations at the same field whilst the number of connectors is kept to a minimum to ensure overall system resilience [START_REF] Bai | Representations for evolutionary algorithms[END_REF]; • One power transformer is planned on each S/S and the technology used, high voltage alternated current (HVAC) or direct current (HVDC), is not compared as it depends on the distance to shore; And this optimization serves for both; • The number of power distribution outlets of the EDU is estimated between 3 and 6 at a given voltage level whereas subsea wellheads are assumed to be rated the same [3]; • The number of S/S is estimated based on the quantity of wellheads, as the cost of the PDU is largely higher than the HV umbilical an estimation based solely in the number of connections to do was used [START_REF] Bai | Representations for evolutionary algorithms[END_REF]; • The initial positioning of the S/S is done via the K-means clustering method whereas the GA populates multiple candidate substation positioning. The size of the solution array is constant which means that all subsea templates are powered by one and only one of the subsea S/S; • The solutions (individuals) share traits (S/S connections) between each other and the GA application is able, not only to foster the better combination of those traits, but also, to short-list the better candidate solutions [22].

A joint approach on both hydraulic and electrical power subsea distribution was not made although it may lead to further improvements. The geography of the seabed and fields are outside the focus of this paper. The objective is an universal canvas.

Process Flow

As said, KM is used as an unsupervised learning method for clustering the field data. This algorithm provides each individual with a candidate position and the distance to each wellhead to which respects to. Please refer to Fig. 1. for a process flow example. Subsea wells are identified in blue whereas the S/S are identified in green. The links are mapped accordingly. This is the baseline for the GA, representing one individual from the sense of the S/S location. For any well and subsea S/S set, the objective is to minimize the total cost of the electrification solution (cumulative cost of the subsea cabling for well-substation link). The cost objective function (1) to be minimized is presented next, where L represents each cable distance, C and CKM the total and unitary cost, respectively:

C = LCABLE-N * CKM (1)
The connections are represented via a single-digit code which states the reference of the respective S/S. That allows for a direct and natural encoding. The same node (position) cannot be completed by neither none nor more than one S/S -an example is shown on Fig. 2. Script-wise, each individual is made of an integer array with a length equal to the number of load centres (wells).

Fig. 2. Example of an individual as used on the GA

The distance is calculated as the shortest straight line between the two-points (S/S and well). Such calculation outputs a similarly-structured array. The decision variable is the total cost for an individual (composed of the cumulative cost of each wellhead/substation tuple). Typically, for this analysis, the restrictions are the spatial representation of the wells, S/S and connection links. All individuals are deemed acceptable and the tuples are selected within the number of wells and S/S proposed.

Characterization of the GA Population

The population size was set so that a balance between both the number of S/S and templates to power is reached. Sensibility tests were made to assess the appropriate setup. The best compromise was found a population equal to the number of wellheads (links) with a medium crossover rate and a very low mutation rate. The iteration rate was kept high during most cases, although fitness is achieved quite early. The population size is the most critical variable although the mutation rate must be kept to a minimum (up to 10%) given the large population size. [START_REF] Hazel | Taking power distribution under the sea: Design, manufacture, and assembly of a Subsea electrical distribution system[END_REF].

Ranking, Crossover and Mutation

The total cost per each individual is calculated by means of the sum of the length of all links combined and the unitary cost. Fig. 3 presents an example of the elite selection by addressing the total cost (per individual) array of the population.

Fig. 3. Example of a cost array as used on the GA

A portion of the remaining population is selected for cross-over. An even number of individuals are selected, and a proper parent-children pair can be selected. The cross-over is done using the sexual reproduction principal with two individuals, done via cross-over pointing, where the two new sub-set strings are swapped and will form the children individuals. An example is presented in Fig. 4. The parents are randomly selected, and each tuple provides one children. For that reason, the number of cross-over iterations matches the number of parents (population size remains unchanged). There is no preference on the parents and the two new individuals from each round proceed for the mutation pool.

Parent A 1 2 2 2 3 3 1 2 1 3 3 2 2 1 3 Parent B 1 3 1 1 2 2 3 3 2 2 2 1 3 3 2 Child A 1 2 2 2 3 3 3 3 1 3 3 2 2 1 3 Child B 1 3 1 1 2 2 1 2 2 2 2 1 3 3 2
The mutation is done to the pool of all individuals with an exception for the top fitness one, thus including the untouched individuals, as well as, the cross-over children. The low rate mutation is applied on a single random point of an individual. One link (S/S reference) is changed in that percentage of individuals. In order not only to keep accurate results and avoid a recalculation of the costs, but more importantly, to reanimate possible important solutions, the "mutated" individual receives a gene (link and distance tuple) from another individual of the pool, also randomly selected. The process is shown in Fig. 5.

Original 1 2 2 2 3 3 1 2 1 3 3 2 2 1 3 Mutated 1 2 2 2 3 3 1 3 1 3 3 2 2 1 3 Donor 1 3 1 1 2 2 3 3 2 2 2 1 3 3 2

Fig. 5. Example of a mutation operation as used on the GA

The distance is already available, taken also from the equivalent position on the distances' array. This avoids the need to recalculate all distances of the current population and optimizes computation time. The standard deviation is used to assess the performance of the GA. This process is iterative and is repeated for a given number of times.

Results

Initial Conditions

The test conditions are presented in Table 2. The mapping is reviewed against the subsea grids and O&G fields known. The cost for the umbilical is indicative and the number of S/S is selected based on maximum number of connections of each unit (PDU), in this case 5 wellheads each. The optimized deployment is shown on Fig. 5. The grey line represents the HV umbilical laid between the wellheads and the respective S/S. The authors acknowledge that the solutions proposed are not perfect and additional measures should be taken to avoid dispersion of the S/S along the field during GA execution.

Cost Minimization

The results of the total distance and cost for each iteration are presented in Figure 6.

Fig. 7. {Cost, distance} per iteration

In the same manner, the standard deviation is presented in Figure 7.

Fig. 8. {Standard distance deviation, cost deviation} per iteration

There is a quick convergence towards the best cost position. However, as already commented, the higher rate of iterations was defined so that a sustained convergence of the overall population was possible. Results are summarized on Table 3. In some scenarios, the same cost was found for two simultaneous positions. Yet, these are a minor occurrence on the overall spectrum of scenarios addressed. The routine is, however, able to identify the individuals which match this best fitness performance.

In fact, using a KM before the GA implementation dramatically reduces the computation time for the optimization as the centroids are calculated (via the Euclidean distance) close to their final positions in the majority of cases. Increasing the threshold was used to accommodate the lack of convergence, in multiple cases, of the standard deviation. Efforts were made to adjust the cross-over and mutation ratios to improve such convergence.

An increase on the mutation rate is known to allow for an optimization of the computation of the solution which supports the claim that associative optimization problems can be served with better results by using mutation rates higher than crossover rates. The increased cross-over ration introduces a very high disturbance on the standard deviation, which should be avoided. On the other hand, the increase mutation factor causes the abrupt convergence of the least cost having a similar deviation conversion as of the initial case scenario.

Conclusions

It was concluded that the paper presents an appropriate solution for an optimized deployment of a subsea infrastructure by means of a KM and GA-based routine. A satisfactory optimization of the total cost was achieved, and an appropriate narrowing of the individual's fitness was also accomplished. The authors also acknowledge that there is enough room to innovate as multiple deployment restrictions have to be accounted for in such models.

The authors trust that a successful tool will certainly include additional optimization tiers (still using the GA or other suitable algorithm) so that topics like the substation positioning vs the HV umbilical, avoidance of unsuitable paths for the links or maximum tie-in length topics are addressed.

Fig. 1 .

 1 Fig. 1. Deployment of wellheads, substations and umbilical (algorithm process)

Fig. 4 .

 4 Fig. 4. Example of a cross-over operation as used on the GA

Fig. 6 .

 6 Fig. 6. Field distribution with substations deployed (Case studies)

Table 1 .

 1 Constrains of the Genetic Algorithm

	Rate designation	Variable	Value	Units
	Population (individuals)	ind_rate	1	[-]
	Iterations of the GA	ite_rate	200	[-]
	Top fitness rate	-	1	[-]
	Crossover rate	cross-rate	40	%
	Mutation rate	mut_rate	10	%

Table 2 .

 2 Key data of the Case Study

	Constrain description	Value	Units
	Map resolution	100	km
	Number of substations	3	[-]
	Number of wellheads	15	[-]
	Number of iterations	1,000	[-]
	Unitary umbilical cost	5,000	USD/km

Table 3 .

 3 Final results

	Variable identification	Value	Units
	Total of Iterations	1,000	-
	Initial Length	258	km
	Final Length	235	km
	Initial Cost	1,280	MUSD
	Final Cost	1,170	MUSD

22 35 12 14 18 19 20 12 23 25 17 15 20 22


Bruno S Miranda
email: miranda.bsm@live.com

Akebo Yamakami

Priscila C B Rampazzo
email: priscila.rampazzo@fca.unicamp.br

A New Approach for Crop Rotation Problem in Farming 4.0

Keywords: Crop Rotation Problem, Farming 4.0, Precision Agriculture, Genetic Algorithm, Farm Management

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

Introduction

Technological advances have been widely shaped agriculture. Innovations in sensor devices and embedded systems have improved tillage and crop yields. Efficiency in resources management and autonomous data acquisition turned out to be essential among farmers. As Industry 4.0 thoughts have pushed ahead manufacturing units, agriculture also moves toward a transformation called Farming 4.0, or also referred to as Agriculture 4.0. In general, important developments in automation have been noticed, related to Precision Agriculture expansion [START_REF] Braun | Farming in the Era of Industrie 4.0. 51st CIRP Conference on Manufacturing System[END_REF] [START_REF] Pereira | A review of the meanings and the implications of the Industry 4.0 concept[END_REF].

Precision agriculture (PA) follows from the integration of crop management and information technologies. The field measurement of crop requirements and the proper supply of these demands are promoted by PA and intend to improve crop production and resource consumption, approaching sustainable ideas and reducing environment impacts [START_REF] Pandey | Agricultural Cyber-Physical System: In-Situ Soil Moisture and Salinity Estimation by Dielectric Mixing[END_REF] [START_REF] Far | Impacts of the precision agricultural technologies in Iran: An analysis experts' perception & their determinants[END_REF]. Industry 4.0 relies on digital technologies, such as Big Data, Internet of Things (IoT) and Artificial Intelligence (AI), to enhance production systems and reach the highest level of cooperation and mobility. In agriculture, connectivity has been expanding allied with IoT and transforming the whole infrastructure. This trend has been noticed in farm equipment: connected tractors, spray devices with route optimization devices and other advanced machines, which are able to gather data and even share information with their manufactures for maintenance purposes [START_REF] Bonneau | Industry 4.0 in agriculture: Focus on IoT aspects. European Commission, Digital Transformation Monitor[END_REF].

Besides farm technologies stand as advanced as industrial applications, production planning remains a challenging task in agriculture. Aside from monoculture, scheduling crops in a large set of fields rely on market outcomes, weed, and pest control or weather forecasts. Searching for strict mathematical approaches remains quite unusual, relying much on the farmer expertise and short-term goals.

According to USDA, Latin American agricultural GDP growth ought to recover strength in a short-term period marked by recession and slow growth, which relates mainly to Brazil and Argentina economic instability [6]. Brazil is one of the greatest food producers in the world. Investments in agrarian researches and development are quite expressive, they reached 1.82% of Brazilian GDP in 2013 [7]. Therefore, agrarian researches present as much potential as industrial ones for Latin America economies.

This research paper acknowledges the importance of CRP in the Farming 4.0 developments. The proposed mathematical model for CRP bases on fertilization management techniques, such as nutrient budgets.

The main contribution of this paper relies on the relevance of the proposed model supporting farm management decisions. It achieves profitable solutions and evaluates the resource balance, providing valuable insights in the crop sequence.

To accomplish a broad range of CRP scenarios, optimization techniques ought to provide solutions for small and sizeable instances of the problem. The proposed evolutionary algorithm presents good results in the large instances where the deterministic method takes large amounts of time. The stochastic approach in this paper combined with the deterministic method ensures that the proposed model can assist agrobusiness management providing solutions for many scenarios.

The subsection 1.1 presents a review about the CRP approaches in the literature. Section 2 establishes the relationship to industrial and service systems. Section 3 presents the proposed model and its details. The optimization techniques are discussed in Section 4 and the attained results in Section 5. The conclusion and further steps of this research presents in Section 6.

Discussion on Previous Works

The CRP has been extensively researched. Distinguished approaches and mathematical models have enriched widely the available literature. The complexity related to this NP-hard problem is the main reason for continuous innovation. Heuristic and metaheuristic approaches are quite common among applied techniques.

A review of the stochastic and deterministic methods applied to the CRP presented in [START_REF] Memmah | Metaheuristics for agricultural land use optimization. A review[END_REF] and pointed out that the deterministic approach has been successful in some specific models, though lacks in efficient over sizeable instances of the problem. Although there is no agreement among researches about the best method, evolutionary and hybrid algorithms are promising, due to the high flexible structure of genetic algorithms.

Concerned about a large amount of consumed resources and the environmental impact of monoculture farming, the research presented in [START_REF] Santos | Crop rotation scheduling with adjacency constraints[END_REF] discussed an optimization model for CRP, based on organic farming concepts. A new column generation algorithm combined with a greedy heuristic was developed and generated solution for this proposed model of the CRP.

The research presented in [START_REF] Aliano | Metaheuristics for a crop rotation problem[END_REF] analyzed hybrid metaheuristic algorithms in search of quality solutions for the CRP. Attaining feasibility presented to be a hard task on the proposed model, to overcome this hurdle, the initial population was generated by a heuristic procedure. The hybrid algorithms with local search and with Simulating Annealing presented good results in this related work.

Multi-objective approaches for the CRP are quite contemporary. A bi-objective approach presents in [START_REF] Aliano | Metodologias de escalarizações para o problema de rotação de culturas biobjetivo[END_REF], it acknowledges that profitability and diversity of crop rotation are conflicted goals and provides a bi-objective model that explores both objectives.

Relationship to Industrial and Service Systems

Based on a vast area of expertise, Operational Research deals with problems regardless of the context in which they arise. Problems in which the objective is to determine, according to one criterion, the best choice within a set of alternatives. The development of the areas of engineering, computing, and economics has been characterized by the increasing use of optimization models as paradigms for representation and resolution of decision-making problems.

Scheduling Problems represent a class of significant decision-making problems in optimization. Numerous companies and service systems face problems regarding task sequencing, which can be caused by improper allocation of resources and poorly defined processes. The areas affected by these problems are diverse: industry, manufacturing, agriculture, process management, transportation, among others. We can optimize processes by performing task sequencing planning, which results in improved production flow control, meeting deadlines, and scheduling tasks to better utilize available resources.

Solutions of this type of problem are not obtained in closed forms. Instead, they are determined by algorithms: a sequence of procedures applied repeatedly to the problem until the best solution is obtained. According to the adopted formulation (fundamentally dependent on the parameters), we can solve the problem exactly in polynomial time, or even deal with NP-class problems. In these cases, it is necessary to abandon the search for an optimal solution to seek a quality solution, through heuristic procedures.

The present research aims to provide a valuable method to analyze information only acquired by the current state of interconnection between physical and cyber world, such as actual crop nutrient requirements. And, although the proposed methodology applies to agrarian management, the concepts of Industry 4.0 are embraced in this research. The perception of this work exceeds agricultural environment by the integration of planning, profit and production requirements.

Crop Rotation Problem

The proposed mathematical model for the CRP is described ahead and bases on mathematical approaches presented in [START_REF] Santos | Crop rotation scheduling with adjacency constraints[END_REF] [START_REF] Aliano | Metaheuristics for a crop rotation problem[END_REF]. The nutrient balance concept presented in the model relates to researches in [START_REF] Watson | A review of farm-scale nutrient budgets for organic farms as a tool for management of soil fertility[END_REF] [13] and bases on the surface nutrient budget idea.

• 𝑁: size of the crop set;

𝑚𝑎𝑥 ∑ ∑ ∑ 𝑎𝑟𝑒𝑎 𝑘 • 𝑙 𝑖 • 𝑝 𝑖 • 𝑥 𝑖𝑗𝑘 𝐿 𝑘=1 𝑀 𝑗=1 𝑁 𝑖=1 -∑ ∑ 𝐹 𝑁 𝛼𝑘 • 𝑐 𝑁 + 𝐹 𝑃 𝛼𝑘 • 𝑐 𝑃 + 𝐹 𝐾 𝛼𝑘 • 𝑐 𝐾 𝐿 𝑘=1 𝛼∈Ω (1) 𝑠. 𝑎 ∑ ∑ ∑ 𝑥 𝑖(𝑗-𝑟)𝑣 𝑣∈𝑆 𝑘 𝑡 𝑖 -1 𝑟=0 𝑖∈𝐹 𝑝 ≤ 𝐿 • (1 -∑ ∑ 𝑥 𝑖(𝑗-𝑟)𝑣 𝑡 𝑖 -1 𝑟=0 𝑖∈𝐹 𝑝 ) , 𝑝 = 1, … , 𝑁 𝑓 , 𝑗 = 1, … , 𝑀, 𝑘 = 1, … , 𝐿 (2)
∑ ∑ ∑ 𝑥 𝑖(𝑗-𝑟)𝑘 𝑣∈𝑆 𝑘 𝑡 𝑖 +𝛽 𝑟=0 𝑖∈𝐹 𝑝 ≤ 1, 𝑝 = 1, … , 𝑁 𝑓 , 𝑗 = 1, … , 𝑀, 𝑘 = 1, … , 𝐿 (3)
∑ ∑ 𝑥 𝑖(𝑗-𝑟)𝑘 𝑡 𝑖 -1 𝑟=0 𝑁 𝑖=1 ≤ 1, 𝑗 = 1, … , 𝑀, 𝑘 = 1, … , 𝐿 (4)
𝐹 𝑁 𝛼𝑘 -∑ ∑ 𝑥 𝑖𝑗𝑘 • 𝑎𝑟𝑒𝑎 𝑘 • 𝑅 𝑁 𝑖 𝛼•𝜃 𝑗=1+(𝛼-1)•𝜃 𝑁 𝑖=1 ≥ 0, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (5)
𝐹 𝑃 𝛼𝑘 -∑ ∑ 𝑥 𝑖𝑗𝑘 • 𝑎𝑟𝑒𝑎 𝑘 • 𝑅 𝑃 𝑖 𝛼•𝜃 𝑗=1+(𝛼-1)•𝜃 𝑁 𝑖=1 ≥ 0, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (6) 𝐹 𝐾 𝛼𝑘 -∑ ∑ 𝑥 𝑖𝑗𝑘 • 𝑎𝑟𝑒𝑎 𝑘 • 𝑅 𝐾 𝑖 𝛼•𝜃 𝑗=1+(𝛼-1)•𝜃 𝑁 𝑖=1 ≥ 0, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (7) ∑ ∑ 𝑎𝑟𝑒𝑎 𝑘 • 𝑝 𝑖 • 𝑥 𝑖𝑗𝑘 𝐿 𝑘=1 𝑀 𝑗=1 ≥ 𝐷 𝑖 , 𝑖 = 1, … , 𝑁 (8)
∑ ∑ 𝑥 𝑖𝑗𝑘 𝑗∉𝐼 𝑖 𝐿 𝑘=1 = 0, 𝑖 = 1, … , 𝑁 (9)
𝑥 𝑖𝑗𝑘 ∈ {0,1}, 𝑖 = 1, … , 𝑁, 𝑗 = 1, … , 𝑀, 𝑘 = 1, … , 𝐿 (10)
F N 𝛼𝑘 = {F N 𝛼𝑘 ∈ 𝑅 + | 𝐹 𝑚𝑎𝑥 ≥ F N 𝛼𝑘 ≥ 𝐹 𝑚𝑖𝑛 }, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (11) F P 𝛼𝑘 = {F P 𝛼𝑘 ∈ 𝑅 + | 𝐹 𝑚𝑎𝑥 ≥ F P 𝛼𝑘 ≥ 𝐹 𝑚𝑖𝑛 }, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (12) F K 𝛼𝑘 = {F K 𝛼𝑘 ∈ 𝑅 + | 𝐹 𝑚𝑎𝑥 ≥ F K 𝛼𝑘 ≥ 𝐹 𝑚𝑖𝑛 }, 𝑘 = 1, … , 𝐿, 𝛼 ∈ Ω (13) Ω = {𝛼 ∈ 𝑁 * |𝛼 • 𝜃 ≤ 𝑀, 𝜃 ∈ 𝑁 * } (14)
The objective function presented in eq. 1 evaluates crop schedule profits and fertilization costs. Besides crop budgets already cover fertilization costs, the objective function also incurs these costs to minimize the number of consumed resources.

The constraint set in eq. 2 prevents adjacent areas to hold the same crop family scheduled in the same period and the constraint set in eq. 3 ensures that crops from the same family shall not be scheduled in sequence on each crop field.

The constraint set in eq. 4 prevents scheduling more than one crop in the same period, in other words, it is a spatial restriction. Fertilization balances are established in equations 5, 6 and 7, following the concept of surface nutrient budget. The production requirements for each crop are evaluated by the constraint set in eq. 8. The constraint set in eq. 9 ensures that crop scheduling happens just in the proper planting period, denying allocation outside this window. The decision variables 𝑥 𝑖𝑗𝑘 are Boolean type and each one represents the schedule of the crop "i" in the period "j" on the field "k" and the fertilization variables (F N 𝛼𝑘 , F P 𝛼𝑘 and F 𝐾 𝛼𝑘 ) are real values.

In order to enlighten the mathematical model and its constraints, a sample solution of the CRP is presented in Table 1.

Table 1.

A crop schedule solution attained by the proposed model for the CRP. The left side presents the crop schedule in a 24 periods interval (each period corresponds to a 15-day interval). The right side presents the related families to the sequence.

Lily

Rose Lily

Lettuce

CROPS FAMILIES Periods

The constraint set in eq. 2 prevents the same family in adjacent crops. This can be verified in the sample solution in Table 1, according to the adjacent plots presented in Table 2.

Table 2. The adjacency among plots for this solution presents on the left side. Each row represents adjacent plots from the row index. Additionally, the right side shows each plot area.

Table 3 shows that the crops from the same family can be schedule again in the same plot after a period defined by 𝛽, which satisfy the constraint set in eq. 3.

Table 3. Partial view of the previous sample solution. It shows how the problem fits the constraint set.

Nutrient constraints in eq. 5, 6 and 7 are quite intuitive and rely on supplying the required inputs from the scheduled crops on the interval 𝜃 in each plot "k". Reaching production requirements follows from constraint set in eq. 8 and scheduling in the proper planting period from eq. 9. The fertilization inputs and the attained production from sample solution are not presented for assumption of easy understandability.

The presented CRP model requires data about the crop features and the market. Usual planting dates were gathered from [START_REF]Vegetables Usual Planting and Harvesting Dates[END_REF], average prices and production statistics based on [START_REF]Vegetables 2017 Summary[END_REF]. Crop budget for the established profitability were mainly taken from [START_REF]Oregon Agricultural Enterprise Budgets[END_REF] and, in some cases, updated with production and average prices from [START_REF]Vegetables 2017 Summary[END_REF]. The crop nutrient requirements were acknowledged from [START_REF] Mohler | Crop rotation on organic farms: a planning manual[END_REF].

𝛽 𝛽 4 Methodology

Evolutionary algorithms are computational procedures for solving problems, resulting from the iterative application of heuristic techniques, which makes these algorithms capable of promoting the search for a solution in a huge space of possibilities, in a very flexible way [START_REF] Knowles | Multiobjective Problem Solving from Nature: From Concepts to Applications[END_REF]. Evolutionary Algorithms work with a population formed by a set of individuals. The coding of the individuals (representation of possible solutions to the problem) is the most important step of the algorithm. The adaptation function allows assigning each element of the search space (individual) a value that is used as a measure of performance. In optimization problems, this function incorporates all aspects of the objective function. During the generations, this population is evaluated. The most suitable ones tend to be selected and can undergo modifications in their characteristics through crossover and mutation operators, generating descendants; finally, the more adapted go to the next generation, resulting in individuals getting more fit, while others tend to disappear. Some characteristics of these methods deserve to be highlighted: flexibility, generality, ability to escape from great locations, ability to deal with complex problems for which it is not possible or difficult to obtain a detailed description, as well as being less susceptible to form or continuity [START_REF] Bäck | Evolutionary Computation 1: Basic Algorithms and Operators[END_REF] [20] [START_REF] Coello | Evolutionary Multiobjective Optimization: Current and Future Challenges[END_REF].

For the development of the approach proposed in this work, we chose the Genetic Algorithm (GA) [START_REF] Michalewicz | Genetic Algorithms + Data Structures = Evolution Programs[END_REF], because it is a classic method, one of the first Evolutionary Algorithms proposed in the literature and easily adapted to any problem, besides being extensively applied in the different areas of engineering. The functionalities of the algorithm proposed in this work were implemented with the concern of respecting the physical and operational characteristics of the problem. The coding technique and the main operators to approach this work are described ahead.

Encoding and Decoding

The binary decision variables in the CRP might appear rather suitable for binary encoded genetic algorithms, but evaluations showed that performance decreases on mid-term and long-term planning. Also, the binary GA requires hybrid approaches to improve performance [START_REF] Santos | Crop rotation scheduling with adjacency constraints[END_REF] [START_REF] Aliano | Metaheuristics for a crop rotation problem[END_REF]. The proposed GA based on the real encoding technique presented ahead.

The left side of Fig. 3 presents the data structure which holds all the feasible combinations of crop, period and plot. The right side presents a chromosome that holds two memory positions, one has a real number and the second one is an integer index related to the previous data structure.

Table 4.

Real-encoded information in the GA. The left side presents a support data structure with all the possible crop allocations and the right side shows the chromosome filled with floating numbers.

Decoding process initializes with sorting the chromosome according to the real numbers. The highest priority returns an index that represents the first crop to be allocated, then the process continuous until there are no more possible allocations. Table 5. Decoding process. The chromosome is sorted according to the floating values. The second column of the chromosome holds the index related to the support data structure. Decoding follows down this priority list until there is no more availability in the schedule.

The fertilization variables have direct correspondence to the chromosome values, apart from the schedule priority list and without requiring any decode technique.

Genetic Operators

The presented GA uses the tournament selection as a technique to sort out individuals from the population. Tournaments happens with "k" random selected individuals; the winner is the highest fitness among them and shall be placed in the mating pool, which will generate a new offspring.

The tested crossover and mutation operators are presented in [23] [24]. A brief description of Laplace Crossover technique and Power Mutation follows ahead.

• Laplace Crossover: 1) Random numbers: 𝑢 𝑖 , 𝑟 𝑖 ∈ [0,1]

2) 𝛽 𝑖 = { 𝑎 -𝑏 • log(𝑢 𝑖 ) , 𝑟 𝑖 ≤ 0.

Overview of the Proposed Genetic Algorithm

A random procedure generates the initial population. At first, the current population is the same as the initial population. Then, the genetic operators (selection, crossover and mutation) produce the so-called new population. The mixed population combines the current population and the new population, sorting individuals according to their fitness. Individuals with high fitness among the mixed population are selected to fill the current population and compose the next generation. The code was developed in C language to reach high efficiency. Fig. 1 summarizes the characteristics of the GA. The coding and decoding generate solutions which already fit constraint set in eq. 2, 3, 4 and 9, decreasing infeasibility workarounds. To increase performance, there is a penalty method applied to the nutrient balance constraints.

Discussion on Results

The proposed model has been evaluated by a deterministic approach and by the presented GA. Glop linear solver, which is a powerful linear optimization solver developed by Google, handled the deterministic evaluations. According to Table 2, achieving exact solutions requires exhaustive computational time even in mid-term planning.

Performance of the GA was tested with 48 periods and with 72 periods in comparison with the Glop solver results. In each case, the GA ran 20 times to attain the presented results, this stochastic procedure requires several executions to validate its proper performance. The initialization process generates 1000 random individuals. Each further generation has 400 new individuals. The current population and the new population generate the mixed population, which is sorted out based on fitness and the best individuals from this population shall be inserted in the current population.

Crossover and mutation probability are 1.00 and 0.15, respectively. The total number of generations is 300. These parameters were selected after tuning and provided the best results for the tested instances.

This research presents the proposed model as potential tool for farm management. In order to ensure the reliability of this mathematical approach, there should be optimization strategies able to manage sizeable instances of the CRP. The proposed GA and the deterministic approach are complementary techniques in this work. The first one turned out to be quite suitable for the large data analyzes, whereas the exact searches may seem extremely high-cost alternatives.

The results in table 6 and table 7 analyzes the proposed in terms of repeatability and stability, which are very important to ensure that the solution attained in the large instances are reasonable. Although the computational time required by the GA exceeds Glop's one, the performance of the Glop solver is undoubtedly better in the small-size problems, besides GA provides sub-optimal solutions, and so, the proposed evolutionary approach is designed to work together with deterministic methods.

Conclusion

The proposed mathematical approach for the CRP enlightens agriculturally sustainable appeal without turning back on profits as the main goal of the problem. The developed algorithm achieved prospect results in a reasonable computational time and provided quality solutions for the large instances of the CRP.

This research acknowledges the potential for further developments related to CRP in sustainability and multi-objective approaches. Field crop information will increase fast and become more reliable due to observation and monitoring techniques which are driven by Precision Agriculture. Consequently, crop nutrient balance holds an important role in further developments. The real-encoding technique presents to be very powerful in the long-term analyzes of the CRP, avoiding exceeding computational time related to the infeasibility workarounds. Establishing the relationship between crop scheduling and nutrient balance in the same model is the main contribution of this research and requires several agricultural cases to acknowledge the extent of this approach on the farm management.

Fitness ($) Gap

Based on the presented results, the proposed mathematical model is reliable and provides quality solutions through the deterministic approach or the proposed GA, each according to the data complexity. The proposed approach for this scheduling problem and the characteristics of the GA presented, such as the real-encoding technique, could suit many applications in the Industry 4.0.

Fig. 1 .

 1 Fig. 1. A flowchart presents the proposed genetic algorithm structure.

  • 𝑁 𝑓 : number of crop families; • 𝐹 𝑝 : the set of crops from the family p; • 𝑀: number of periods; • 𝐿: number of plots; • 𝑎𝑟𝑒𝑎 𝑘 : available tillage area of each plot k (acre); • 𝑙 𝑖 : profitability of crop 'i' ($); • 𝑡 𝑖 : crop 'i' production cycle, from planting date to the harvesting; • 𝑝 𝑖 : crop 'i' average production per acre; • 𝐼 𝑖 : crop planting interval {𝐼 𝑖 1 , … , 𝐼 𝑖 𝑆 𝑘 : adjacent plots of the plot k; • 𝐹 𝑁 𝛼𝑘 , 𝐹 𝑃 𝛼𝑘 , 𝐹 𝐾 𝛼𝑘 : Quantity of nitrogen, phosphorus and potassium fertilizer applied in the plot k over the interval 𝛼; • 𝑅 𝑁 𝑖 , 𝑅 𝑃 𝑖 , 𝑅 𝐾 𝑖 : requirement of nitrogen, phosphorus and potassium per area unit for crop 'i'; • 𝐹 𝑚𝑖𝑛 , 𝐹

𝑛 };

• 𝐷 𝑖 : demand for crop 'i' (units / M periods);

• 𝑚𝑎𝑥 : Fertilization limits; • 𝑐 𝑁 , 𝑐 𝑃 , 𝑐 𝐾 : fertilization costs ($ / unit); • 𝛽: sequence crop restriction factor; • 𝜃: interval of fertilization balance.

Table 6 .

 6 Results from executions in a 60-crop database over 48 periods (2-year schedule) and 7 plots.

Table 7 .

 7 Results from executions in a 60-crop database over 72 periods (3-year schedule) and 7 plots.


Rui Azevedo

Luís B Palma

Hermínio Duarte-Ramos
email: hduarteramos@gmail.com

Paulo Gil

Intelligent HCI Device for Assistive Technology

Keywords: Human-in-the-Loop (HuIL) Control, Assistive Technology (AT), Microelectromechanical Systems (MEMS), Microcontrollers

This paper presents a new intelligent Human-Computer Interaction (HCI) device for Assistive Technology. The developed device can be used as a mouse or as a gamepad, moving a part of the human body, typically the head, for hands-free computer access tasks. The state-of-the-art hardware uses an Advanced RISC Machine 32-bit microcontroller and a microelectromechanical 9-degree motion sensor, including a System-in-Package (SiP) accelerometer, gyroscope and magnetometer. The hardware/software device with a human-inthe-loop controller can be identified as a Cyber-Physical intelligent system to be incorporated in the "Industry 4.0" trend. Experimental results reveal that the embedded controller of the HID device allows the improvement of the user's performance, decreasing the effort and the execution time of the hands-free computer tasks.

Introduction

The main goal of this work is to test an intelligent Cyber-Physical System (CPS) device, developed in [START_REF] Antunes | Controlo em Sistemas de Interação Humano-Máquina[END_REF], with several microelectromechanical sensors, to allow hands-free computer access for persons with disabilities. The simplified architecture of the designed Human Interface Device (HID) is presented in Fig. 1. The human operator performs head movements, with the HID module, to use the computer, replacing the traditional mouse, keyboard and gamepad/joystick. The HID module dimensions are only: 42  21  9 mm.

State-of-the-art hands-free human-computer assistive devices nowadays in the market are the Quha Zono gyroscopic mouse [START_REF]Quha Assistive Technology Product Catalog[END_REF] and the GlassOuse assistive device [START_REF] Assistive | [END_REF]. There are also other non-gyroscopic commercial solutions, not so accurate, since are based on vision, eye-gaze and speech recognition technologies, such as the Tobii Dynavox PCEye Plus [START_REF] Dynavox | PCEye Plus, eye tracker for individuals facing mobility challenges[END_REF] and the AbleNet TrackerPro [START_REF] Trackerpro | [END_REF] assistive devices. Recent advances in brain-computer interaction and human-computer accuracy, using microelectrodes implanted into the brain [START_REF] Nuyujukian | Cortical control of a tablet computer by people with paralysis[END_REF], have the disadvantages and potential risks of being intrusive experimental techniques. The assistive devices mentioned are mainly based on open-loop electronic systems.

The proposed approach in this work is focused on a closed-loop assistance scheme with error feedback, to assist the operator and reduce human effort. It is based in the Human Adaptive Mechatronics (HAM) concept [START_REF] Laurikkala | Predicting operator's cognitive and motion skills from joystick inputs[END_REF], [START_REF] Suzuki | Assist control and its tuning method for haptic system[END_REF], and consists in adding an embedded electronic assistance-controller, tuned from a simplified mathematical model that represents the main characteristics and handicaps of each individual user with disabilities. Hence, the main research question of the PhD thesis [START_REF] Antunes | Controlo em Sistemas de Interação Humano-Máquina[END_REF], in which the proposed prototype is described in detail, is: "Can performance of a system integrating the human being be improved, if an assistance-controller is used, that takes into account the behaviour and constraints of the human operator and the machine?", and the research hypothesis, related to this part of the PhD work, is: "If following (Point-To-Point) tasks are performed, then linear modelling is a feasible method for the design of assistance control systems".

The new HID module, with an "intention" feedback PID controller, tuned by a nature-inspired metaheuristic algorithm, is able to sense the head/body movements. The developed hardware/software allows a human user to perform, hands-free, the computer's cursor movements, the gamepad/joystick movements, the virtual keyboard and the mouse/gamepad/joystick buttons activation, the single click and double click mouse configurations, and also the "drag and drop" and "scroll" functionalities.

The electronic device incorporates the 32-bit ARM Cortex-M4 NXP-Freescale MK20DX256VLH7 microcontroller and the 9-degree MEMS motion sensor TDK InvenSense MPU-9250 [START_REF]InvenSense MPU-9250 Product Specification Revision[END_REF], [START_REF]MPU-9250 Nine-Axis (Gyro+Accelerometer+Compass) MEMS Motion Tracking Devices[END_REF]. The field of application for the developed device includes the assistance for several motor and sensorial impairments, resulting from tetraplegia, quadriplegia, multiple sclerosis, cerebral palsy, amyotrophic lateral sclerosis, cranioencephalic and vertebral-medullary injuries, muscular dystrophy, the carpal tunnel or the Rett syndrome, where a human user has reduced or no control of his hands, and is commonly unable to use a traditional Human-Computer Interface.

Relationship to Industrial and Service Systems

The fourth Industrial Revolution results nowadays from the integration of innovative production systems and information technologies. The concept of Additive Manufacturing (AM) [START_REF] Meboldt | Industrializing Additive Manufacturing -Proceedings of Additive Manufacturing in Products and Applications -AMPA[END_REF] is considered the new paradigm of the "Industry 4.0" and the "Economy 4.0" trends.

The developed hardware/software can be considered as an intelligent Cyber-Physical System (CPS) with a dedicated human-in-the-loop assistance architecture, based on the human operator "intention", providing several applications for Assistive Technology and allowing the increase of the user's autonomy and productivity in computer access for service, industrial systems, active learning and entertainment, complying the Web Content Accessibility Guidelines (WCAG). The aim of this work is to make a valid contribution, adopting CPS, to a more inclusive society.

Assistance-Controller Project

According to previous works [START_REF] Suzuki | Assist control and its tuning method for haptic system[END_REF], one of the simplest human controller model is a Proportional and Derivative (PD) controller, defined as

𝐻(𝑠) = 𝐾 dh 𝑠+𝐾 ph 𝑇 N 𝑠+1 e -𝜏𝑠 . ( 1
)
where 𝐾 dh and 𝐾 ph are the derivative and the proportional gains of the human brain controller, 𝑇 N is a neuromuscular time constant and 𝜏 is a reaction time-delay. These parameters change according to human status and the operation. The experiments conducted by Suzuki et al. [START_REF] Suzuki | Assist control and its tuning method for haptic system[END_REF] revealed that the overcoming progress from learning could be considered a slow and time-variant process. Hence, the time-constant 𝜏 can be discarded for identification purposes, and the simplified human model becomes

𝐻(𝑠) ≈ 𝐾 dh 𝑠+𝐾 ph 𝑇 N 𝑠+1
.

(

) 2
The human controller 𝐻(𝑠) transfer function [START_REF]Quha Assistive Technology Product Catalog[END_REF], can be rewritten in the form

𝐻(𝑠) ≈ 𝐾 1 𝑠+𝑎 h 𝑠+𝑏 h 𝑎 h , 𝑏 h , 𝐾 1 ∈ ℝ + . ( 3
)
𝑠 = -𝑎 h and 𝑠 = -𝑏 h are the zero and the pole of 𝐻(𝑠), and 𝐾 1 a system constant. The proportional gain 𝐾 ph and the derivative gain 𝐾 dh of the human PD controller are defined by

𝐾 ph = 𝑎 h 𝐾 1 𝑏 h . ( 4
)
𝐾 dh = 𝐾 1 𝑏 h . ( 5
)
The human neuromuscular time-constant 𝑇 N is obtained by

𝑇 N = 1 𝑏 h . ( 6
)
With the 𝐻(𝑠) (3) static gain 𝐾 0 computed by

𝐾 0 = 𝑎 h 𝐾 1 𝑏 h = 𝐾 ph . ( 7
)
A simplified block diagram of the human-in-the-loop (HuIL) system, formed by the computer access hardware that features the cursor movements (𝑚 0 𝑠 ⁄ ), for each axis, and the 𝐻(𝑠) human user model is presented in Fig. 2, with the Operating System (OS) Enhanced Pointer Precision (EPP) option disabled. .

Resulting, from (3) and ( 8)

𝐺 x,y (𝑠) = 𝑚 0 𝐾 1 (𝑠+𝑎 h ) 𝑠 2 +(𝑚 0 𝐾 1 +𝑏 h )𝑠+𝑚 0 𝐾 1 𝑎 h . ( 9
)
Which leads, for each axis, to a second-order stable linear system, formed by two negative real poles and a zero, admitting the condition

2𝑚 0 𝐾 1 (2𝑎 h -𝑏 h ) ≤ 𝑏 h 2 + (𝑚 0 𝐾 1 ) 2 .
The resulting 𝐺 x,y (𝑠) system can be identified, for each axis, from an Auto-Regressive Linear Model with Exogenous Input (ARX)1 discrete time structure, with 𝑛𝑎 = 2, 𝑛𝑏 = 2 and 𝑛𝑑 = 1, neglecting the delay, which corresponds to

𝐺 x,y (𝑠) = 1+𝑇 z 𝑠 (1+𝑇 p 1 𝑠)(1+𝑇 p 2 𝑠) . ( 10
) With 𝑎 h = 1 𝑇 z . ( 11
)
𝑏 h = 𝑇 p 1 +𝑇 p 2 -𝑇 z 𝑇 p 1 𝑇 p 2 . ( 12
)
𝐾 1 = 𝑇 z 𝑚 0 𝑇 p 1 𝑇 p 2 . ( 13
)
Resulting, from (3), ( 11), ( 12) and ( 13)

𝐻 x,y (𝑠) = 𝑇 z 𝑚 0 𝑇 p 1 𝑇 P 2 ( 𝑠+ 1 𝑇 z 𝑇 p 1 +𝑇 p 2 -𝑇 z 𝑇 p 1 𝑇 p 2 +𝑠 ) . ( 14
)
The parameter 𝑚 0 = 0.1458 of ( 8) is experimentally obtained, being proportional to the computer's cursor speed, moving in one screen direction (X) with the electronic HID device configured for mouse operating mode. The 𝑚 0 constant corresponds to the inverse of the time duration for the cursor to reach the half screen distance using one axis (from normalized positions 1 to 0, or normalized positions -1 to 0). This value represents the effortless horizontal head rotation angle, affecting the proportional gain of the computer's cursor dynamics and the also the resulting closed-loop system transfer function 𝐿 HCR , obtained in (17). The cursor speed was settled to correspond to an (effortless) horizontal head rotation (yaw) of ±10°= ± 𝜋 18 ⁄ rad, without EPP. It was experimentally verified that when 𝑚 0 is too high, it could lead to instability. Nevertheless, this effect could still be compensated, by adjusting directly from the software, the computer's mouse speed. On the other hand, if 𝑚 0 is too low, the Point-To-Point (PTP) time response increases, resulting in a worse performance.

The block diagram of the closed-loop computer access assistance system (Fig. 3), includes a Proportional-Integral-Derivative (PID) controller, tuned by the Firefly Algorithm (FA), described in section 3.1. The control architecture is based on a cascade control scheme, that uses the operator's cursor "intention". The PID controller 𝐶 PID with anti-windup is tuned by the nature-inspired optimization FA algorithm, which takes into account the estimated human model 𝐻(𝑠), an expert operator specification and a set of constraints. 16) From ( 14), ( 15) and ( 16), the resulting transfer function 𝐿 HCR (𝑠), for each axis, is obtained by

𝐿 HCR x,y (𝑠) = 𝐾 p 𝑇 z 𝑇 p 1 𝑇 p 2 ( 𝑠+ 1 𝑇 z 𝑇 p 1 +𝑇 p 2 -𝑇 z 𝑇 p 1 𝑇 p 2 +𝑠 )(1+ 1 𝑇 i 𝑠 +𝑇 d 𝑠) 𝑠+𝐾 p ( 𝑇 z 𝑇 p 1 𝑇 p 2 ( 𝑠+ 1 𝑇 z 𝑇 p 1 +𝑇 p 2 -𝑇 z 𝑇 p 1 𝑇 p 2 +𝑠 )+𝑚 0 )(1+ 1 𝑇 i 𝑠 +𝑇 d 𝑠)
.

(17)

Controller Tuning with the Firefly Algorithm

The goal for the PID controller tuning is to find the values of 𝐾 p , 𝐾 i and 𝐾 d that minimize the specification for 𝐿 HCR (𝑠), given by

𝑚𝑖𝑛{0.3|𝑜𝑣𝑒𝑟𝑠ℎ -10| + 0.5|𝑟𝑖𝑠_𝑡𝑖𝑚𝑒 -3| + 0.2|𝑠𝑒𝑡_𝑡𝑖𝑚𝑒 -𝑟𝑖𝑠_𝑡𝑖𝑚𝑒 -2|}
which characterize the step response of 𝐿 HCR (𝑠) with the (expert) operator A, without EPP. The variables 𝑜𝑣𝑒𝑟𝑠ℎ, rise_time and set_time stand respectively for percentage overshoot, rise time and settling time.

The gain constraints of the PID assistance-controller in the three-dimensional search space (𝑑 = 3) take into account that the proportional and the integral gains 𝐾 p and 𝐾 i are not greater than 0.5 and the derivative gain 𝐾 d is less than 0.02, so that the step response of the human-computer control system does not become too fast (which could generate instability), and to allow to minimize the rise time and the settling time.

As the attractiveness 𝛽 of a given firefly in the search space is proportional to its luminous intensity, recognized in the environment by the adjacent fireflies, 𝛽 will be obtained by

𝛽 = 𝛽 0 e -𝛾𝑟 2 . ( 18
)
With 𝛾 representing the light absorption coefficient.

At distance 𝑟 = 0, the attractiveness of the firefly is given by 𝛽 0 . The Cartesian distance 𝑟 𝑖𝑗 in the 𝑑-space between two fireflies 𝑖 and 𝑗, located at 𝒙 𝑖 and 𝒙 𝑗 , is defined by

𝑟 𝑖𝑗 = ‖𝒙 𝑖 -𝒙 𝑗 ‖ . ( 19
)
That is

𝑟 𝑖𝑗 = √∑ (𝑥 𝑖,𝑘 -𝒙 𝑗,𝑘 ) 2 𝑑 𝑘=1 . ( 20
)
When a given firefly 𝑖, at the search space location 𝒙 𝑖 , is attracted to another brighter firefly 𝑗, its new 𝒙 𝑖 ′ position in the 𝑑-Euclidean space is computed, at iteration 𝑖𝑡𝑒𝑟, by 𝒙 𝑖 ′ = 𝒙 𝑖 + 𝛽 0 e -𝛾𝒓 𝑖𝑗 2 (𝒙 𝑗 -𝒙 𝑖 ) + 𝛼 0 𝜐 𝑖𝑡𝑒𝑟 𝐳 𝑖 .

∝ 0 ∈ [0,1] is the randomization parameter, 𝒛 𝑖 is a 𝑑-dimensional array of random numbers with uniform distribution, and 𝜐 ∈ [0,1] is the mutation coefficient employed. The Proportional-Integral-Derivative (PID) controller tuning method is described in the Algorithm 1 [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF], [START_REF] Yang | Cuckoo Search and Firefly Algorithm[END_REF], taking into account the light absorption 𝛽 0 = 1.0, the randomization parameter 𝛼 = 0.2 and the mutation coefficient 𝜐 = 0.98.

Results and Critical View

The obtained experimental results concerning the PTP tasks time duration 𝑇 PTP for human participants A to E are presented in Table 1, without and with the assistancecontroller, using the Mahony et al. sensorial fusion methodology [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF], [START_REF] Mahony | Complementary filter design on the special orthogonal group SO(3)[END_REF]. Best results are depicted in green.

It can be confirmed that the assistance-controllers improve the human-computer performance in the PTP tasks with the developed hands-free HID device, resulting in a decrease of 44 % in the average PTP time duration.

Conclusions and Further Work

The hardware/software device with the user assistance-controller improves the handsfree human-computer performance, allowing the replacement of computer access interfaces that require manual actions/movements. The relevance of the human hands can be proven by the Motor Homunculus model (Fig. 6), initially proposed in 1937 by Penfield & Boldrey [START_REF] Penfield | Somatic motor and sensory representation in the cerebral cortex of man as studied by electrical stimulation[END_REF], consisting of a representation of how the human body would be, if the size of its parts grew in proportion to the area of the brain related with the correspondent motor functions in the central sulcus of the frontal cortex. The controller of the developed hardware/software device can be configured, in a customized way, for the specific needs of each human user, allowing the adjustment of the PID parameters and also their storage in the EEPROM of the HID microcontroller, using a C# configuration application, also built. The ability to customize the computer access device, tailored to the behaviour and the relevant capabilities of each human user, modelled by a simplified mathematical transfer function, represents a novel trend in Assistive Technology (AT).

The hardware/software developed in [START_REF] Antunes | Controlo em Sistemas de Interação Humano-Máquina[END_REF] for Assistive Technology, which can be shown working, may in future be used for other relevant areas, such as biomedical engineering, rehabilitation and sports medicine, physiotherapy and interface systems engineering.

Fig. 1 .

 1 Fig. 1. HID module placed in the head (left). Simplified architecture of the Human-Computer Interface device in USB connection mode (right).

Fig. 2 .

 2 Fig. 2. Block diagram of the Human-Computer Interface (HCI), in "mouse" mode.The transfer function 𝐺 x,y (𝑠) of the resulting closed-loop system, that aggregates the human operator and the HID device, is

Fig. 3 .

 3 Fig. 3. Cascade control architecture, for each axis, of the hands-free computer access system.The transfer function 𝐿 HCR (𝑠) of the resulting closed-loop system, for each screen axis (X and Y), is obtained by

(

Fig. 5 .

 5 Fig. 5. PTP target signals (𝑟 X , 𝑟 Y ) and correspondent responses (𝑋, 𝑌), for X and Y axis of human participant B: without controller (top); with the controller (bottom).

Fig. 6 .

 6 Fig. 6. The Motor Homunculus 3D model.

The MATLAB ident command from the System Identification toolbox, provides a Graphical User Interface (GUI), allowing the identification of this type of ARX structure by Least Squares (LS), returning the 𝑇 P 1 , 𝑇 P

and 𝑇 z parameters of the equivalent closed-loop resulting model 𝐺 x,y (𝑠)[START_REF] Meboldt | Industrializing Additive Manufacturing -Proceedings of Additive Manufacturing in Products and Applications -AMPA[END_REF], for each screen axis, that aggregates the operator dynamics and the HID device.

An example of the Firefly Algorithm (FA) evolution, with 20 fireflies and 10 iterations, is presented in Fig. 4, taking into account the human operator A 𝐻(𝑠) model ( 14), with 𝑇 p 1 = 0.777, 𝑇 p 2 = 0.874 and 𝑇 z = -0.245.

Intelligent HCI Device for Assistive Technology 161

The proposed nature-inspired optimization methodology ensures improved performance in the hands-free task movements of the computer's cursor, as shown in the experimental results.

Discussion of Results and Critical View

A total of 500 Point-To-Point (PTP) random sequences were performed using the headmounted HID device, connected to the computer via USB cable, without and with the Proportional-Integral-Derivative (PID) controllers. Each PTP task, with the cursor and the EPP option enabled, consists in following a random PTP target in the computer screen, at distance, in pixels (px), 𝐴 = {300px, 400px, 500px, 600px, 700px}, and with the target diameter 𝑊 = 100 px, whose position on the screen changes every 10 seconds.

The PID controller, with anti-windup, is tuned from the Firefly Algorithm. The choice for the number of fireflies, iterations and attempts took into account a good convergence-speed compromise. The FA algorithm results revealed convergence with 20 fireflies, 10 iterations and 5 attempts, resulting in an average time duration of 5 minutes and 30 seconds, using a Skylake dual-core processor Intel® Core i5-6300U, running at 2,40 GHz, and the simplified mathematical model for each human participant (A to E), consisting of an equivalent Auto-Regressive Model with Exogenous (ARX) structure, with 𝑛𝑎 = 2, 𝑛𝑏 = 2 and 𝑛𝑑 = 1, neglecting the delay ( 14), with an operating frequency of 50 Hz (𝑇 𝑠 = 20 ms), the filter coefficient 𝑁 PID = 10 and the anti-windup time constant 𝑇 t = 1 s.

The discrete implementation of the real-time, anti-windup, PID controller structure into the embedded microcontroller of the Human Interface Device (HID) module was programmed in C programming language. Examples of acquired signals during Point-To-Point (PTP) random tasks, performed by the human participant B, using the handsfree HID developed device, without and with the PID assistance controller, are depicted in Fig. 5.