Jad Darrous

Shadi Ibrahim

Understanding the Performance of Erasure Codes in Hadoop Distributed File System

Keywords: Erasure coding, HDFS, Performance evaluation

Replication has been successfully employed and practiced to ensure high data availability in large-scale distributed storage systems. However, with the relentless growth of generated and collected data, replication has become expensive not only in terms of storage cost but also in terms of network cost and hardware cost. Traditionally, erasure coding (EC) is employed as a cost-efficient alternative to replication when high access latency to the data can be tolerated. However, with the continuous reduction in its CPU overhead, EC is performed on the critical path of data access. For instance, EC has been integrated into the last major release of Hadoop Distributed File System (HDFS) which is the primary storage backend for data analytic frameworks (e.g., Hadoop, Spark, etc.). In this work, we measure and compare the performance of data accesses in HDFS under both replication and EC. Our analysis indicates that EC is a feasible solution for dataintensive applications and it can outperform replication in many scenarios. Furthermore, we demonstrate that it is the block placement algorithm in HDFS that mostly impacts the performance under EC.

Introduction

Replication has been successfully employed and practiced to ensure high data availability in large-scale distributed storage systems [START_REF] Ghemawat | The Google File System[END_REF][START_REF] Memishi | Fault Tolerance in MapReduce: A Survey[END_REF][START_REF] Shvachko | The Hadoop Distributed File System[END_REF]. Moreover, replication can be leveraged to improve data access performance under high load. While important, replication becomes very expensive [START_REF] Hadoop | HDFS Erasure Coding[END_REF][START_REF] Rashmi | A "Hitchhiker's" Guide to Fast and Efficient Data Reconstruction in Erasure-Coded Data Centers[END_REF][START_REF] Sathiamoorthy | XORing Elephants: Novel Erasure Codes for Big Data[END_REF][START_REF] Zhang | Introduction to HDFS Erasure Coding in Apache Hadoop[END_REF], not only in terms of storage cost (given the rapid growth in data size -IDS estimates data to grow from 64.2 ZettaByte (ZB) in 2020 to 181 ZB by 2025 [START_REF] Reinsel | The World Keeps Creating More Data -Now, What Do We Do with It All[END_REF]) but also in terms of network cost (half of all cross-rack traffic is attributed to replication traffic in production clusters at Facebook and Microsoft [START_REF] Chowdhury | Leveraging Endpoint Flexibility in Data-Intensive Clusters[END_REF]) and hardware cost (expensive high-speed storage devices, i.e., SSDs and DRAMs, are increasingly deployed in storage systems to meet low latency data-intensive applications [START_REF] Kim | Towards SLO Complying SSDs through OPS Isolation[END_REF]).

While already popular for archived data in peer-to-peer systems [START_REF] Kubiatowicz | OceanStore: An Architecture for Global-Scale Persistent Storage[END_REF][START_REF] Rodrigues | High Availability in DHTs: Erasure Coding vs. Replication[END_REF] and cold data [START_REF] Haeberlen | Glacier: Highly Durable, Decentralized Storage despite Massive Correlated Failures[END_REF][START_REF] Subramanian Muralidhar | F4: Facebook's Warm BLOB Storage System[END_REF], in recent years, erasure coding (EC) has been progressively performed on the critical path of data access [START_REF] Rashmi | EC-Cache: Load-Balanced, Low-Latency Cluster Caching with Online Erasure Coding[END_REF][START_REF] Zhang | Efficient and Available In-Memory KV-Store with Hybrid Erasure Coding and Replication[END_REF]. For example, EC has been employed in in-memory storage systems on cached * Work done while at Inria, IMT Atlantique, LS2N (hot) data [START_REF] Rashmi | EC-Cache: Load-Balanced, Low-Latency Cluster Caching with Online Erasure Coding[END_REF]. The main reasons of such an adoption are that (1) EC can provide the same fault tolerance guarantee as replication but with much lower storage cost [START_REF] Rodrigues | High Availability in DHTs: Erasure Coding vs. Replication[END_REF], and (2) the cost of data encoding and decoding is progressively reduced, thanks to the release of storage acceleration libraries, e.g., ISA-L [START_REF]Intel Intelligent Storage Acceleration Library Homepage[END_REF], that allow EC operations to run at CPU speed (e.g., encoding throughput using only a single core is 5.3 GB/s for Intel Xeon Processor E5-2650 v4 [21]).

As a natural result of such an adoption, EC has been integrated in Hadoop Distributed File System since its last major release [START_REF] Hadoop | HDFS Erasure Coding[END_REF] (i.e., HDFS 3.0.0). HDFS is the primary storage back-end for data analytic frameworks (e.g., Hadoop [START_REF]Apache Hadoop[END_REF], Spark [START_REF]Apache Spark[END_REF], Flink [START_REF]Apache Flink[END_REF], etc.). Several studies have been conducted to understand the performance of MapReduce applications under EC (with contiguous and striped data layouts [START_REF] Darrous | Enabling Data Processing under Erasure Coding in the Fog[END_REF][START_REF] Darrous | Is it Time to Revisit Erasure Coding in Data-Intensive Clusters?[END_REF][START_REF] Yao | EC-Shuffle: Dynamic Erasure Coding Optimization for Efficient and Reliable Shuffle in Spark[END_REF][START_REF] Zhang | Does erasure coding have a role to play in my data center?[END_REF]), reduce the data reconstruction times under failures [START_REF] Li | Degraded-First Scheduling for MapReduce in Erasure-Coded Storage Clusters[END_REF], and improve data locality [START_REF] Li | On Data Parallelism of Erasure Coding in Distributed Storage Systems[END_REF][START_REF] Li | Parallelism-Aware Locally Repairable Code for Distributed Storage Systems[END_REF]. For example, Darrous et al. [START_REF] Darrous | Is it Time to Revisit Erasure Coding in Data-Intensive Clusters?[END_REF] provide an in-depth study on the impact of EC with striped data layout on the task runtimes and MapReduce application performances under different system configurations (e.g., HDD, DRAM, 1 Gbps and 10 Gbps network, etc.). Although important, the aforementioned studies target MapReduce applications -where the clients (i.e., the tasks) reside inside Hadoop cluster and are processing the data -and do not consider multi-tenant storage systems.

In an effort to complement existing efforts on understanding the performance of data-intensive applications under EC, in this paper, we study the performance of (concurrent) data accesses (write and read applications) under erasure coding and replication in a complementary and contrast approach. We experimentally study their behaviors through extensive micro-benchmarks on the Grid'5000 [START_REF] Balouek | Adding Virtualization Capabilities to the Grid'5000 Testbed[END_REF] testbed. We report our findings pertaining to the performance of write and read applications when varying the size of data and the number of concurrent accesses. Our results show that the write throughput under EC approaches that under replication and even outperforms it under high concurrency. Moreover, EC -by leveraging parallel reads from multiple disks -can deliver a 4.95x higher throughput than replication for a single read. However, under high concurrent reads, the performance of EC is impacted by the imbalanced load across nodes. We hope that the insights drawn from this paper will enable a better understanding of the performance of data-intensive applications under EC and motivate further research in adopting and optimizing EC in distributed storage systems.

The remainder of this paper is organized as follows. In Section 2, we present how EC is implemented in HDFS. The experimental methodology is explained in Section 3. Section 4 presents the experiments under write operation while Section 5 presents the different sets of experiments under read operation. Section 6 discusses the related work. Finally, Section 7 concludes this paper.

Erasure codes in HDFS

EC was officially integrated into HDFS in 2018, in the third major release of Hadoop ecosystem [START_REF] Hadoop | HDFS Erasure Coding[END_REF]. Hadoop community has adopted the XOR and Reed-Solomon erasure codes (RS) with striped layout (i.e., one logical block is represented physically by (𝑛) chunks -usually distributed on multiple disks -and then used to compute (𝑘) parity chunks. Accordingly, any (𝑛) out of (𝑛 + 𝑘) chunks is sufficient to rebuild the original data block). EC with striped layout is efficient for small files [START_REF] Hadoop | HDFS Erasure Coding[END_REF] which are prominent in data-intensive clusters [START_REF] Fan | DiskReduce: RAID for Data-Intensive Scalable Computing[END_REF][START_REF] Niazi | Size Matters: Improving the Performance of Small Files in Hadoop[END_REF], and requires less memory overhead when encoding and decoding data. To reduce the metadata overhead, every 𝑛 blocks -belonging to the same file -are grouped into an EC group. Blocks of the same EC group are usually placed on the same set of nodes. EC policy is defined by a scheme and the size of striping cell. The scheme consists of the number of data and parity chunks alongside the code algorithms (XOR and Reed-Solomon are supported). The striping cell size determines the granularity of data access and buffer size (1 MB by default). The default policy under EC is "RS-6-3-1024k", which means that Reed-Solomon (RS) codes are used with 6 data chunks and 3 parity chunks and the encoding/decoding operations are performed with a striping cell of 1 MB.

Methodology Overview

We conducted a set of experiments to assess the impact of data access pattern and concurrent data access when HDFS operates under replication (REP) and erasure coding (EC).

Testbed

Our experiments were conducted on the French scientific testbed Grid'5000 [START_REF] Balouek | Adding Virtualization Capabilities to the Grid'5000 Testbed[END_REF] at the site of Nantes. Two clusters -Econome and Ecotype -are used for the experiments with 21 and 40 machines, respectively. Each machine is equipped with two Intel Xeon E5-2660 8-core processors, 64 GB of main memory, and one disk drive (HDD) at 7.2k RPM with 1 TB. The machines of each cluster are connected by 10 Gbps Ethernet network. The two TORs switches of both clusters are connected with four 40 Gbps links. The machines run 64-bit Debian stretch Linux with Java 8 and Hadoop 3.0.0 installed. All the experiments have been done in isolation on the testbed, with no interference originated from other users. Econome runs HDFS (one node hosts the NameNode (NN) and the remaining 20 nodes act as DataNodes (DNs)) while Ecotype hosts the clients. To exclude the impact of the local disks at the client side, we store the dataset in the main memory.

HDFS configuration

HDFS block size is set to 256 MB (similar to [START_REF] Dinu | RCMP: Enabling Efficient Recomputation Based Failure Resilience for Big Data Analytics[END_REF][START_REF] Yildiz | Enabling fast failure recovery in shared Hadoop clusters: Towards failure-aware scheduling[END_REF]) and the replication factor is set to 3 (default value). For EC, we use the default EC policy in HDFS, i.e., 𝑅𝑆 (6, 3) scheme with a cell size of 1 MB.

Benchmarks

We used the read and write operations to measure the performance of HDFS. These operations are issued from the clients using the HDFS commands hadoop fs -get and hadoop fs -put, respectively. The test files are generated by the dd command from /dev/urandom as an input source.

Metrics

The throughput at the client side is the main metric used to measure the performance of read/write operations from/to HDFS. For one client, it is the amount of data read/written per second. In the case of concurrent reads and writes (by multiple clients), the average throughput per client is shown alongside the standard deviation of the clients' throughput.

We also profile disk and network I/O of the nodes using the python library psutil [6] version 5.4.8.

Cost of Adding New Data

Write dataflow in HDFS

Why studying the write performance under EC? Data is generated at an extreme rate. To benefit from this data, it is usually stored -and then analyzed -in data-intensive clusters. Recent studies revealed that hundreds of terabytes of data are ingested every day in data-intensive clusters [START_REF] Chowdhury | Leveraging Endpoint Flexibility in Data-Intensive Clusters[END_REF]. Those data, known as data inputs, vary significantly in their sizes. For example, traces from production Hadoop cluster [START_REF]CMU Hadoop traces[END_REF] reveal that the size of input data varies from 3 GB to 13 TB. Furthermore, these input data are populated frequently by multiple concurrent users and applications. For example, one computation requires 150 pipelined jobs to complete [START_REF] Dinu | RCMP: Enabling Efficient Recomputation Based Failure Resilience for Big Data Analytics[END_REF]; hence, the output of a job is used as an input for the later one. This motivates us to study the cost of adding data under erasure coding compared to replication. We further study the impact of concurrency on the performances of both erasure coding and replication. What is the impact of EC on write performance? Writes are pipelined under replication (i.e., the client sends the data block to the first node which, in turn, pipelines it to the second and then the third). However, under EC, the client encodes the data and then writes the data and parity chunks directly to all the nodes in parallel. This parallelism can provide higher throughput but results in more data sent by the client stressing the link between the clients and the HDFS cluster. Regardless of the HDFS storage policy, the client splits the file into blocks equal to the HDFS configured block size. Under replication, and for each block, the client contacts the NN to get a list of DNs equal to the replication factor (e.g., 3 nodes for 3-way replication) to write the block to. The client streams each block to the first DN in the corresponding list which in turn pipelines the data to the second one and so on.

On the other hand, the write dataflow under EC is a bit different; for 𝑅𝑆 (6, 3), after splitting the file into blocks, every 6 blocks are grouped into an EC group. For each EC group, the client contacts the NN to get a list of 9 DNs. To encode this group of blocks, the blocks are encoded sequentially; each block is split into cells (e.g., 1 MB), then the client encodes every 6 cells to generate 3 parity cells. The client sends these 9 cells to the 9 DNs, and then continues to do the same with the remaining cells, and repeats the same process for the remaining blocks in the group. Then, do the same for the remaining EC groups. In conclusion, three (number of replicas) DNs will be contributing to a write operation at a time under REP, while 9 DNs (𝑛 + 𝑘) will be simultaneously writing data under EC. It is important to mention that blocks are written sequentially to HDFS (the second block is sent once the first block is completely stored in a DN, persisted in disk or buffered in the memory). Moreover, when the data is stored completely in all DNs (at least buffered in the memory), the write operation is considered successful.

Results of single write

Figure 1a shows the write throughput of one client with different file sizes. The write throughput is 1.11x to 1.3x higher under REP compared to EC when increasing the file sizes from 256 MB to 20 GB, respectively.

When writing 20 GB file, REP achieves 363 MB/s write throughput, while EC obtains a throughput of 277 MB/s which is 76% of the throughput under REP. This can be explained due to the larger amount of data which is sent from the client to HDFS in EC (i.e., 30 GB accounts for the original and parity data) compared to REP (i.e., 20 GB) and the "low cost" data pipelining under REP. First, data stays in the buffers of the DNs when received from the client or pipelined from other DNs, hence, it is not synchronized to disks directly. As a result, data transfer time -to HDFS -strongly depends on the amount of data sent to the cluster. Given that the throughputs of inter-cluster traffic (network traffic between the clients and the HDFS cluster) are 444 MB/s and 399 MB/s under EC and REP respectively, the file under EC takes 35% more time to be transferred to HDFS compared to REP (the difference in transfer times between EC and REP is clearer when increasing the file size). Second, (buffered) data can be pipelined as soon as they reach the first DN because the intra-cluster (data transferred inside the HDFS cluster) throughput is almost two times the inter-cluster throughput (i.e., intra-cluster throughput is 801 MB/s), thus replicating (transferring) data inside the cluster does not introduce an extra cost. Finally, more data is persisted to disk under REP compared to EC (i.e., 21 GB under EC while it is 29 GB under REP); which in turn slightly slows down the write operation under REP.

Note that, by default, Hadoop client utilizes only one thread to send the data and therefore the inter-cluster throughput may not be fully utilized. In case of REP, data is communicated between two nodes only and thus the inter-cluster throughput is limited to 399 MB/s. In case of EC, despite that the client is sending data to 9 DNs, inter-cluster throughput is limited to 444 MB/s because only one core is responsible for encoding and sending data (in case of 20 GB, we observe that one core is always at 100% utilization). Observation 1. Unlike under replication, more data (i.e., parity data) goes from the client to HDFS cluster under EC. This extra parity data (0.5x of the original data in our configuration) results in a loss of throughput compared to replication as data are pipelined (replicated) inside the cluster with minimal cost.

2. Disk, inter-network and intra-network throughput for concurrent writing of 5 GB files by 40 clients.

Results of concurrent writes

When increasing the number of concurrent clients to 5, EC starts to approach the throughput of REP (Figure 1b). Interestingly, with 40 clients, EC has a throughput 2 times that of REP for files equal and larger than 5 GB as depicted in Figure 1c.

In case of writing 5 GB files by 40 concurrent clients, REP achieves 32 MB/s write throughput, while EC obtains a throughput of 83 MB/s which is 2.59x higher than the throughput under REP. Disks are saturated under both EC and REP (the average disk throughput is 110 MB/s and 124 MB/s, respectively). As more data is persisted to disks under REP compared to EC (i.e., 464.4 GB under REP while it is 160.5 GB under EC), writing data to disks under EC is almost 2.5 times faster. On the other hand, the throughput of the intra-cluster network is low (i.e., 183 MB/s to 233 MB/s when the files are larger or equal to 5 GB), thus write performance will further degrade under REP. The low disk and intra-cluster performances will limit the arrival rate of data to HDFS cluster under REP (the arrival rate is 1162 MB/s under REP compared to 4535 MB/s under EC). As shown in Figure 2, the network performance (inter and intra) drops once the buffers are filled and data starts to be synchronized to disks while still not completely replicated (transferred) to other DNs.

To further explain the impact of concurrency on the performance of write under both EC and REP, Figure 1d shows the aggregated write throughput when writing a data set of 100 GB while varying the number of clients and therefore the files sizes. When the number of concurrent clients is 5, REP achieves a throughput of 1410 MB/s while it is 1340 MB/s under EC. The aggregated throughputs increase by 1.94x and 3.25x when increasing the number of clients to 10 and 20 clients under EC, respectively. This increase is mainly due to the increase in the inter-cluster throughputs; On the other hand, the aggregated throughputs increase 1.2x and 1.5x under replication when increasing the number of clients to 10 and 20 clients, respectively. Here, the high cost of data pipelining limits the scalability of HDFS under REP when increasing the number of clients. Specifically, buffers under replication are filled faster (original data and replicated data) compared to EC. Thus, the arrival rate of data to a DN is limited by the disk throughput (125 MB/s under EC and 128 MB/s under REP). However, as a DN under replication is receiving data from the clients and from other DNs, this will reduce the inter-cluster throughput and prolong the writing time under REP compared to EC. Observation 2. Under concurrent writes, in contrast to under EC, the performance of write workloads under replication is limited by the intra-cluster transfer and disk contention (as in total, the amount of written data is larger under REP compared to EC). Hence, due to the high network and memory cost of data pipelining, EC can even outperform replication -under heavy concurrent writes. This gap is more obvious for large files (e.g., in our experiments, the throughput under EC is at least 2x the throughput achieved under replication when 20 clients are writing 10 GB file each (Figure 1c).

Implications

Replication is not only costly in terms of storage requirements but also contributes to the problems of oversubscribed networks [START_REF] Dinu | RCMP: Enabling Efficient Recomputation Based Failure Resilience for Big Data Analytics[END_REF] and poor disk performance [START_REF] Dinu | RCMP: Enabling Efficient Recomputation Based Failure Resilience for Big Data Analytics[END_REF] in production clusters (e.g., traces from production clusters at Facebook and Microsoft pointed out that replication results in almost half of all cross-rack traffic [START_REF] Chowdhury | Leveraging Endpoint Flexibility in Data-Intensive Clusters[END_REF]). Applying EC as an alternative to replication does not only reduce storage cost and disk overhead but also results in lower network traffic (up to 50%). In conclusion, given its performance, EC is feasible for write applications, more importantly, EC sustains high throughput for write operations under high concurrency.

Reading Data under EC

Why studying the read performance under EC? Reading input data is the first step for any data-intensive application. Many studies have discussed the importance of optimizing the read step (i.e., reading input data) to improve the job execution time -mainly through targeting high data locality [START_REF] Chen | Design insights for MapReduce from diverse production workloads[END_REF][START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF][START_REF] Ibrahim | Maestro: Replica-Aware Map Scheduling for MapReduce[END_REF][START_REF] Zaharia | Delay Scheduling: A Simple Technique for Achieving Locality and Fairness in Cluster Scheduling[END_REF]. Therefore, it is important to understand read performance under EC. Typically, multiple jobs are running concurrently (to achieve higher cluster utilization) therefore we need to evaluate concurrent reads. Also, it is possible that multiple jobs (could belong to different users) use the same files as input, thus we study concurrent reads of the same file under both storage policies. What is the impact of EC on read performance? Clients reading data under replication can contact any replica and retrieve all the data from only one replica, this can provide load balance between the nodes hosting the replicas under concurrent reads. On the other hand, to read a block under EC, the client needs to contact the data nodes and read the data block in parallel. This might result in higher read throughput if the network bandwidth is higher than disk read bandwidth, but may cause an imbalance in the load as nodes hosting parity chunks may serve fewer requests.

Read dataflow in HDFS

To read a file in HDFS, the client contacts the NN to get the addresses of the DNs hosting each block of the file. For each block, the DNs are ordered by the network distance to the client and are randomized for the nodes that have the same distance. The client then reads the blocks in order i.e., sequentially, one by one. Reading one block of data (without any failure) has the same cost under both EC and REP, as the same amount of data (equal to the block size) will be read from disk and transferred over the network. However, under EC the data block can be read in parallel from multiple DNs (6 DNs in our configuration). This parallelization results in better performance especially when the network bandwidth is higher than the disk bandwidth [START_REF] Ganesh Ananthanarayanan | Disk-Locality in Datacenter Computing Considered Irrelevant[END_REF].

Results of single read

The read throughput for different file sizes by one client is depicted in Figure 3a. As we expect, EC performs better than REP. Specifically, the read throughput is 1.4x to 4.95x higher under EC compared to REP when increasing the file size from 256 MB to 20 GB, respectively. When reading 20 GB file, REP achieves 116 MB/s read throughput, while EC obtains a throughput of 575 MB/s which is 4.95x higher than the throughput under REP. The reason behind this is that an EC block is read from 6 nodes, while under REP, a block is read from a single node. Therefore -as the performance of the read operation is limited by the disk -6 disks are leveraged in parallel under EC, while a single disk is used under REP. Hence, the throughput of the read workload is strongly co-related with the performance of the disks. Therefore, the performance gap between EC and REP is clearer when increasing the file sizes, especially when the disk throughput under EC increases.

Under replication, a block will be read sequentially from the disk, hence, the disk will be exploited efficiently (i.e., the throughput of the disk is between 85 MB/s and 100 MB/s most of the time). On the other hand, under EC, a block will be read in parallel from multiple nodes (6 nodes in our configuration), therefore, each node reads a fraction of the block (one data chunk of 43 MB) and thus cannot leverage the full throughput of the disk. As we can see in Table 1, when the size of the file is 256 MB, the disk throughput is limited by the physical block size (43 MB/s) under EC. However, when the client reads all the blocks of an EC group (6 blocks in our configuration), the disk throughput of an individual node is fully utilized as 6 chunks belong to different blocks will be read sequentially from the disk. For instance, when the file size is 20 GB, the average disk throughputs (of active disks) for both EC and replication are 83.4 MB/s and 86.6 MB/s (the disk throughput is higher than 80 MB/s most of the time), respectively. Hence, this explains the 4.95x performance difference.

Observation 3. Compared to 3-way replication, a block (256 MB) is scattered on more DNs under EC. Therefore, reading a file under EC has a clear advantage compared to under REP, as multiple disks are leveraged in parallel.

Results of concurrent reads

Reading distinct files. When increasing the number of concurrent clients to 5 -each reading a distinct file, the gap between EC and REP starts to decline as shown in Figure 3b. Interestingly, with 40 concurrent clients, REP achieves slightly higher throughput (clearer when the file size is 1 GB) as shown in Figure 3c. Moreover, we notice that the performance of read operations exhibits noticeable variation in-between different clients under both REP and EC.

When 5 clients are reading 5 files, each has a size of 20 GB, REP achieves 97 MB/s read throughput, while EC obtains a throughput of 202 MB/s which is 2.08x higher than that under REP. The average disk throughput of active DNs is 60 MB/s and 86 MB/s under EC and REP, respectively. The average disk throughput under EC decreases -compared to the one in case of 1 client (i.e., 83.4 MB/s) -which is due to data skew. As shown in Figure 4, although most of the DNs are serving data, some of them (around 6 DNs) serve more clients a time which in turn results in exploiting the full capacity of the disk (i.e., 125 MB/s) on those nodes. On the other hand, the majority of nodes were serving 1 client most of the time, thus low utilization of the disk. We observe that disk throughputs are above 100 MB/s during 50% and 18% of the read workload time under REP and EC, respectively. Hence, this, in addition to the increase in the number of contributing nodes (active disks), explains the decrease in the performance difference from 4.95x to 2.08x between EC and REP.

When 40 clients are reading 40 files, each has a size of 20 GB, we observe that HDFS suffers from obvious load imbalance under both REP and EC. Under REP, at least one disk was not active during 25% of the read workload (see Figure 3d). Under EC, some DNs still serve more clients compared to other DNs which in turn results in high contention on the disks and thus low disk throughputs; disk throughputs are above 100 MB/s during 40% and 2% under REP and EC, respectively. As a result, the aggregated disk throughput under EC is 1352 MB/s and is 1350 MB/s under REP, which explains the 2% difference. In addition, this load imbalance results in high variation in the latency of the read workload in-between clients, therefore, clients which are served by nodes under heavy load need more time to read their files.

Reading the same file. Figure 5a shows the average read throughput per client when 5 clients are reading the same file concurrently. With 20 GB file, the average read throughput is 506 MB/s under EC while it is 122 MB/s under REP, thus 4.14x faster under EC. To explain the difference, we inspect the amount of read data and sent data at the cluster level. Under EC, the cluster has an aggregated disk throughput of 537.6 MB/s and network throughput of 2840 MB/s: each block under EC is read once from the same nodes, cached and then sent to the 5 clients; hence, network throughput is 5.2x higher than the disk throughput. However, under REP, the cluster disk throughput is 322.5 MB/s and the network has a throughput of 652.8 MB/s, 1.7x higher than the disk. Here, a block is served by the 3 DNs at the same time: as a block is available on 3 DNs, the DNs receive requests -from different clients -to read the block, accordingly, the 3 DNs will read the data from the disk, cache it and then sent it to 5 clients (one DN serves 1 or 2 clients). This explains why the total amount of read bytes is 52.1 GB (2.6x the size of the file) under REP and 20 GB under EC.

When increasing the number of concurrent clients to 40, we notice that clients maintain their average throughput under REP (119.2 MB/s) while it drops to 195.3 MB/s under EC, as shown in Figure 5b. Under REP, each block is served by 3 different DNs, so a block will be served by a new set of DNs (3 DNs) while the previous block is still being served to clients (new block is requested once a client out of the 40 clients finishes fetching the current block). Consequently, disk throughput increases under REP from 322 MB/s with 5 clients to 367.6 MB/s and thus the network throughput increases to 5079 MB/s. On the other hand, 6 DNs will be continuously serving the 40 clients: each DN is responsible for one original chunk. This will saturate the network bandwidth of the DNs and cause variations in-between clients. Moreover, as the reading of a new block can start when all the chunks are received by at least one client, the same 6 DNs will be serving two blocks (belong to the same EC group) simultaneously and this will reduce the number of disk read requests for the new block and at the same time reduce the amount of data sent to clients through network and further increase the variation. Even worse, the performance degradation at the disk-level and the network-level (and the resulted variation in the throughput in-between clients) will be amplified when DNs are serving two successive EC groups. We observe that the number of disk requests drops by almost 50% (from 50 to 24) when increasing the number of clients from 5 to 40 clients. Consequently, the disk throughput drops from 537.6 MB/s to 229.2 MB/s-the disk throughputs of active DNs are less than 50 MB/s for almost 80% of the read workload time. This, in addition to the contention at the network due to the high concurrency, results in a network throughput of only 8048 MB/s and a variation of 5.7% in the throughput in-between clients. Figure 5c depicts the read throughput of the same file when increasing the number of clients. REP maintains the same throughput per client, however, we can clearly see how the throughput under EC drops with increasing the number of concurrent clients. Table 2 shows the cluster-level network and disk throughput. REP maintains a stable disk throughput when increasing the number of clients and relative increase of network throughput with the number of clients. However, under EC, the disk throughput drops while the network throughput increases sub-linearly.

Observation 4. The goal of presenting data as EC groups is to reduce metadata overhead in the NN, but it results in a high imbalance in data distribution across DNs compared to replication. This, in turn, causes stragglers when high concurrent reads are performed to either distinct files or the same file, and therefore, it reduces the advantage of parallel chunks reads. Specifically, DNs which are serving data continuously to multiple clients will exhibit low disk throughput due to high contention on disk (distinct files) or high contention at the network (same file).

Implications

The advantage of striped block layout under EC is to enable parallel data reads. However, the use of EC groups (to reduce metadata overhead) reduces this advantage under concurrent reads (by different clients). Because it introduces imbalance in the data access across DNs which, in turn, causes stragglers. Therefore, mitigating stragglers under concurrent access is essential. While late binding can be employed (as in EC-Cache [START_REF] Rashmi | EC-Cache: Load-Balanced, Low-Latency Cluster Caching with Online Erasure Coding[END_REF]), it comes with the cost of extra network transfer. Moreover, intelligent data access can be used (as in EC-Store [START_REF] Michael Abebe | EC-Store: Bridging the Gap between Storage and Latency in Distributed Erasure Coded Systems[END_REF]), but this method requires historical analysis of data access. On the other hand, reads under EC benefit more from OS caches than replication as data blocks are always read from the same node, which might be beneficial for iterative applications. In conclusion, the previous experiments demonstrate the feasibility and the advantages of using EC for read operations in large-scale clusters.

Related Work

Zhang et al. [START_REF] Zhang | Does erasure coding have a role to play in my data center?[END_REF] show the performance benefits of EC when the intermediate and output data of MapReduce applications are encoded compared to 3-ways replication. Yao et al. [START_REF] Yao | EC-Shuffle: Dynamic Erasure Coding Optimization for Efficient and Reliable Shuffle in Spark[END_REF] demonstrate that encoding the intermediate data of MapReduce like applications can provide faster recovery after failures especially for multi-stage applications. Darrous et al. [START_REF] Darrous | Is it Time to Revisit Erasure Coding in Data-Intensive Clusters?[END_REF] study the impact of EC with striped data layout on the performance of Sort, Wordcount and Kmeans applications under different system configurations (e.g., failures, different back-end storage devices and network configurations, etc.). Among their findings, they also show that disk contention caused by chunks distribution affects the performance of data-intensive applications under EC. Li et al. [START_REF] Li | On Data Parallelism of Erasure Coding in Distributed Storage Systems[END_REF][START_REF] Li | Parallelism-Aware Locally Repairable Code for Distributed Storage Systems[END_REF] propose to collocate original data chunks and parity chunks on the same nodes to improve data locality (increase the number of local map tasks). This work is complementary to the aforementioned studies in that it specifically targets understanding the performance of (concurrent) write and read operations when the clients reside outside Hadoop cluster and provides a more detailed analysis of the network and disk performance under EC.

Conclusion

Efficient data access in large-scale storage systems is an important problem with the current rate of data generation. Erasure codes are increasingly deployed in many storage systems as cost-efficient alternative to replication. In this paper, we experimentally evaluate the performance of data read and write operations in HDFS under both replication and erasure coding. Our findings can be summarized as follows: (1) Write operations under replication have higher throughput than under EC for a single client. However, when increasing the number of concurrent clients, the throughput under EC approaches that under replication and even outperforms it under high concurrency. In addition, applying EC as an alternative to replication does not only reduce storage cost and disk overhead but also results in lower network traffic.

(2) When reading data, EC can leverage parallel reads from multiple disks and deliver a 4.95x higher throughput than replication for a single client. However, under high concurrency, the performance of EC is impacted by the imbalanced load across nodes, which is caused by the distribution of chunks and the design of EC group.

As future work, we plan to implement a new data placement strategy that targets balancing the data load across nodes and study its impact on the performance of read/write operations and MapReduce applications in shared Hadoop clusters.

Figure 1 .

 1 Figure 1. Write under REP and EC.

Figure 3 .

 3 Figure 3. Read distinct files under REP and EC.

Figure 4 .

 4 Figure 4. Disk throughput of the DNs during the concurrent read of 5 files by 5 clients each of 5 GB under EC.

Figure 5 .

 5 Figure 5. Read the same file under REP and EC.

Table 1 .

 1 Disk read throughput (MB/s) of active DNs in case of single client read.

		File size mean 25% 50% 75%	max
		256 MB	35.42 26.62 42.00 42.75	43.00
		1 GB	71.28 57.00 71.25 85.50	86.00
	EC	5 GB	74.59 65.51 85.33 85.51	85.67
		10 GB	83.76 85.35 85.45 85.61 102.80
		20 GB	83.39 77.44 85.39 86.85 102.80
		256 MB	85.33 85.33 85.33 85.33	85.33
		1 GB	85.42 85.33 85.33 85.42	85.67
	REP	5 GB	85.19 85.46 85.67 85.67 102.40
		10 GB	88.44 85.33 85.55 89.40 102.80
		20 GB	86.63 83.19 85.50 92.86 102.60

Table 2 .

 2 Average aggregated cluster network and disk throughput (MB/s) in case of concurrent read of the same 20 GB file by the number of clients.

	Number of	Network through-	Disk through-
	clients	put (MB/s)		put (MB/s)
		EC	REP	EC	REP
	5	2840.3 652.8	537.6	322.5
	10	5836.8 1331.5	552.8	374.2
	20	7009.7 2648.2	384.1	376.1
	40	8048.7 5079.1	229.2	367.6

Acknowledgments

This work was supported in part by the ANR KerStream project (ANR-16-CE25-0014-01), the Stack/Apollo connect talent project, and the Inria Project Lab program Discovery (see https:// beyondtheclouds.github.io). Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER, and several Universities as well as other organizations (see http://www.grid5000.fr/).


Radita Liem

Shadi Ibrahim

Revisit Data Partitioning in Data-intensive workflows

As data volumes increase (e.g., in 2020 alone, 64.2 ZettaByte of data was generated and replicated [START_REF] Reinsel | Worldwide global datasphere forecast, 2021-2025[END_REF]), data-intensive workflows are now executed across multiple platforms from the Edge to Cloud and HPC, which is known as computing continuum. The intention behind this execution model is to exploit data locality and reduce data movement by executing data close to their sources while taking into consideration the computation capacities of different platforms.

Data-intensive workflows usually consist of multiple computation stages where data flows among them. For example, MapReduce jobs consist of two consecutive stages (Map and Reduce) and the output of the map stage is transferred and used as an input for the Reduce stage [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF], [START_REF] Jin | The MapReduce Programming Model and Implementations[END_REF]. The completion time of a stage strongly depends on the finishing time of the last task in this stage. Accordingly, much efforts have focused on reducing the variation in task execution times in the same stage by launching speculative copies of slow tasks [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF] or evenly partitioning data across tasks to mitigate data skew [START_REF] Ibrahim | Leen: Locality/fairness-aware key partitioning for mapreduce in the cloud[END_REF], [START_REF] Gufler | Handling data skew in mapreduce[END_REF], [START_REF] Ibrahim | Handling partitioning skew in mapreduce using leen[END_REF], [START_REF] He | Ds 2: Handling data skew using data stealings over high-speed networks[END_REF].

With the proliferation of data-intensive workflows such as IoT applications, machine learning applications, and deep learning applications; data skew remains a bottleneck in data analytic frameworks. Many research efforts have been dedicated to mitigating data skew by evenly distributing data across tasks while considering data locality [START_REF] Ibrahim | Leen: Locality/fairness-aware key partitioning for mapreduce in the cloud[END_REF], [START_REF] Ibrahim | Handling partitioning skew in mapreduce using leen[END_REF] and relying on high-speed networks [START_REF] He | Ds 2: Handling data skew using data stealings over high-speed networks[END_REF], or by partitioning the data based on the computation capacities of the nodes [START_REF] Chen | Libra: Lightweight data skew mitigation in mapreduce[END_REF]. However, most of these efforts focus on only two-stage applications and don't consider network and I/O heterogeneity within or across platforms, making them impractical when running multi-stage data-intensive workflows on heterogeneous and shared environments.

As shown in Figure 1, PageRank application exhibits a severe skew in partitioned data (i.e., shuffled data) in stage 3 and stage 4 of the application (bottom chart). This in turn results in a noticeable variation in task execution times. The variation is represented by the gap between the average and maximum task execution times (top chart). We can also observe that the variation is more obvious when the degree of parallelism is set to 24, 32, and 40. Here, it is important to note that, the data assigned to each task (i.e., task input) includes the partitioned data and the data retrieved from the storage memory (cache). § Work done while at Inria, Univ. Rennes, CNRS, IRISA, France. In this work in progress, we will showcase a comprehensive analysis of the current state-of-the-art solutions for data skew mitigation in several environments. Our experiments and evaluation comprise several data-intensive workflows running on Spark using the Grid'5000 testbed [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF]. The data-intensive workflows vary from a highly optimized WordCount application, an iterative application like PageRank, to an SQL-based decision support system benchmark, TPC-H with various sizes and configurations.

Fig. 1 .

 1 Fig.1. The skew in partitioned data (i.e., data shuffled from previous stage and read by each task as a part of its input) and its impact on task execution times: We report the total size of shuffled data and task execution times in stages 3 and 4 when running PageRank application with 2.8 GB of data on Spark cluster of 4 nodes, deployed in Rennes site of Grid'5000. The number of stages is set to 6 stages.

ACKNOWLEDGMENTS

This work is supported by the ANR KerStream project (ANR-16-CE25-0014-01). Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER, and several Universities as well as other organizations (see http://www.grid5000.fr/).


Walid Amehri

Alexandre Kruszewski

FEM-based Reachable Workspace Estimation of Soft Robots using an Interval Analysis approach

Keywords: Soft robotics, Reachable workspace, Finite-element method, Interval analysis

Considering a soft robot with installed bounded actuators, this paper studies the reachable workspace estimation problem for a pre-chosen point of interest. To this aim, the method of finite-element has been used to obtain the static model of the studied robot and two methods are proposed (forward one which is based on Newton iterative method and forward-backward one which is based on interval analysis techniques) to study this estimation problem. Various simulations with different configuration scenarios are provided to show the effectiveness of the forward-backward approach which can largely reduce the exploration time, compared to the forward method.

Introduction

Rigid robots utilize stiff materials to produce enough stiffness and are commonly used in the industry to increase productivity, which however suffer from safety problems when collaborating with human beings [START_REF] Cetinkunt | Performance limitations of joint variable-feedback controllers due to manipulator structural flexibility[END_REF]. Therefore, researchers in the field of robotics have sought a new design of robots that can provide more functionality, flexibility, and safety in order to cope with new robotic applications [START_REF] Yim | Design and rolling locomotion of a magnetically actuated soft capsule endoscope[END_REF][START_REF] Payne | An implantable extracardiac soft robotic device for the failing heart: Mechanical coupling and synchronization[END_REF]. Being inspired from natural elements [START_REF] Iida | Soft robotics: challenges and perspectives[END_REF], soft robots represents the perfect instrument to operate in uncertain and dynamic environments [START_REF] Shepherd | Multigait soft robot[END_REF] owing to their low stiffness and natural compliance. At the same time, as pointed out by [START_REF] Rus | Design, fabrication and control of soft robots[END_REF], the classical procedure of designing a soft robot needs to be guided and optimized via the evaluation of their reachable workspace, while such an important subject has been rarely examined in this community.

In the field of rigid robotics, different approaches were explored to treat the workspace estimation problem for different types of rigid robots. An optimization method was presented in [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF] to determine the boundaries of the workspace of rigid-link manipulators. A continuation method based on bifurcation theory has been stated in [START_REF] Haug | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] to estimate the workspace boundaries of different rigid-link manipulators. In [START_REF] Abdel-Malek | Analytical boundary of the workspace for general3-dof mechanisms[END_REF], a method based upon analytical analysis of Jacobian was developed to determine mechanism's singular behavior of the investigated rigid manipulator.

Also, the wrench-feasible workspace of parallel rigid robots driven by cables has been analyzed in [START_REF] Gouttefarde | Interval-analysis-based determination of the wrenchfeasible workspace of parallel cable-driven robots[END_REF] by applying an interval analysis bisection method. Such a method has been then extended to analyze the workspace of cable-suspended parallel rigid robots in [START_REF] Merlet | On the workspace of suspended cable-driven parallel robots[END_REF][START_REF] Berti | Workspace analysis of redundant cable-suspended parallel robots[END_REF]. Although the interval analysis based approaches have been successfully implemented in determining the workspace of several types of rigid robots, however they cannot be directly extended to analyze soft robots because of the complexity of the modeling of soft robots. In fact, the inputs/outputs relation (i.e., the actuation vector/the end-effector position) might be explicitly defined through its direct kinematic model for rigid manipulator, however such a relation for soft manipulator is generally defined in an implicit way through a complex static model (which will be explained in Section 3) for the purpose of modeling soft material deformation.

In the literature, different types of modeling techniques have been proposed to characterize the flexibility of soft robots. Roughly speaking, the PCC (Piece-wise Constant Curvature) modeling approach [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] is broadly employed by the soft robotic community, however the assumption of the constant curvature may not be always valid, more specifically when external loads are applied on the soft robot. The PCS (Piecewise Constant Strain) model [START_REF] Renda | Discrete cosserat approach for multisection soft manipulator dynamics[END_REF], where a finite set of strain vectors describes the soft robot, considers the material proprieties and the geometric non-linearity, which in fact can be seen as an extension of PCC approach. However, like the PCC method, the PCS approach can only be applied to model uniformly slender shapes of soft robots. FEM (Finite Element Method) is suitable to model general geometric form of soft robots [START_REF] Duriez | Control of elastic soft robots based on real-time finite element method[END_REF], and can provide an accurate model for soft robots, but with a relatively higher dimensional model.

In order to estimate the reachable workspace of soft robots with general form, this paper adopts the FEM method to describe the input-output relation of investigated soft robots. In this paper, we solve the problem of workspace estimation for the studied soft robots through 2 steps: firstly the static model used to describe soft robots is deduced by using FEM [16, p. 428], and secondly 2 methods (classical forward method based on Newton iterative approach and novel one based on interval analysis [START_REF] Jaulin | Interval analysis[END_REF]) are proposed to solve the reachable workspace estimation problem for a pre-chosen point of interest. Compared to the classical forward method (which consists of discretizing the actuators' space), the efficiency of the novel approach (named as forward-backward method in this paper) has been validated via different configuration scenarios of soft robots. The preliminary result of this work was presented in the conference [START_REF] Amehri | Fem based workspace estimation for soft robots: A forwardbackward interval analysis approach[END_REF] and the full version presented in this paper contains new results for FEM static model with external forces, self-contained results on interval analysis, novel uniform spatial grid discretization strategy, detailed algorithm with correct stop condition and comprehensive configuration scenarios of simulation to emphasize the effectiveness of the proposed forward-backward approach.

Problem Statement

Pick-and-Place operation represents the most basic task that a soft robots' end-effector must be able to perform. But, for such an operation, a soft robot in the first place must be able to reach the position of the object to be handled. In other words, the reachable workspace of the designed structure of the soft robot must include the position of that object, as indicated by the blue point in Fig. 1. Consistently, considering a particular configuration of a soft robot that is controlled via bounded installed actuators (i.e., cables), this paper investigates the evaluation of its reachable workspace.

We would like to emphasize that the cable-driven actuation is only taken as an example in this paper to illustrate the proposed method when estimating the reachable workspace of soft robots, since it is broadly employed in the community of soft robotics. However, since this paper uses the FEM modeling method, and that FEM was applied to different actuation manners, e.g. pneumatic actuation [START_REF] Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF] and fluidic actuation [START_REF] Rodríguez | Real-time simulation of hydraulic components for interactive control of soft robots[END_REF], therefore the proposed workspace estimation approach is broadly applicable to various types of actuations. In addition, using a different method of actuation will only influence the deduced mathematical FEM model, and the proposed workspace estimation approach remains applicable for other types of actuations.

To tackle the workspace problem, we denote q E ∈ R 3 as the displacement vector of the end-effector with respect to its position vector p E ∈ R 3 . Then, its reachable workspace can be defined as follows:

Definition 1. Consider a soft robot modeled by q = Φ(u, q, q), and controlled via installed bounded actuators u ∈ U, the reachable workspace W E of the end-effector p E = g(q) is a subspace of R 3 , that can be defined as:

W E = {p E | ∃u ∈ U, q s.t. Φ(u, q, 0) = 0, p E = g(q)}
According to the above definition, the reachable workspace W E contains all possible equilibrium positions for which p E can reach and stay there. We would like to emphasize that the 'stay' property is important and equivalent to the space of all feasible equilibrium points.

Generally speaking, the evaluation of a soft robot's reachable workspace can be physically realized provided that we have the real prototype of the investigated soft robot. But, for the purpose of guiding and optimizing the design of a soft robot before its actual design, it is more reasonable to evaluate its workspace in a virtual environment. To achieve this, two steps are followed: the first one consists of using FEM to establish the mathematical model of soft robots, and the second one consists of estimating the reachable workspace using two different methods.

Static Model of Soft Robots

In order to estimate the reachable workspace of a soft robot, we need first to deduce the designed soft robot's static model, and for this we adopt the Finite Element Method. This modeling approach firstly consists of spatially discretizing the deformable domain of the structure into smaller finite elements n e , through the creation of a corresponding mesh which contains a finite number of degrees of freedom n p to interpolate the behavior of the deformable domain by measuring the variation of the associated nodal displacements. In the literature, various geometries are proposed to discretize the spatial domain of a given structure, such as prisms, beams, hexahedrons and tetrahedrons [21, p. 732]. For the purpose of generality and convenience for element matrix integration, this paper selects the linear tetrahedron element as the appropriate meshing structure (Fig. 2 illustrates the implementation of the selected linear tetrahedron mesh on the trunk-like soft robot).

Since we are investigating the reachable workspace of soft robots which contains all equilibrium points, then according to Definition 1, the investigated soft robot's motion [16, p. 428] can be described by the following static model:

K N (q) = H T (q)u + F (q) (1)
where q ∈ R 3np is the nodal displacement vector, K N (q) ∈ R 3np is the vector of internal forces, H T (q) ∈ R 3np×nu is the direction of the actuators' forces matrix, u ∈ R nu is the control input vector, n u is the dimension of the input vector, and F (q) ∈ R 3np is the vector of external forces applied to the soft robot.

In the deduced static model (1), the term u represents the applied force to the soft robot which might be generated via different types of actuators. As we have explained before, FEM has been widely used to describe the material deformation subject to different types of actuators, including cable tension [START_REF] Zhang | Visual servoing control of soft robots based on finite element model[END_REF], cable displacement [START_REF] Zhang | Visual servoing control of soft robots based on finite element model[END_REF], pneumatic actuation [START_REF] Bieze | Finite element method-based kinematics and closed-loop control of soft, continuum manipulators[END_REF], fluidic actuation [START_REF] Rodríguez | Real-time simulation of hydraulic components for interactive control of soft robots[END_REF], etc. The only difference is to determine the vector function H T (q) according to different actuation methods, for the purpose of describing how the associated force u is generated by those actuators (which might be tension, displacement, pressure, etc.), and how u is applied to the soft robot. In this sense, the result presented in this paper is quite general, and is not limited only to cable-driven actuator. We use the cable-driven actuator only as an illustrative example since it is one of the widely used techniques in the community of soft robotics. Also, due to the fact that FEM uses only Newton's second law to describe the robot's motion, therefore only the forces vector can be involved in the final FEM model [START_REF] Cetinkunt | Performance limitations of joint variable-feedback controllers due to manipulator structural flexibility[END_REF]. This explains as well why FEM model contains only the nodes' displacement, but not their rotations. Assumption 1. The vector of internal forces of the structure K N (q) can be approximated, in accordance with [16, p. 430], as follows:

K N (q) = K(q) q (2)
where K(q) represents the assembly of each element's stiffness matrix:

K(q) = ne e=1 K e (q)
and the stiffness matrix [16, p. 429] of an element e is dependent of the elasticity matrix ∇, the strain matrix B e (q) and the volume V e (q) of the element e: K e (q) = Ve(q) B T e (q) ∇B e (q)dv(q).

Estimation Approaches of Reachable Workspace

In accordance with the workspace definition (see Definition 1), and with Assumption 1 being applied to the static model ( 1), the workspace of a soft robot can finally be defined by the following set of equations:

K(q) q = H T (q)u + F (q) q E = Cq p E = q E + p (0) E ( 3
)
where q E is the displacement vector of the end-effector with respect to its known initial position vector p (0) E , p E represents the position vector of the end-effector, and C is a selection matrix associated to the end-effector node coordinates.

Based on the above model, the estimation of the reachable workspace W E is then equivalent to determine, for all relevant inputs u ∈ U, the bound of the end-effector p E in (3). However, under the FEM framework, the difficulty of such a task lies in the following two reasons: firstly the nonlinearity of K(q)q in (2) and its dependence on material properties, and secondly the lacking of an analytic implicit relation such that p E = φ(u), due to the nonlinearity of the mapping φ : U → W E .

In the following, two different approaches, a forward one and a forward-backward one, will be proposed. Both of them use different numerical schemes to approximate the mapping φ in order to estimate the reachable workspace W E from the bounded input u ∈ U.

Forward method

The most trivial approach to determine the reachable workspace of a soft robot is to iteratively approximate the map φ through the discretization of the input bound U (as illustrated by Fig. 3).

To achieve this, given the information of the end-effector's position vector p (j-1)

E and the vector of nodal displacements q (j-1) (for the case j = 1, p (0) E and q (0) are respectively the initial position of the end-effector and the initial nodal displacement vectors of the investigated soft robot), we apply Newton-Raphson method to obtain the next nodal displacement q (j) corresponding to its input value u ∈ U d , with U d being the discretized set of U, which can be realized as follows: q (j) = q (j-1) -Ψ q (u, q (j-1) ) Ψ(u, q (j-1) ), ∀u ∈ U d where Ψ q (u, q) is the gradient of Ψ(u, q) with respect to q, and Ψ(u, q (j-1) ) = K(q (j-1) )q (j-1) -H T (q (j-1) )u-F (q (j-1) ). Then, the corresponding end-effector position vector can be deduced by:

p (j) E = Cq (j) + p (0)
E Finally, as depicted by Fig. 3, the estimation of the reachable workspace W E can be achieved by iterating the above process for all u ∈ U d .

In this paper, such an approach is named as the forward method because it is based on using the input vector u to estimate the end-effector's reachable workspace in a forward manner. The computation complexity of this method represents its primary limitation. Since u ∈ U ⊂ R nu , and suppose that the bound for each entry of u is discretized into n d subsets, then we will have n nu d possible combinations of input for u. In this sense, the computation complexity for the forward method will be increased if the number of inputs n u becomes bigger, and when the discretization precision is increased, it will explode exponentially, even if n u is small.

We can clearly notice that the dimension of the inputs n u , which may vary depending on various configurations of soft robots, represents the base of this method's computation complexity. However, the end-effector space dimension is invariant, always equal to R 3 , and does not depend on the configuration of the investigated soft robot. In accordance with this reflection, it is therefore logical to discretize the end-effector's space with constant dimension, as opposed to discretizing U with variant dimension, for the purpose of decreasing the computation complexity of the reachable workspace estimation problem. In this paper, we name such a method as the iterative forward-backward method, which will be explained in the following sections.

Forward-Backward method

Interval analysis introduction

The proposed forward-backward approach uses interval analysis techniques in order to estimate the workspace of soft robots. Therefore, a brief introduction of the interval analysis techniques is presented in the following. More details can be found in [START_REF] Jaulin | Interval analysis[END_REF].

Consider two real intervals

X 1 = [X 1 , X 1 ] and X 2 = [X 2 , X 2 ]
, then according to [START_REF] Jaulin | Interval analysis[END_REF], the product of X 1 and X 2 can be calculated by the following:

X 1 X 2 = min C(X 1 , X 2 ) , max C(X 1 , X 2 ) with C(X 1 , X 2 ) = {X 1 X 2 , X 1 X 2 , X 1 X 2 , X 1 X 2 }.
The product of a real matrix M ∈ R by a real interval X = [X , X ] with appropriate dimension is calculated as follows:

M X = M X , M X , M > 0 M X , M X , M < 0
However, due to the wrapping effect ( AX ⊃ {Ax | x ∈ X } ) [START_REF] Jaulin | Interval analysis[END_REF], this technique leads to an over-estimation of the resulted space which can be explained by considering the following example:

A = 1 1 0 1 , X = [-1, 0] [1, 2] then AX = [0, 2] [1, 2]
which implies that (0, 2) T belongs to AX , whereas it does not belong to the actual value set B = {Ax | x ∈ X }, as illustrated by Fig. 4. Therefore, we need to eliminate the resulted over-estimation in order to obtain the true space, and this will be implemented in the proposed forward-backward approach and explained in the following.

Forward-Backward approach introduction

In what follows, we present a brief overview of how the proposed forward-backward approach operates.

On the opposite to the forward approach which consists of discretizing the bounded input space, the forward-backward approach is based on the idea of discretizing the end-effector's space which is always constant, typically 3 (the dimension of position) or 6 (with orientation), regardless of how many actuators are installed. However, the process of discretization of the end-effector's space requires firstly the knowledge of its bound. For this, we estimate this boundary, at the (j -1)th iteration, via a linear approximation between the corresponding displacement vector q and the input vector u: K(q (j-1) )δq (j) = H T (q (j-1) )δu (j-1) + F (q (j-1) ) with δq (j)

E = Cδq (j) , p (j) E = δq (j) E + p (j-1) E , and
δq (j) = q -q (j-1) δu (j-1) = u -u (j-1)
where δq (j) and δu (j-1) represent respectively the displacement vector and the input vector with reference to the prior soft robot configuration.

Roughly speaking, based on the above approximation, we firstly implement interval analysis techniques in order to forwardly estimate the end-effector's workspace W (j) E for the related δq (j) (as illustrated by Step 1 in Fig. 5). Secondly, since we now have the information of the workspace W (j) E ⊂ R 3 , then we adopt a prescribed precision to discretize it. After that, we use the current end-effector's position p (j-1) E as a base to determine a feasible small displacement neighborhood around it (noted as W (j) Es in Step 2 of Fig. 5). Thirdly, in order to guarantee that the robot can reach this feasible small displacement neighborhood W (j)

Es , we backwardly calculate the equivalent input vector u (j) (as depicted in the domain U in Step 3 of Fig. 5). In the fourth step of the method, we calculate the new input bound U (j) using the deduced value of u (j) (illustrated as Step 4 in Fig. 5). Next, based on the obtained value of U (j) , we can forwardly compute the end-effector's new bound W (j+1)

E . If W (j)
Es and W (j+1) E both contain the feasible small displacement neighborhood of the end-effector, it implies that this neighborhood has already been visited. Such a conclusion allows us to efficiently and accurately shrink the over-estimated workspace caused by interval analysis (see Section 5 for the details). Finally, the proposed approach can explore the whole reachable workspace W E by iterating the above procedure.

In comparison with the forward approach, the forward-backward method provides the following two benefits:

1. Since the end-effector space is fixed to 3 for the position case, the forward-backward approach can significantly decrease the dimension of the discretized space in the forward method;

2. The forward approach cannot reduce the number of iterations, whereas, when discretizing the endeffector's space using a specified precision (with n d sub-spaces), the forward-backward method can eliminate the wrong sub-spaces caused by the interval estimation technique.

5 Forward-Backward method Realization

In this section, we present the details of how the forward-backward method is realized (refer to Fig. 5).

Step 1: Forward Over-estimation of W (j)

E using interval analysis

In the first step of the forward-backward approach, interval analysis techniques are used for the purpose of forwardly estimating the end-effector's workspace. To do this, we firstly need to linearly approximate the workspace defined by (3) with the knowledge of the end-effector's position p (j-1) E , of the displacement vector q (j-1) , and of the corresponding input vector u (j-1) ∈ U (j-1) , as follows: K(q (j-1) )δq (j) = H T (q (j-1) )δu (j-1) + F (q (j-1) )

δq (j) E = Cδq (j) p (j) E = δq (j) E + p (j-1) E (4) with δq (j) = q -q (j-1)
and δu (j-1) = u -u (j-1)

Next, suppose that the actuators are well installed in a manner that the investigated soft robot is always controllable [START_REF] Zheng | Controllability pre-verification of silicone soft robots based on finite-element method[END_REF], then the matrix K(q (j-1) ) is always invertible for any given value q (j-1) . Let us now introduce the matrix A as:

A(q (j-1) ) = CK -1 (q (j-1) )H T (q (j-1) )

and define the matrix A + by

A + = max{0, A}
where the operation max{•} is applied element-wise [START_REF] Zheng | Design of interval observer for a class of uncertain unobservable nonlinear systems[END_REF], and

A -= A -A +
With the above necessary elements, for the related displacement δq (j) , the over-estimated workspace W (j) E corresponding to the input bound U (j-1) can be calculated via the following theorem: Theorem 1. Given a soft robot's controllable configuration, with the input bound U (j-1) defined as

U (j-1) = [δu (j-1) , δu (j-1) ]
then the end-effector workspace W (j) E is over-estimated as

W (j) E = p (j) E ∈ R 3 |p (j) E ∈ p (j) E , p (j) E (5)
with p (j) E = A + δu (j-1) + A -δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1) E and p (j) E = A + δu (j-1) + A -δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1) E .

Proof. Using the linear approximation defined in (4), the nodal displacement vector can be directly deduced as: δq (j) = K -1 (q (j-1) ) H T (q (j-1) )δu (j-1) + F (q (j-1) )

With the defined matrix A, we can then use the above equation to write the end-effector displacement in function of the matrix A, which yields:

p (j) E = Cδq (j) + p (j-1) E = CK -1 (q (j-1
) ) H T (q (j-1) )δu (j-1) + F (q (j-1) ) + p (j-1) E = CK -1 (q (j-1) )H T (q (j-1) )δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1) E = A(q (j-1) )δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1) E Next, according to Lemma 3 in [START_REF] Zheng | Design of interval observer for a class of uncertain unobservable nonlinear systems[END_REF], the following over-estimation is obtained:

A + δu (j-1) + A -δu (j-1) ≤ Aδu ≤ A + δu (j-1) + A -δu (j-1)
By simply defining p (j) E = A + δu (j-1) + A -δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1) E and p (j) E = A + δu (j-1) + A -δu (j-1) + CK -1 (q (j-1) )F (q (j-1) ) + p (j-1)

E that leads to:

p (j) E ≤ p (j) E ≤ p (j)
E Finally, we obtain the over-estimated workspace W (j) E defined by (5).

Step 2: Feasible small displacement neighborhood W (j) Es computation

The purpose of the second step consists of determining the feasible small displacement neighborhood W (j)

Es around the current end-effector position p (j-1) E . To do this, we propose two strategies to discretize the over-estimated workspace W (j)

E calculated via

Step 1 with a specified precision to obtain a related discretized space W (j) E d

. The first discretization strategy is based on the principle that linear approximations in the FEM model are only valid within ±5% of the displacement, thus the first strategy consists of using the linear approximation of K N (q) for varying the value of p (j-1) E in each of the three Cartesian directions (x, y and z) with a ratio of ±5%. Denote l (j-1) = p (j-1) E 2 which represents the Euclidean distance between p (j-1) E and the origin of the earth frame, then such a scheme enables us to obtain the following discretized space W (j)

E d : W (j) E d = p (j) E ∈ W (j) E |p (j) E = p (j-1) E + 5%l (j-1) Γ, for Γ ∈ Z 3 (6)
The second discretization strategy is simpler and it consists of discretizing the space W (j) E d by using a globally uniform spatial grid with a globally small pre-defined precision s E d , i.e.,

W (j) E d = p (j) E ∈ W (j) E |p (j) E = p (j-1) E + s E d Γ, for Γ ∈ Z 3
In addition, at the point p (j-1)

E

, we can define its η-neighborhood as

S η (p (j-1) E ) = pE ∈ R 3 | p (j) E -p (j-1) E ∞ = η
where η is equal either to 5%l (j-1) or s E d , depending on the chosen discretization strategy. Consequently, we can finally deduce the feasible small neighborhood around p (j-1) E as follows:

W (j) Es = W (j) E d ∩ S η (p (j-1) E ) (7)
The illustration presented by Fig. 6 clearly depicts the relation between W (j) E , the discretization of its space W (j) E d , and the relative feasible small displacement neighborhood space W (j) Es .

Step 3: Backward configuration estimation

After having calculated the feasible small neighborhood W (j) Es via [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF], the third step consists of computing the feasible input vector u (j) via a backward mechanism, to guarantee that any point p(j) E ∈ W (j) Es \ p (j-1) E might be accessible by the end-effector of the soft robot (see Fig. 7). For this, with the values of p (j-1) E and q (j-1) , the associated deformation vector of the end-effector can be then calculated as: q(j)

E = p(j) E -p (j-1)

E

Consequently we can calculate the equivalent input and displacement vector by solving the optimization problem formulated below: [u (j) , q (j) ] = min

u,q q E -q(j) E 2 2 s.t. u ∈ U K(q)q = H T (q)u + F (q) q E = Cq (8)
Solving the backward mechanism (8) yields the feasible input vector u (j) with its related feasible displacement vector q (j) , which will eventually drive the soft robot's end-effector to reach p(j) E . Finally, we deduce the end-effector position associated with the optimal displacement vector q (j) as

p (j) E = Cq (j) + p (j-1) E .
It is important to state that, for each iteration (j), the inverse problem ( 8) has a unique solution, because this optimization problem is formulated inside a small neighborhood of q(j)

E , and we always solve this problem with reference to the prior iteration (j -1).

Step 4: Update input bound U (j) In the fourth step of the forward-backward approach, with the information of U (j-1) , we calculate the new input bound U (j) related to the input vector u (j) obtained in the third step.

To this aim, let us define the following diagonal n u × n u matrix:

diag(u (j) ) = diag u (j) 1 , • • • , u (j)
nu where u (j) i

(1 ≤ i ≤ n u ) represents the ith related value of u (j) , then the new input bound can be updated by the following equation:

U (j) = U (j-1) -diag(u (j) ) × 1 nu×2 (9)
with 1 nu×2 being an n u × 2 dimensional matrix whose entries are equal to 1. The presented updating scheme is illustrated in Fig. 8.

Step

5: Condition of stop/iteration

After having obtained p (j-1) E , q (j-1) , u (j-1) and U (j-1) , by applying the above presented forward-backward method, we might get p (j) E , q (j) , u (j) and U (j) . Iteratively, we can finally obtain the whole reachable workspace of the pre-chosen point of interest.

For a particular configuration p (j-1)

E

, assume that the corresponding over-estimated workspace W (j) E has been calculated via Step 1, then its related discretized space W (j) E d can be computed via Step 2, from which, according to [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF], an associated feasible small displacement neighborhood W (j) Es is deduced. An important property of Step 2 is that it allows us to eliminate certain non-feasible neighborhoods around p (j-1) E , which results in decreasing the amount of the feasible small displacement neighborhoods W (j) Es to be explored for the next iteration of the forward-backward approach, as illustrated in Fig. 6 by verifying the condition [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF]. On the other hand, given the information of W (j-1) Es and the newly calculated feasible small displacement neighborhood W (j) Es from [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF], an additional property which allows us to further shrink the space to be explored is that some of the points belonging to W (j) Es would probably have been explored by the former iterations, noted as W (j-1) E ∪ where W (j-1)

E ∪ = ∪ j-1 k=1 W (k)

Es

The above elimination procedure is effectuated in Step 5, as depicted in Fig. 9, by verifying the following condition:

W (j) Es ← S η (p (j-1) E ) \ (W (j-1) E ∪ ∩ S η (p (j-1) E )) (10)
Combining these two conditions (eq. ( 7) in Step 2 and eq. ( 10) in Step 5) allows us to shrink the space to be explored efficiently and accurately, which will result in significantly reducing the computation complexity of the proposed forward-backward approach (as shown in the simulation results Section 6).

Obviously, for any given point pE belonging to the union of the jth iteration, i.e., pE ∈ W (j) Es , if the corresponding end-effector position p(j) E belongs to the union of the (j -1)th iteration (as shown in Fig. 10), i.e., pE ∈ W (j-1) E ∪ , then we can draw the conclusion that the neighborhood p(j) E has already been visited. The above observation can be expressed as the following check-condition:

∃j ∈ N, s.t. pE ∈ W (j-1) E ∪ , ∀p E ∈ W (j) Es or equivalently ∃j ∈ N, s.t. W (j) E ∪ = W (j-1) E ∪
Consequently, we deduce that all feasible small displacement neighborhoods were explored. At this level, we can obtain the end-effector's reachable workspace via the following equation:

W E = W (j-1) E ∪ .

Realization of the algorithm

In the table below, a brief algorithm is presented in order to describe the main steps of the forward-backward approach.

Algorithm 1 Calculate {p E , q, u, U, W Es , W E ∪ } (j) , W E Require: Γ, s E d , {p E , q, u, U, W Es , W E ∪ } (j-1)
while W (j) Es / ∈ W (j-1)

E ∪ do Stop Condition. W (j) E ← f 1 ({p E , q, U} (j-1) )
Step 1: Forward estimation of W (j) E , eq. ( 5)

W (j) E d ← f 2 (p (j-1) E , Γ, s E d )
Step 2: Discretized space W (j) E d

, eq. ( 6)

W (j) Es ← W (j) E d ∩ S η (p (j-1) E )
Step 2: Feasible small displacement neighborhood W (j) Es , eq. ( 7) p(j)

E ← W (j) Es \ p (j-1) E Exclusion. q(j) E ← p(j) E -p (j-1) E {p E , q, u} (j) ← f 3 (q (j)
E , q (j-1) )

Step 3: Backward estimation of u (j) , eq. ( 8)

U (j) ← f 4 (u (j) , U (j-1) )
Step 4: New input bound U (j) , eq. ( 9)

W (j) Es ← S η (p (j-1) E ) \ (W (j-1) E ∪ ∩ S η (p (j-1) E ))
Step 5: Dismiss explored feasible neighborhood, eq. ( 10)

W (j) E ∪ ← W (j-1) E ∪ ∪ W (j) Es Append. {p E , q, u, U, W E ∪ } (j-1) ← {p E , q, u, U, W E ∪ } (j) Update Values. j ← j + 1 Next iteration. end while W E ← W (j-1) E ∪
The Estimated Reachable Workspace.

Comprehensive Simulation Results

This section presents the comparative results between the forward and forward-backward approaches for 3 different configuration scenarios for a trunk-like soft robot with a Young's modulus E = 1.8 × 10 7 P a, and a Poisson ratio ν = 0.45. The simulations were implemented on an Intel Xeon(R) with a 16-GB RAM and a 3.50 GHz processor.

We consider those different scenarios where the trunk-like soft robot is actuated by different cable routing configurations, then we implement both the forward and forward-backward approach in order to estimate the reachable workspace W E for each scenario. The results of both approaches are then superposed and illustrated by three-view drawing (Oxy, Ozy, and Ozx plans) for each scenario. For the forward approach, we propose to discretize the actuation force vector with a discretization step size of 1 Newton(unit) in order to respect the linear approximation of the FEM model and ensure small displacements.

6.1 Scenario 1: Trunk-like soft robot actuated by two cables.

In the first scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by two cables (as shown by Fig. 11a). The applied force is of the following range:

U (0) = U = [0 100] × [0 100].
Following the procedure presented in Section 4.1, we firstly apply the forward approach with a discretization number equal n d = 100, and we obtain the corresponding estimation of W E ( illustrated by blue colored points in Figs. 11b, 11c and11d) for this soft robot's configuration. In term of operations' computation complexity, the forward approach requires n nu d = 100 2 operations, and in term of time's computation complexity, it took approximately 1463 seconds to obtain the full estimation of this scenario workspace.

Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 5. Given the value of U (0) we forwardly estimate W (1) E using in Step. 1 of Fig. 5, next, we calculate a feasible small displacement neighborhood W (1) Es as indicated in Step. 2 of Fig. 5. After, as explained in Step. 3 of Fig. 5, we use the backward mechanism to compute the feasible input vector u (1) of each configuration belonging to W (1) Es . Finally, we compute U (1) as shown in Step. 4 of Fig. 5. We iterate the above process for the next iterations until the stop condition is verified, as explained in Step. 5 and illustrated by Fig. 10.

For this scenario, the forward-backward approach took 37 complete iterations and found a total of 635 feasible points in the workspace with a computation time of 757 seconds. Those obtained feasible points are depicted by red colored points in Figs. 11b, 11c and11d, and the union of all small neighborhood of those feasible points gives the estimation of workspace W E , which is represented by the gray zone with green contour. It is clear to see from the three-view drawing that the workspace estimation result obtained via forward-backward method coincides with that obtained via forward approach.

Moreover, for a complete 37 iterations it is expected to have W E d = 6 37 = 6.18 × 10 28 points to explore. However, what we actually explored is a total of W E d = 1680 points, and a final total of W Es = 635 feasible points. This significant reduction was partially due to the elimination process of inappropriate configurations from Step. 2 (See Fig. 6), but mainly due to the elimination of all feasible configurations that were already explored during Step. 5 (See Fig. 9).

Consequently, the proposed forward-backward approach can yield comparable workspace estimation precision compared to forward method, and at the same time reduce both the number of operations and the computational time required for the workspace estimation of this soft robot's configuration.

6.2 Scenario 2: Trunk-like soft robot actuated by three cables.

In this scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by three cables (as shown by Fig. 12a). The applied force is of the following range:

U (0) = U = [0 100] × [0 100] × [0 150].
First, we apply the forward approach with a discretization precision of n d = 100, and we achieve the following workspace W E estimation (illustrated by blue colored points in Figs. 12b, 12c, and12d) of this particular configuration. In term of operations' computation complexity, the forward approach requires n nu d = 100 3 operations, and in term of time's computation complexity, it will approximately require 1.46×10 5 seconds to obtain the full estimation of this scenario workspace.

Next, we apply the forward-backward interval analysis approach using the equations presented in Section 5. For this scenario, the forward-backward approach took 44 complete iterations and calculated a total of 1941 feasible points in the workspace with a computation time of 2063 seconds. Those achieved feasible points are illustrated with red colored points in Figs. 12b, 12c and12d, and the union of all small neighborhood 7 Discussion and Perspectives

In this paper, we proposed a forward-backward approach, based on interval analysis, to efficiently estimate the workspace of soft robots. Our approach uses FEM (3) as the modeling method, where the nodal displacement of elements q is used to describe the deformation of material. As a consequence, such an approach cannot directly catch the orientation information of the end-effector, and that is the reason why this paper discusses only the workspace of the end-effector's position p E in (3), by selecting C ∈ R 3×3np as a selection associated to the end-effector node coordinates.

However, we would like to emphasize that such an approach might be adapted to analyze the orientation workspace of the end-effector by simply modifying the above matrix C. Precisely, denote q E 1 , q E 2 and q E 3 as 3 independent nodal positions around the end-effector, which can be used to define a unique plane of the end-effector. Then, we can always properly choose the matrix C ∈ R 6×3np such that q E now is a function of those 3 points and contains both position and orientation information of the defined plane. It is evident that such a modification of matrix C will not influence the proposed approach to estimate the workspace of soft robots, since what we need to do is to adapt the matrix C in the workspace set [START_REF] Payne | An implantable extracardiac soft robotic device for the failing heart: Mechanical coupling and synchronization[END_REF].

In addition, the physical buckling phenomenon can occur under special combined action of cables with certain tensions, i.e., under special configuration with specific applied inputs. In other words, buckling phenomenon might not happen if the used numerical method can avoid those special configurations when calculating the workspace. In general, for the classical forward method, this phenomenon might occur which however might be avoided in the proposed forward-backward method, due to the backward mechanism (Step 3) applied to the feasible small displacement neighborhood around the current configuration (Step 2). Precisely, suppose that the current configuration (note the combined cable tension as u B , the nodal displacement as q B , and the end-effector position as p B ) is quite close to yield buckling. The proposed forward-backward method will firstly determine a feasible small displacement neighborhood around p B (Step 2), based on which it will backwardly calculate the associated actuation action around u B which can lead the end-effector of the soft robot to reach the points in the determined feasible small displacement neighborhood. In other words, since those equilibrium points are all located inside this small neighborhood, i.e., there exist feasible paths to link those points, thus the proposed forward-backward approach can always assess the workspace under this situation, even if it is quite close to buckling. This can be seen as another advantage of the proposed forward-backward method, compared to the classical forward one.

Conclusions

Using the FEM modeling approach, we have presented a classic method, namely the forward approach, which consists of discretizing the actuators' space in order to achieve the reachable workspace estimation of soft robots. In order to reduce the computation complexity of this forward method, we have proposed a novel method, named as the forward-backward approach, which consists of discretizing the end-effector's space in order to achieve an efficient estimation of the reachable workspace of soft robots.

The proposed approach was validated via various soft robot's configurations by using the forward method as a comparative reference. Through those simulations, we have shown the efficiency of the proposed forwardbackward method in reducing the complexity of both operations and time required to achieve the full estimation of the reachable workspace of soft robots, comparing to the basic forward method which explodes exponentially when the actuator's dimension increases.

Figures

Figure 1 :

 1 Figure 1: Scheme of the investigated soft robot.

Figure 2 :

 2 Figure 2: Geometric model and the related FEM model.

Figure 3 :

 3 Figure 3: Newton method based forward approach.

Figure 4 :

 4 Figure 4: Over-estimation introduced by the wrapping effect.

Figure 5 :

 5 Figure 5: Interval analysis based forward-backward method.

Figure 6 :

 6 Figure 6: Elimination of non-feasible neighborhoods.

Figure 7 :

 7 Figure 7: Estimation of input u (j) .

Figure 8 :

 8 Figure 8: Update the boundedness of input U (j) .

Figure 9 :

 9 Figure 9: Dismiss pre-explored feasible displacement neighborhoods.

Figure 10 :

 10 Figure 10: Iterative procedure and final result.

  (a) Soft Robot Structure -Two Cables. (b) Oxy plane view of W E . (c) Ozy plane view of W E . (d) Ozx plane view of W E .

Figure 11 :

 11 Figure 11: Scenario 1 -W E estimation via forward approach (blue points) and forward-backward (gray area) approach.

  (a) Soft Robot Structure -Three Cables. (b) Oxy plane view of W E . (c) Ozy plane view of W E . (d) Ozx plane view of W E .

Figure 12 :

 12 Figure 12: Scenario 2 -W E estimation via forward approach (blue points) and forward-backward (gray area) approach.

Figure 13 :

 13 Figure 13: Scenario 3 -W E estimation via forward approach (blue points) and forward-backward (gray area) approach.

Figure 14 :

 14 Figure 14: Operations and time computation complexity of the forward and forward-backward approach for workspace estimation of soft robots.

Table 1 :

 1 Operations and time computational complexity: forward vs forward-backward.

		Operations	Time (seconds)
	Scenario / Approach Forward	Forward-Backward	Forward	Forward-Backward
	First (n u = 2)	10 4	635	1463	757
	Second (n u = 3)	10 6	1941	≈ 1.46 × 10 5	2063
	Third (n u = 4)	10 8	3667	≈ 1.46 × 10 7	3733

of those feasible points yield the estimated workspace W E , which is represented by the gray zone within the green contour. Using the three-view illustration, we can clearly observe that the estimated workspace calculated via the proposed forward-backward approach conforms with that obtained through the forward approach.

Furthermore, for a complete 44 iterations it is expected to get W E d = 6 44 = 1.17 × 10 34 points to investigate. Yet, the total of points that we actually explored is W E d = 5466, with a final total of W Es = 1941 feasible points. This substantial reduction was slightly due to Step. 2 (refer to Fig. 6) elimination process of irrelevant configurations, but considerably due to the exclusion of already explored feasible configurations throughout Step. 5 (See Fig. 9).

Subsequently, the proposed forward-backward approach can achieve comparable workspace estimation precision compared to the forward method, and simultaneously reduce both the number of operations and the computational time required for the workspace estimation of this soft robot's configuration.

Scenario 3: Trunk-like soft robot actuated by four cables.

In this scenario, we consider a cable routing configuration where the trunk-like soft robot is actuated by four cables (as shown by Fig. 13a). The applied force is of the following range:

First, we apply the forward approach with a discretization precision equal n d = 100, and we obtain the corresponding estimation of W E (illustrated by blue colored points in Figs. 13b, 13c and13d) for this soft robot's configuration. In term of operations' computation complexity, the forward approach requires n nu d = 100 4 operations, and in term of time's computation complexity, it will approximately take 1.46 × 10 7 seconds to achieve the full estimation of this scenario workspace.

Next, we apply the proposed forward-backward interval analysis approach by following the procedure presented in Section 5. For this scenario, the forward-backward approach took 42 complete iterations and computed a total 3667 feasible points in the workspace with a computation time of 3733 seconds. Those obtained feasible points are depicted by red colored points in Figs. 13b, 13c and13d, and the union of all small neighborhood of those feasible points gives the estimation of the workspace W E , which is represented by the gray zone with green contour. It is clear to see from the three-view drawing that the workspace estimation result obtained via forward-backward method coincides with that obtained via forward approach.

In addition, for a complete 42 iterations it is expected to have W E d = 6 42 = 4.81 × 10 32 points to explore. However, what we actually explored is a total of W E d = 10626 points, and a final total of W Es = 3667 feasible points. This significant reduction was partially due to the elimination process of inappropriate configurations from Step. 2 (See Fig. 6), but mainly due to the elimination of all feasible configurations that were already explored during Step. 5 (See Fig. 9). Based on those results, the same conclusion can be drawn, i.e., the proposed forward-backward approach can yield comparable workspace estimation precision compared to forward method, and at the same time reduce both the number of operations and the computational time required for the workspace estimation of this soft robot's configuration.

Notes on Computation Complexity

For each particular soft robot's configuration, we compared the workspace estimation results obtained via the forward and forward-backward approach. According to the simulation results, we observe that the proposed forward-backward approach has, compared to the forward approach, significantly reduced the computation complexity in both operations and time of the workspace estimation of soft robots.

In the following table (see Table 1), we recapitulate the operations and time computation complexity of each scenario, where the computation complexity of the forward approach for both the second and the third scenario is just approximated since they required a very large time.

Visually, both the operations and time computation complexity can be illustrated respectively by Fig. 14a and Fig. 14b, where we can clearly observe the influence of increasing the inputs' dimension on the forward approach which makes the computation complexity explodes exponentially in both operations and time, whereas the computation complexity of the proposed forward-backward approach remains almost linearly stable.