Masih Hani
email: masih.hani@insa-strasbourg.fr

Denis Cavallucci

Prioritizing the initial problem choice in the inventive design process: Proposal of an FMEA-AHP-based method

Keywords: Inventive design, Complex problem, TRIZ, FMEA, AHP

Today, companies are competing to improve the innovation cycle time in order to gain pioneering advantages in the market. To achieve this goal, these companies could apply the TRIZ-based systematic approaches such as Inventive the Design Methodology (IDM). Nevertheless, one of the criticisms often leveled is that these approaches are time-consuming due to the demands for building a complete map to clarify a problem situation in the initial analysis phase of the inventive design process. As a result, the Inverse Problem Graph (IPG) method has been proposed to increase the agility of the process by starting the problem formulation from an initial problem located in the lower level of a problem situation. The application of the IPG method decreases problem formulation time by identifying a limited number of problems. However, this method does not provide any tool for prioritizing the initial problems in the process. The purpose of this article is to integrate a Failure Mode Effect Analysis (FMEA) -Analytic Hierarchy Process (AHP) based method into the IPG method to prioritize the initial problems in the rst phase of the design process. The proposal capability is nally tested through its application in a case study on the Lattice Structure.

Introduction

In today's fast-changing business environment, many companies are seeking effective means to reduce their innovation cycle time [1]. Among them, it is possible to mention the TRIZ-based systematic approaches, such as the Inventive Design Methodology (IDM), which help companies to reduce the amount of time needed to obtain the optimal solutions [2], [START_REF] Cavallucci | Designing the Inventive Way in the Innovation Era[END_REF]. Nevertheless, one of the drawbacks of these approaches is that they use the methods such as Problem Graph, Network of Problems, which require a complete map of the problem situation. This makes the inventive design process very timeconsuming [START_REF] Hani | A Method to Formulate Problem in Initial Analysis of Inventive Design[END_REF]. To solve this problem, Hani et al. [START_REF] Hani | Contribution to TRIZ in Combining Lean and Inventive Design Method[END_REF], [START_REF] Hani | A Method to Formulate Problem in Initial Analysis of Inventive Design[END_REF] proposed the Inverse Problem Graph (IPG) by applying Lean principles [6]. The feature that differentiates IPG from other existing methods is the way used to formulate a complex problem. Tools such as Network of Problems require looking for the effect of a general problem located at the upper level of a problem situation. In contrast, IPG allows designers to look for the causes of an initial problem at the lower level, helping to improve the agility of the inventive design process. However, IPG does not provide a method for prioritizing the existing initial problems by considering project criteria. Hence, it is necessary to integrate this method with approaches such as FMEA to select the most critical initial problem.

FMEA is a widely applied engineering technique to improve a system by identifying and eliminating potential failures [START_REF] Wang | A risk evaluation and prioritization method for FMEA with prospect theory and Choquet integral[END_REF]. The original version of this technique uses the Risk Priority Number (RPN) method to prioritize the identi ed failures [8]. However, the RPN method has been criticized in the literature due to some drawbacks [START_REF] Liu | Improving Risk Evaluation in FMEA With Cloud Model and Hierarchical TOPSIS Method[END_REF]. Therefore, many authors have proposed to integrate the FMEA approach with alternative methods such as the Analytical Hierarchy Process (AHP) [START_REF] Lo | A novel multiple-criteria decision-making-based FMEA model for risk assessment[END_REF].

The Analytic Hierarchy Process is a mathematical method for prioritizing a given set of alternatives by considering the decision criteria [START_REF] Kurttila | Utilizing the analytic hierarchy process (AHP) in SWOT analysis-a hybrid method and its application to a forest-certi cation case[END_REF]. Since AHP's introduction, this method has been applied to many types of decision problems, such as the justi cation of new manufacturing technology or the evaluation of strategic alternatives [START_REF] Jiaqin | An AHP decision model for facility location selection[END_REF]. This paper aims to integrate the FMEA and AHP methods into the Inverse Problem Graph process. In this integration, FMEA helps to identify the initial problems, and AHP contributes to prioritizing them. The main contribution of the paper is to introduce a methodology for formulating the most crucial contradiction in each iteration of the initial analysis of inventive design.

The rest of the paper is organized as follows: In the rst part, we present the literature relating to the IDM framework, a description of the FMEA method and its integration to improve TRIZ, and a review of the AHP method. Then, the second part shows the structure of the proposal and its application in the Lattice Structure case study. In the last part, we present a comparison between our proposal and the classical method, a discussion, and the conclusion of the paper.

State Of The Art

Inventive Design Methodology (IDM) framework

Due to the development of IPG in the framework of Inventive Design Methodology, we need to provide a brief description of this framework. IDM is a systematic TRIZ-based approach introduced to complement the TRIZ's body of knowledge with other theories such as graph theory [START_REF] Cavallucci | From TRIZ to Inventive Design Method (IDM): towards a formalization of Inventive Practices in R&D Departments[END_REF]. The IDM's framework consists of the following phases [START_REF] Zanni-Merk | Use of formal ontologies as a foundation for inventive design studies[END_REF]:

1. Initial Analysis Phase: This phase begins with the collection of all the knowledge related to the problem situation. Then, this collected knowledge is converted into a graphical model using tools such as Inverse Problem Graph or Problem Graph [START_REF] Cavallucci | Designing the Inventive Way in the Innovation Era[END_REF].

Contradiction Formulation Phase:

The aim of this phase is to formulate the contradictions used as inputs for the application of TRIZ tools [START_REF] Cavallucci | From TRIZ to Inventive Design Method (IDM): towards a formalization of Inventive Practices in R&D Departments[END_REF], [START_REF] Chibane | Innovation beyond optimization: Application to cutting tool design[END_REF]. to select the most appropriate ones [START_REF] Cavallucci | From TRIZ to Inventive Design Method (IDM): towards a formalization of Inventive Practices in R&D Departments[END_REF].

Solution Concept Synthesis

In this paper, we will focus on the initial analysis phase (phase 1) of this TRIZ-based framework. In the next section, we will provide a description of the FMEA method. Furthermore, we will review its closest studies to the rst phase of IDM.

FMEA method and its integration to improve TRIZbased processes

Failure Mode and Effect Analysis (FMEA) is a widely applied methodology for identifying and prioritizing critical failures in a system [START_REF] Liu | Improving Risk Evaluation in FMEA With Cloud Model and Hierarchical TOPSIS Method[END_REF], [START_REF] Rhee | Using cost based FMEA to enhance reliability and serviceability[END_REF]. Since the introduction of FMEA in 1960, FMEA has been applied in various elds such as the chemical, aeronautical, automotive, marine and nuclear industries [START_REF] Wang | A risk evaluation and prioritization method for FMEA with prospect theory and Choquet integral[END_REF].

Several authors have proposed to integrate FMEA into TRIZ. Among these integrations, Russo et al. [18] developed a tool by integrating FMEA, TRIZ function analysis, and Anticipatory Failure Determination (AFD) to identify the cause of failure. Besides, Regazzoni et al. integrated FMEA with TRIZ function analysis and Su-Fields models. This integration aimed to reduce failure occurrence by increasing the capability of anticipating problems and technical solutions [START_REF] Regazzoni | TRIZ tools to enhance risk management[END_REF]. Moreover, Mzougui et al. proposed to integrate the advantages of the TRIZ AFD method into FMEA. The objective of this integration was to improve the failure identi cation process [START_REF] Mzougui | Proposition of a modi ed FMEA to improve reliability of product[END_REF]. In addition, Sprea co et al. proposed to enhance the FMEA by integrating functional analysis, Film Maker, perturbed functional analysis, and subversion analysis [START_REF] Sprea Co | Case: Can TRIZ Functional Analysis Improve FMEA?[END_REF]. But closer to our research, Hakim et al. proposed to apply FMEA to prioritize the output contradictions of the Network of Problems (NoP) [START_REF] Hakim | Reduction Defect in Sewing Work Stations by Integrating OTSM-TRIZ and FMEA[END_REF]. This application can reduce the number of contradictions to solve in subsequent phases of inventive design. Furthermore, it can help designers to focus on problems with a higher degree of importance. However, this proposal cannot help to decrease the number of formulated problems and contradictions in the initial analysis phase of the inventive design.

The traditional FMEA approach prioritizes the failures utilizing the Risk Priority Number (RPN), which is calculated by multiplying the scores of risk factors, namely, detection (D), severity (S), and occurrence (O) as shown in Eq. 1 [START_REF] Sawhney | A modi ed FMEA approach to enhance reliability of lean systems[END_REF].

1

The classical RPN formula is straightforward and very simple to understand [START_REF] Ciani | A Critical Comparison of Alternative Risk Priority Numbers in Failure Modes, Effects, and Criticality Analysis[END_REF]. However, the RPN formula only consider three factors to prioritize failures, which may neglect other important factors such as time, cost, urgency, etc. [START_REF] Mzougui | Proposition of a modi ed FMEA to improve reliability of product[END_REF]. Besides, it does not take into consideration the possible weights for the importance of its risk factors [START_REF] Lo | A novel multiple-criteria decision-making-based FMEA model for risk assessment[END_REF]. To overcome the mentioned limitations, several studies in the literature have developed various models [START_REF] Lo | A novel multiple-criteria decision-making-based FMEA model for risk assessment[END_REF]. The analytic hierarchy process (AHP) is one of these models.

Analytical Hierarchy Process (AHP)

AHP is an approach that was developed to support multi-criteria decision making [START_REF] Wang | An integrated AHP-DEA methodology for bridge risk assessment q[END_REF]. Satty introduced this approach to determine a ranking of the choices when all the decision criteria are considered [26]. The process of the AHP method consists of the following steps [START_REF] Dağdeviren | Decision making in equipment selection: an integrated approach with AHP and PROMETHEE[END_REF]:

(1) In the rst step, a hierarchical structure is established [28]. To do so, AHP decomposes a complex decision problem into a hierarchy of interrelated decision alternatives and criteria [START_REF] Dağdeviren | Decision making in equipment selection: an integrated approach with AHP and PROMETHEE[END_REF]. This hierarchy consists of at least three levels: the decision alternatives at the bottom, the criteria in the middle of the structure, and the objective located at the top level [START_REF] Wang | An integrated AHP-DEA methodology for bridge risk assessment q[END_REF].

(2) The second step is to construct the pairwise comparison matrix to indicate the importance of criteria and alternatives [START_REF] Onut | Transshipment site selection using the AHP and TOPSIS approaches under fuzzy environment[END_REF]. The construction of this matrix is performed according to a numerical rating

RP N = S × O × D
including nine ranking scales, as shown in Table 1 [START_REF] Dağdeviren | Decision making in equipment selection: an integrated approach with AHP and PROMETHEE[END_REF].

(3) In the third step, the pairwise matrices should be normalized to calculate the priority weight of alternatives and criteria [START_REF] Onut | Transshipment site selection using the AHP and TOPSIS approaches under fuzzy environment[END_REF].

(4) The fourth step of AHP is de ned to check the consistency of the pairwise matrices output [START_REF] Dağdeviren | Decision making in equipment selection: an integrated approach with AHP and PROMETHEE[END_REF]. To do this, the consistency index (CI) is rst determined by applying the following equation [START_REF] Onut | Transshipment site selection using the AHP and TOPSIS approaches under fuzzy environment[END_REF]:

2
Where λ max is the principal Eigenvalue, and n is the number of compared elements in the comparison matrix.

Then, the consistency ratio (CR) is obtained as the ratio between the CI and the random index (RI), as shown in Eq. 3 [28]. This ratio is used to check the consistency of the comparison matrix.

3

In this equation, the value of RI relates to the dimension of the matrix, as shown in Table 1. The acceptable value of CR must be less than 0.1. Otherwise, the matrix construction process and its evaluation must be repeated [START_REF] Lin | Using AHP and TOPSIS approaches in customerdriven product design process[END_REF].

Proposed Method: Integration Of (Fmea-ahp) Into The Inverse Problem Graph

This proposed methodology consists of three phases to show the most important and critical contradiction in the rst iteration of the Inverse Problem Graph. Figure 1 illustrates these phases. In the rst phase, the application of the FMEA method helps to determine the objective, possible initial problems, and the project criteria. The second phase relates to prioritizing the initial problems by applying

C I = λ max -n n -1 CR = CI RI
the AHP method. Then, in the third phase, the inverse problem graph of the selected initial problem should be constructed.

In the rst phase, the initial problems and the criteria used to evaluate them are determined. This phase begins with the de nition of the project objective. Next, the initial problems are identi ed as the failure modes in the FMEA. In the nal step of this phase, the criteria are determined as the risk factors of the FMEA method.

The second phase is de ned to prioritize the initial problems by applying the AHP method. This phase consists of the following steps:

Step 1: Develop a hierarchical structure for the objective, initial problems, and criteria: In this step, the hierarchical structure related to the determined criteria, the identi ed initial problems, and the de ned objective should be developed, as shown in Fig. 2.

Step 2: Create pairwise comparison matrices for the initial problems: The second step of this phase relates to creating the pairwise comparison matrices of the initial problems by considering each criterion.

In this step, a numerical rating with nine ranking scales can be used, as shown in Table 2. Let A in Eq. 4 represent n x n pairwise contradiction matrix related to the initial problems.

A = (a ij ) nxn = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 a 12 … a 1n a 21 = 1 … a 2n ⋮ ⋮ ⋮ ⋮ a n1 = a n2 = … 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1 a 12 1 a 1n 1 a 2n
Step 3: Calculate relative priorities for the initial problem matrices: In this step, the priorities related to each initial problem in the matrices must be calculated. For this purpose, each set of column values is summed. Subsequently, each value in the matrices is divided by the sum of its column. Finally, the average of row values is calculated to obtain the priorities of the matrices.

Step 4: Calculate the consistency of initial problem matrices: It is necessary to check the consistency of the pairwise matrices related to the initial problems in the fourth step. This step includes the following sub steps: If the consistency ratio is less than the threshold of 0.1, the pairwise comparison matrix has an acceptable consistency; otherwise, the process of constructing the matrix must be repeated.

Step 5: Create a comparison matrix for the criteria: This step relates to creating the comparison matrix for the determined criteria in the rst phase with respect to the objective of the project. Let B, as shown in Eq. 7, represent m x m pairwise contradiction matrix related to the criteria.

7

Step 6: Calculate relative priorities for the criteria matrix: In this step, we calculate the priorities related to each criterion in the matrix, as we should accomplish in Step 3 of the process for the initial problem

C I I P = λ max -n n -1 C R I P = C I I P RI B = (b ij ) mxm = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 b 12 … b 1m b 21 = 1 … b 2m ⋮ ⋮ ⋮ ⋮ b m1 = b m2 = … 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1 b 12 1 b 1m 1 b 2m
matrices.

Step 7: Calculate the consistency of the criteria matrix: It is necessary to check the consistency of the pairwise matrix related to the criteria in step seven by applying the three following sub-steps: Here also, if the consistency ratio is less than the threshold of 0.1, the pairwise comparison matrix has an acceptable consistency; otherwise, the process of constructing the matrix should be repeated.

Step 8: Calculate the priority of each initial problem: In this step, the priorities of the initial problems with respect to the objective of the project should be calculated.

Step 9: Select the most important initial problem of the project: The last step of the second phase is to select the most important initial problem taking into account the calculated priorities.

After prioritizing the initial problems and selecting the most important one, it is time to identify the contradictions in the third phase through the application of the IPG method presented by [START_REF] Hani | Contribution to TRIZ in Combining Lean and Inventive Design Method[END_REF]. The capability of our proposal will be illustrated in the next section using a Lattice Structure case study.

Application Of

The Proposed Method To The Case Study "lattice Structure"

In this section, we applied the proposed method "FMEA-AHP based method" to the Lattice Structure (LS) case study to evaluate its applicability. Lattice Structures (LS) are open cell structures constructed by one or more reproducible unit cells [START_REF] Maconachie | SLM lattice structures: Properties, performance, applications and challenges[END_REF]. These structures are types of cellular solids that are classi ed into two-and three-dimensional structures [START_REF] Tao | Design of lattice structure for additive manufacturing[END_REF]. The superior properties of LSs make them applicable in various engineering elds such as heat transfer, energy absorption, and biological bone craft [START_REF] Ma | Lattice structures of Cu-Cr-Zr copper alloy by selective laser melting: Microstructures, mechanical properties and energy absorption[END_REF]. In this

C I C = λ max -n n -1 C R C = C I C RI
case study, we formulate the contradiction related to the lattice structure by considering the priority of the initial problems, and we propose a solution to solve the contradiction. To do this, we rst prioritize the initial problems through the rst and second phases of the proposal. We then apply the Inverse Problem Graph to construct the contradiction. At the end of the case study, we apply the TRIZ tools to propose the solution. In what follows, we explain the application of our proposal to this case study in the following phases.

In the rst phase of the process, the designer should de ne the objective at the beginning of the project. In our case study, the "Improvement of Lattice Structure" was selected as the objective. Then, we analyze the situation to identify the problems that prevent the objective from being achieved. Therefore, we determined "IP 1 : The structure does not have a good strength.", "IP 2 : The energy absorption of the structure has been reduced.", and "IP 3 : The structure is heavy." as the existing initial problems in this case study. At the end of this phase, we determined three criteria "C 1 : Strategic", "C 2 : Financial", and "C 3 :

Urgency".

In second phase, designer developed the hierarchical structure of the objective, initial problems, and criteria. As shown in Fig. 3, this structure has three levels for prioritizing the problems. The identi ed project objective is located at the rst level. The criteria and the problems are located on the second and third levels of hierarchy, respectively.

After forming the hierarchical structure of the lattice structure, the initial problems were prioritized by considering the criteria through the steps in the second phase. To do this, we rst constructed the pairwise matrices based on the project criteria. Then, the priority related to each initial problem according to the criterion of its matrix was calculated. In the following, we checked the consistency of the pairwise matrices. For instance, as shown in Table 3, we rst created a pairwise matrix (reciprocal matrix) for the initial problems de ned with respect to the criterion "C 1 : Strategic". To construct this matrix, we applied the scale given in Table 2. Once the pairwise matrix of the initial problems based on the strategic criterion was created, we converted the fractional values into decimal values. Next, we normalized the matrix. To this end, we divided each element of the matrix by the sum of its column. Table 4 shows the normalized matrix of the initial problems related to the rst criterion. After constructing the normalized matrix, we calculated the priorities related to the initial problems with respect to the rst criterion. To do this, we calculated the average value of all the elements existing in each row of the matrix, as illustrated in Table 5.

Subsequently, we needed to calculate the consistency of the matrix to check whether the calculated values were correct or not. To do this, we calculated the principal Eigen value (λ max = 3.11) related to the created matrix. Then, the Consistency Index ( ) and Random Index (RI=0.58, Table 2) were obtained. Finally, we calculated the Consistency Ratio ( ) of the matrix, which was less than 0.1. Hence, we accepted the created pairwise matrix.

RI

Similarly, we constructed the pairwise matrices of the initial problems related to the other de ned criteria de ned in the "Lattice Structure" case study. Tables 6 and7 present these matrices. After constructing the pairwise matrix for the initial problems, it was also essential to evaluate criteria of the project by forming their corresponding matrix. Therefore, we rst assigned weights to the criteria by applying the scale given in Table 1. Then, we calculated the relative priorities and the consistency of the matrix. Table 8 shows the pairwise matrix related to the criteria. At the end of the second phase of the process, the overall priority of each initial problem with respect to the project objective was calculated. To do this, we summed the products of the criterion priority by the priority of its initial problem. For example, the overall priority of the rst initial problem was calculated:

. Figure 4 shows a graphical representation of the priorities of the criteria and initial problems. Furthermore, Table 9 demonstrates the overall priorities related to the initial problems of the case study with respect to the objective of the project. As a result, we selected "IP 2: The energy absorption of the structure has been reduced." as the most important initial problem. In the third phase, we formulated the contradiction related to the selected intital problem by using IPG. To do this, we rst identi ed "The structure has a high volume fraction." and "The lattice structure has solid truss lattice" as the causes that reduce the energy absorption of the structure. Next, we ranked the identi ed causes, and we selected "The structure has a high volume fraction" as the most important cause that has an impact on the initial problem. In what follows, we determined the type of the selected cause, which was a Harmful-Useful problem. Therefore, we converted it to the partial solution "To increase the volume fraction of the structure". This means that we converted its structure into the structure of the partial solution. Then, we identi ed the problem, "The structure does not have good strength" as the cause of the partial solution. Subsequently, we constructed the IPG of the rst iteration, as shown in Fig. 5.

In the following, we extracted the illustrated contradiction from the graph, as shown in Fig. 6. In this gure, there is a contradiction between energy absorption of the structure and its strength. This means that attempting to the increase energy absorption leads to a decrease in the strength of the structure. Subsequently, we should allocate the parameters to the problems and the partial solution of the extracted contradiction. As shown in Fig. 8, we assigned "Reliability of the structure" as the evaluation parameter to the initial problem "The energy absorption of the structure has been reduced". The second evaluation parameter "Strength of the structure" was allocated to the cause of the partial solution "The structure does not have good strength". We assigned also "Increase in volume fraction" as an action parameter to the partial solution. At the end of this phase, we constructed the poly-contradiction table by applying the assigned parameters in the previous step to formulate the contradiction, as illustrated in Fig. 7. The contradiction states that the structure has to be fabricated at a high-volume fraction to be rigid. In the same system, it should be fabricated at a low-volume fraction to increase the mechanical energy absorption capacity.

Increase of volume fraction

No Yes

Strength of the structure ϑ Λ

Reliability of the structure (energy absorption) Λ ϑ

Reliability of the structure: The capability of object to perform its function in a predictable manner.

Strength of the structure: The ability of the structure to resist changes in response to force.

Figure 7. The poly-contradiction of the Lattice Structure case study.

In the last phase of the case study, we used the separation principles [START_REF] Savransky | Engineering of Creativity: Introduction to TRIZ Methodology of Inventive Problem Solving[END_REF], introduced in the frame of TRIZ theory, to solve the formulated contradiction. These principles are: 1) Separation in time, 2) Separation based on a condition, 3) Separation in space, and 4) Separation between the parts and the whole. Among them, we selected the principle "Separation between the parts and the whole".

To implement this principle, we used the method "Merging-Interacting". This method states that when several individual elements have property A, they can be uni ed through a eld, mediator, or uni ed arrangement to have property B. Based on this method, we proposed to merge small cells to form a larger cell of lattice structure with the same geometry, as shown in Fig. 8. In this solution, our proposal has both high and low volume fractions. At the macro level, it presents a relatively low volume fraction, which contributes to high energy absorption. On the other hand, it has a high-volume fraction at the micro level, which allows us to construct a structure with good stiffness.

To evaluate the performance of the proposed solution, we carried out a numerical simulation using Abaqus FEA Software. For this purpose, we built three models, including Solid structure, Conventional lattice structure, and our solution. In this section, we compared Network of Problem (NoP) and our FMEA -AHP based method. This comparison allows the readers of this research to be more aware of their differences. To make this comparison, we extracted information from the article [START_REF] Khomenko | A framework for OTSM TRIZ-based computer support to be used in complex problem management[END_REF], in which the authors applied NoP to the "Biomass Power Plant" case study to present its network problem situation. As Table 10 and Fig. 9 (a)

illustrate, the problem "Biomass power plant should be improved" was initially selected as the general one, located at the highest level of the network. Subsequently, designers constructed the chain of problems and partial solutions to explore several problems at the lowest level of the network. Figure 9 (b) shows the chain related to one of these identi ed problems (Problem 68). As Fig. 10 illustrates, we applied FMEA-AHP based method to the "Biomass Power Plant" case study. To perform this application, we rst identi ed the initial problems related to the biomass power plant. Then, our methodology is used to prioritize them according to the determined criteria for the case study. "High ue gas speed and temperature destroys the post combustion chamber at the turn point" was selected as the most important initial problem. At the end, we determined the causes of the initial problem and formulated the contradiction by applying the IPG method.

Discussion

In this study, we proposed a new method, integrating FMEA and AHP into the Inverse Problem Graph, to prioritize the initial problems in the initial analysis phase of inventive design. The contribution of this proposal to the inventive design process is re ected in several aspects. First, this method contributes to prioritizing initial problems at the beginning of the project, thus eliminating time spent on less important problems from the process. As the second contribution, it allows designers to formulate contradictions based on the priorities of the initial problems. This saves designers time in resolving the most appropriate contradiction. Thirdly, we can mention that this method allows the designers to explore the root contradictions by integrating Inverse Problem Graph's steps into its process.

The analysis of our case study also displays several limitations. As the rst one, we can mention that the mathematical calculations of such a method are very time-consuming. This drawback can signi cantly decrease the agility of the inventive design process. As the second one, the current process does not provide any technique to identify initial problems related to a complex system.

Conclusion

In this paper, we introduced a new method for prioritizing initial problems in the initial analysis of the inventive design process. To develop our method, we rst reviewed some documents on the FMEA approach to highlight its drawbacks and advantages. Our reviews revealed that this approach uses RPN, considering only three factors to prioritize failures. To overcome this drawback of FMEA, several methods are presented in the literature. Among them, there is AHP that is developed to support multi-criteria decision making. Next, we proposed to integrate FMEA, AHP methods into the Inverse Problem Graph to prioritize the initial problems in the initial analysis phase of the inventive design process. Finally, we tested the capability of our approach by applying it to the Lattice structure case study. Based on this application, we realized that our proposal can help formulate the most important problem in each iteration of the process. It can also improve the time allocated to the initial analysis.

Further research is necessary to appreciate our proposal. Primarily, we need to develop software to reduce the computational procedure during the problem formulation in the initial analysis phase. Secondly, it is necessary to integrate the proposal process with other existing methods to improve its agility. And third, there is an opportunity to use approaches that help designers to identify the initial problems in complex systems.

Figures

Figure 1 The different phases of FMEA-AHP based method

A three-level hierarchical structure for the objective, criterion, and the initial problems

The hierarchical structure of the Lattice Structure case study

Model of the proposed lattice structure solution with a high-volume fraction at micro level, but a low volume fraction at macro level

  Phase: Designers could apply TRIZ methods to solve technical or physical contradictions in this phase [16]. 4. Solution Concept Selection Phase: In the last phase, external experts evaluate the solution concepts

Sub-step 1 : 5 Sub-step 3 :

 153 Calculate the principal Eigen value (λ max ): The principal Eigen value should be calculated. Sub-step 2: Calculate the consistency index (CI IP ): The consistency index is determined using the following equation: Calculate the consistency ratio (CR IP ): The consistency ratio is calculated as the ratio of CI IP to the random index (RI), as shown in the following equation: 6

Sub-step 1 :

 1 Calculate the principal Eigen value (λ max ): The principal Eigen value should be calculated. Sub-step 2: Calculate the consistency index (CI C ): The consistency index is determined through the following equation: 8 Sub-step 3: Calculate the consistency ratio (CR C ): The consistency ratio is calculated as the ratio of CI C to the random index (RI), as shown in the following equation: 9

Figure 4 .

 4 Figure 4. Graphical representation of the priorities related to the problems and criteria

(( 0 Figure 6 .

 06 Figure 6. The allocation of parameters to the extracted contradiction

Table 1 The

 1

					value of the random index
	n	1	2	3	4	5	6	7	8	9	10
	R.I. 0	0	0.58	0.90	1.12	1.24	1.32	1.41	1.45	1.49

AHP can take into account a large number of qualitative and quantitative criteria to prioritize the failure modes

[26]

. Furthermore, it has the ability to address the relative importance of the risk factors

[8]

.

Table 2

 2 The 1-9 scales for pairwise comparison in AHP

	Importance intensity	De nition
	1	Equal importance
	3	Moderate importance of one over another
	5	Strong importance of one over another
	7	Very strong importance of one over another
	9	Extreme importance of one over another
	2,4,6,8	Intermediate values
	Reciprocals	Reciprocals for inverse comparison
	4

Table 3

 3 Reciprocal matrix for the de ned initial problems with respect to the criterion "Strategic"

	C 1 : Strategic	IP 1 : Structure does	IP 2 : The energy absorption	IP 3 : The
		not have a good	of the structure has been	structure is
		strength.	reduced.	heavy.
	IP 1 : Structure does not have	1	1 3	5
	a good strength.
	IP 2 : The energy absorption	3	1	7
	of the structure has been
	reduced.
	IP 3 : The structure is heavy.	1 5	1 7	1

Table 4

 4 Normalized matrix of the initial problems related to the criterion "Strategic"

	C 1 : Strategic	IP 1 : Structure does	IP 2 : The energy absorption	IP 3 : The
		not have a good	of the structure has been	structure is
		strength.	reduced.	heavy.
	IP 1 : Structure does not have	0.238	0.224	0.384
	a good strength.
	IP 2 : The energy absorption	0.714	0.680	0.538
	of the structure has been
	reduced.
	IP 3 : The structure is heavy.	0.047	0.095	0.076

Table 5

 5

	shows the nal matrix of the initial problems related to the rst
			criterion.
	Strategic C 1	IP 1	IP 2	Strategic C 1	IP 1
	IP 1	0.238	0.224	IP 1	0.238
	IP 2	0.714	0.680	IP 2	0.714
	IP 3	0.047	0.095	IP 3	0.047
	Consistency	λ

max = 3.11, CI = 0.055, CR = 0.095 < 0.1 (Acceptable)

Table 5 :

 5 Pairwise matrix for the initial problems considering the C 1 criterion.

		λ max -3
	CI IP =	3-1	= 0.055
			CI IP
			CR IP =	= 0.095

Table 6

 6 Pairwise matrix for the initial problems considering the C 2 criterion.

	Financial C 2	IP 1	IP 2	IP 3	Priority
	IP 1	0.19	0.183	0.285	0.219
	IP 2	0.761	0.734	0.642	0.713
	IP 3	0.047	0.081	0.071	0.066
	Consistency	λ max = 3.06, CI = 0.03, CR = 0.053 < 0.1 (Acceptable)

Table 7

 7 Pairwise matrix for the initial problems considering the C 3 criterion.

	Urgency C 3	IP 1	IP 2	IP 3	Priority
	IP 1	0.157	0.272	0.148	0.192
	IP 2	0.052	0.090	0.106	0.080
	IP 3	0.789	0.636	0.744	0.721
	Consistency	λ max = 3.11, CI = 0.055, CR = 0.095 < 0.1 (Acceptable)

Table 8

 8 Pairwise matrix for the de ned case study criteria.

	Criteria	C 1	C 2	C 3	Priority
	C 1	0.652	0.692	0.555	0.633
	C 2	0.217	0.230	0.333	0.260
	C 3	0.130	0.076	0.111	0.106
	Consistency	λ

max = 3.06, CI = 0.034, CR = 0.059 < 0.1 (Acceptable)

Table 9

 9 Prioritization of the initial problems with respect to the objective

	Initial problems C 1	C 2	C	Overall priority	Priority
	IP 1	0.178	0.056	0.020	0.255	2
	IP 2	0.407	0.185	0.008	0.600	1
	IP 3	0.046	0.017	0.076	0.139	3

  Table 10 shows a comparison between the output results of the simulation tests. In the table, we also calculated the ratios of energy absorption and young's modulus to the target values considered for each model. These ratios for our solution, conventional lattice structure,

	and solid structure are respectively (ratio of energy absorption to considered target value =
	0.017 0.015	= 1.196 > 1	, ratio of young's modulus to considered target value =	0.509 0.4	= 1.274 > 1	), (0.593,
	17.548), and (0.264, 975.556). As the table illustrates, the only model that can satisfy both target values
	is our solution with ratios greater than one.

Table 10

 10 Comparison of our solution with a solid and normal lattice structure.

	Solid	Conventional lattice	Proposed
	structure	structure	solution

Table 10

 10 List of partial solutions and problems related to the "Biomass Power Plant" topic

	PS	escription of Partial Solution (PS)	PS or	escription of Partial Solution (PS) or
	or	or Problem (PB)	PB N°D Problem (PB)
	PB
	N°D
	1.Pb	Power plant should be improved.	13.Ps	Heat exchanger is used.
	…	…	…
	9.Pb	Biogas should be cooled before	66.Pb	Special system is needed to burn solid
		cleaning by modern technological		particles.
		devices.
	…	…	67.Ps	Special form of this system is needed in order
				to change direction of ue gas.
	…	…	68.Pb	High speed, high temperature of ue gases
				destroy the post combustion chamber at the
				turn point.

The IPG related to the selected initial problem Figure 6 The allocation of parameters to the extracted contradiction Figure 7 The poly-contradiction the Lattice Structure case study.

Network of Problems and its difference with the proposal [START_REF] Khomenko | A framework for OTSM TRIZ-based computer support to be used in complex problem management[END_REF] Figure 10 Application of FMEA-AHP based method to the "Biomass Power Plant" case study


N Börlin

A Murtiyoso
email: arnadi.murtiyoso@insa-strasbourg.fr

P Grussenmeyer
email: pierre.grussenmeyer@insa-strasbourg.fr

EFFICIENT COMPUTATION OF POSTERIOR COVARIANCE IN BUNDLE ADJUSTMENT IN DBAT FOR PROJECTS WITH LARGE NUMBER OF OBJECT POINTS

Keywords: Bundle adjustment, Quality control, Posterior covariance, Software, Photogrammetry

One of the major quality control parameters in bundle adjustment are the posterior estimates of the covariance of the estimated parameters. Posterior covariance computations have been part of the open source Damped Bundle Adjustment Toolbox in Matlab (DBAT) since its first public release. However, for large projects, the computation of especially the posterior covariances of object points have been time consuming. The non-zero structure of the normal matrix depends on the ordering of the parameters to be estimated. For some algorithms, the ordering of the parameters highly affect the computational effort needed to compute the results. If the parameters are ordered to have the object points first, the non-zero structure of the normal matrix forms an arrowhead. In this paper, the legacy DBAT posterior computation algorithm was compared to three other algorithms: The Classic algorithm based on the reduced normal equation, the Sparse Inverse algorithm by Takahashi, and the novel Inverse Cholesky algorithm. The Inverse Cholesky algorithm computes the explicit inverse of the Cholesky factor of the normal matrix in arrowhead ordering. The algorithms were applied to normal matrices of ten data sets of different types and sizes. The project sizes ranged from 21 images and 100 object points to over 900 images and 400,000 object points. Both self-calibration and non-self-calibration cases were investigated. The results suggest that the Inverse Cholesky algorithm is the fastest for projects up to about 300 images. For larger projects, the Classic algorithm is faster. Compared to the legacy DBAT implementation, the Inverse Cholesky algorithm provides a performance increase by one to two orders of magnitude. The largest data set was processed in about three minutes on a five year old workstation. The legacy and Inverse Cholesky algorithms were implemented in Matlab. The Classic and Sparse Inverse algorithms included code written in C. For a general toolbox as DBAT, a pure Matlab implementation is advantageous, as it removes any dependencies on, e.g., compilers. However, for a specific lab with mostly large projects, compiling and using the classic algorithm will most likely give the best performance. Nevertheless, the Inverse Cholesky algorithm is a significant addition to DBAT as it enables a relatively rapid computation of more statistical metrics, further reinforcing its application for reprocessing bundle adjustment results of black-box solutions.

INTRODUCTION 1.1 Background

One of the major quality control parameters in bundle adjustment are the posterior estimates of the covariance of the estimated parameters. In photogrammetric projects oriented towards robust measuring purposes, this is not only essential but very important for quality control. Indeed, the covariance values indicate the quality of the bundle adjustment result and thus its feasibility for use as topographic and surveying products. A common practice would be to compare these posterior covariance values with the project requirements, which often times depends on the required map scale. A photogrammetric mission may be considered acceptable if and only if its quality parameters, including the covariance values, are within the project's tolerance.

A more in-depth analysis to the geometrical qualities of a photogrammetric project is also sometimes necessary, for example when errors or deviation from prior values are detected. This becomes more and more important with the increasing * Corresponding author use of UAV (Unmanned Aerial Vehicle) platforms, colloquially known as drones. Contrary to classical metric photogrammetric survey, the use of consumer-grade cameras in drone imagery is less reliable and therefore requires more attention as far as the quality is concerned. This in-depth analysis may be performed by various means, including by computing the correlation rate between the different calculated parameters; a metric which is also derived from the covariance matrix.

Since geometric quality is of the utmost importance in these spatial applications, the computation of covariance is very interesting for quality control purposes. However, for large projects, the computation of especially the posterior covariances of object points (OP) can be a time consuming process.

The Damped Bundle Adjustment Toolbox (DBAT) for Matlab was conceived as an open source toolbox to perform bundle adjustment computations [START_REF] References Börlin | Bundle Adjustment with and without Damping[END_REF]. It has been used for several applications, most notably to reprocess the bundle adjustment results of commercial solutions [START_REF] Murtiyoso | Open Source and Independent Methods for Bundle Adjustment Assessment in Close-Range UAV Photogrammetry[END_REF]. The main highlight of DBAT has always been its modularity [START_REF] Börlin | Flexible Photogrammetric Computations using Modular Bundle Adjustment[END_REF] and openness, with the possibility to access bundle adjustment metrics contrary to the black-box nature of several commercial software.

Related work

Various schemes that exploit the sparsity of the normal matrix to speed up the bundle adjustment computations were presented already in the early days of Bundle Adjustment [START_REF] Brown | Inversion of very large matrices encountered in large scale problems of photogrammetry and photographic astrometry[END_REF][START_REF] Brown | The Bundle Adjustment -Progress and prospects[END_REF] and are now part of standard photogrammetric textbooks (e.g. [START_REF] Kraus | Photogrammetry[END_REF][START_REF] Wolf | Elements of Photogrammetry with Application in GIS[END_REF][START_REF] Luhmann | Close-Range Photogrammetry and 3D Imaging[END_REF][START_REF] Förstner | Photogrammetric Computer Vision[END_REF].

In most cases, the presentations are focused on speeding up the computations during the bundle adjustment iterations, where the linearised equation has a single right-hand side. In contrast, the computation of the posterior covariance has a large number of right-hand sides, and said algorithm cannot be immediately applied. Literature that discusses algorithms for posterior covariance computations include [START_REF] Triggs | Bundle adjustment -modern synthesis[END_REF]; [START_REF] Wolf | Elements of Photogrammetry with Application in GIS[END_REF]; [START_REF] Förstner | Photogrammetric Computer Vision[END_REF]; [START_REF] Ila | Fast incremental bundle adjustment with covariance recovery[END_REF]. However, the details of how the sparsity can be exploited are rarely given.

The method by [START_REF] Takahashi | Formation of a Sparse Bus Impedance Matrix and its Application to Short Circuit Study[END_REF] can be used to compute subsets of the inverse of a sparse matrix. The sparse inverse (SI) method was developed for electrical networks, but can be applied to any sparse matrix problem.

Recently, [START_REF] Kallin | Computation of posterior covariances of object points in bundle adjustment[END_REF] presented an algorithm that computed the explicit inverse of the Cholesky factorisation of the normal matrix using a combination of sparse and dense computations. Depending on the density of a matrix, i.e., the fraction of non-zeros, dense computations can have performance advantage over sparse computations due to their more efficient use of the memory hierarchy [START_REF] Duff | An Overview of the Sparse Basic Linear Algebra Subprograms: The New Standard from the BLAS Technical Forum[END_REF]; [START_REF] Goto | High-Performance Implementation of the Level-3 BLAS[END_REF]. Furthermore, the dense algorithms are easier to parallellise due to their predictable memory access patterns. The inverse Cholesky algorithm (called VDSIBlock in [START_REF] Kallin | Computation of posterior covariances of object points in bundle adjustment[END_REF]) was shown to outperform the SI algorithm on small-tomedium-sized data sets.

Aim

The aim of this paper is to extend the investigation of the Inverse Cholesky algorithm to larger data sets and to compare the results with the current DBAT implementation and classic photogrammetric algorithms.

THE STRUCTURE OF THE NORMAL MATRIX

Non-zero pattern

The normal matrix has a characteristic pattern of non-zero elements. If the parameters are ordered with the external orientation (EO) parameters first, the normal matrix may be blocked as

N =
Nee Neo Noe Noo .

(1)

The internal block structure is visualised in Figure 1. The Nee and Noo blocks are block diagonal with 6-by-6 and 3-by-3 blocks, respectively. The non-zero coefficients of each block relate only to the EO or OP parameters of a single camera or point, respectively. In contrast, the Neo block contain 3-by-6 blocks that connect images to object points.

N oo N oo N oo N ee N eo N oe

Arrowhead permutation

If the ordering of the unknown parameters is changed to have the object points first, the blocks of the normal matrix are swapped to

N = Noo Noe Neo Nee . ( 2
)
With this ordering, the non-zero structure of the matrix has the shape of an arrowhead, where the large Noo block form the shaft of the arrow, and the remaining block form the "tip" of the arrow at the lower and right borders (Figure 2).

Self-calibration

For a self-calibration project, the normal matrix also include blocks for the internal orientation (IO) parameters. With the IO parameters first, the blocking becomes

N =   Nii Nie Nio Nei Nee Neo Noi Noe Noo   , (3)
where the new blocks have the width corresponding to the number of IO parameters and are typically fully dense. The corresponding arrowhead ordering is

N =   Noo Noe Noi Neo Nee Nei Nio Nie Nii   . (4)
See Figure 3.

N ee N ee N oo N oo N oe N eo Nii Nei Noi Nie Nio N oo N oo N ee N ee N eo N oe Nii Nio Nie Noi Nei
Figure 3: A self-calibration normal matrix has an added border of dense blocks. The blocks have a thickness equal to the number of estimated IO parameters. In the standard ordering (left), the IO blocks form the left and upper border of the matrix. In the arrowhead ordering, the IO blocks instead form the lower and right border.

Algorithm 1 The classic algorithm for computing posterior covariance for object points. The sparse and dense steps are computed in Matlab using sparse and dense linear algebra routines, respectively. The C code step is computed by compiled C-code called from Matlab ("mex" code).

Require: N has standard ordering 1: A ← Nee.

Extract blocks

2: B ← Neo. 3: C ← Noo. 4: E ← D -1 Sparse 5: F ← BEB T Sparse 6: Cee ← F -1 Dense 7: G ← BE Sparse 8: H ← G T CeeG C code, 3-by-3 diagonal blocks only 9: Coo ← G + H Sparse 3. ALGORITHMS

Classic

The classic algorithm, as described in Wolf et al. (2013, Ch. 17-12) uses the normal matrix N in standard ordering (eq. ( 1)):

N = Nee Neo Noe Noo = A B B T D , (5)
where m and n are the number of images and object points, respectively.

Given the same partitioning of the posterior covariance

N -1 = C = Cee Ceo Coe Coo , (6)
the matrix Coo contain the posterior covariance for the object points. Mathematically, the matrix Coo can be computed as follows (cf. (Wolf et al., 2013, eq. (17-34-36))):

Cee = (A -BD -1 B T ) -1 , (7) Coo = D -1 + D -1 B T CeeBD -1 . ( 8
)
The actual computation is given in Algorithm 1.

Inverse Cholesky

The Inverse Cholesky algorithm uses the arrowhead ordering (eq. ( 2))

N = Noo N T eo Neo Nee , (9)
L 11 L 11

L 22 L 22 L 21 0 K 11 K 11 K 22 K 22 K 21 0 Figure 4:
The non-zero patterns of the Cholesky factor L (left) forms half an arrowhead. The L11 block is block diagonal with 3-by-3 lower triangular blocks. Due to the arrowhead shape of the normal matrix, the fill-in during the Cholesky factorisation is restricted to the lower triangular L22 block. The L21 block does not experience any fill-in at all. The non-zero pattern of the Cholesky inverse K = L -1 (right) is similar to that of L, except that the K21 and K22 block may experience fill-in. The K22 block will typically be close to half full.

The normal matrix is Cholesky factorised into

LL T = N with L = L11 0 L21 L22 . ( 10
)
The non-zero structure of L forms half an arrowhead (Figure 4, left). The L11 block will be 3-by-3 block diagonal with lower triangular blocks. Due to the arrowhead structure of the normal matrix, there will be no fill-in in the L21 block during the Cholesky factorisation. Instead, the fill-in is restricted to the lower triangular L22 block.

The inverse K of the Cholesky factor is

K = L -1 = K11 0 K21 K22 . ( 11
)
This leads to the following equations:

K11 = L -1 11 , ( 12
) K21 = -L -1 22 L21K11, ( 13
) K22 = L -1 22 . ( 14
)
The K11 block can be computed using sparse linear algebra and will have the same sparsity pattern as L11. In the computation of K21, K11 is known to be sparse, whereas the density of L21 and L22 may vary. The K21 block can expect a fill-in compared to L22. The K22 block will typically be close to half full. See Figure 4, right.

Using the Cholesky factor, the posterior covariance matrix can be computed as

C = N -1 = (LL T ) -1 = K T K (15a) = K T 11 K11 + K T 21 K21 K T 21 K22 K T 22 K21 K T 22 K22 . ( 15b
)
For notational convenience, let P = K11, Q = K21. Then,

Coo = P T P + Q T Q. ( 16
)
Note that the K22 block is not needed. The structure of the 3-by-3 diagonal blocks of Coo are

Coo =            c11 c12 c13 c21 c22 c23 c31 c32 c33 c44 c45 c46 c54 c55 c56 c64 c65 c66 . . .            =            k T 1 k1 k T 1 k2 k T 1 k3 k T 2 k1 k T 2 k2 k T 2 k3 k T 3 k1 k T 3 k2 k T 3 k3 k T 4 k4 k T 4 k5 k T 4 k6 k T 5 k4 k T 5 k5 k T 5 k6 k T 6 k4 k T 6 k5 k T 6 k6 . . .            , ( 17
)
where ki are column vectors of K. Thus, the diagonal elements of Coo are the inner products of ki with themselves. Furthermore, the (2,1) element within each block is the inner product of ki and ki+1, and similar for the other off-diagonal elements. Due to symmetry, only the sub-diagonal elements need to be computed.

Given that

ki = pi qi , (18)
the inner product can be computed as

k T i kj = p T i pj + q T i qj, ( 19
)
where pi is sparse with at most three non-zeros, but qi may be dense. Thus, the computation may be split into a combination of sparse and dense linear algebra.

The computation of the diagonal elements can be performed efficiently in Matlab using a combination of the Hadamard (element-wise) product and the colsum function. The call v=colsum(A) returns a vector such that vi holds the sum of the elements of the i:th column of A. Also, let the function extract(A,d) return every third column of A starting with column d.

The actual computation is given in Algorithm 2. The diagonal elements are collected in one 3n-vector d, where n is the number of object points. The sub-diagonal elements are collected in three n-vectors sij. Finally, the call C=assemble(d, s21, s31, s32) assembles the vectors to a blockdiagonal matrix C.

Sparse inverse

The SparseInv (SI) algorithm [START_REF] Förstner | Photogrammetric Computer Vision[END_REF][START_REF] Kallin | Computation of posterior covariances of object points in bundle adjustment[END_REF] is a C implementation of the [START_REF] Takahashi | Formation of a Sparse Bus Impedance Matrix and its Application to Short Circuit Study[END_REF] algorithm [START_REF] Davis | Sparseinv: Sparse inverse subset[END_REF]. The SI algorithm computes a subset of the non-zero elements of the inverse C = N -1 .

Legacy DBAT

The posterior covariance computation in DBAT up until version 0.9.1 used the inverse Cholesky approach but looped over the object points.

If the normal matrix is in arrowhead ordering and the identity matrix is column blocked as

I = E1 E2 ... , (20
W ← L -1 E Sparse triangular solve 5: U ← W T W Sparse 6:
Make U block diagonal.

7:

Insert U in C at the appropriate place. 8: end for where each block column is three columns wide, the block Cii corresponding to object point i becomes (cf. eq. ( 15a))

Cii = E T i CEi = E T i N -1 Ei = E T i K T KEi = W T i Wi, (21)
where Wi = KEi.

In an attempt to use memory more efficiently, the implemented algorithm loops over blocks of 30 points. The algorithm is given as Algorithm 3.

Total variance only

The diagonal elements of Coo can be used to compute the total variance (trace) of each point. The necessary modification of the classic algorithm is to modify step 8 to compute the diagonal elements only. In the inverse Cholesky algorithm, steps 8-15 are omitted and the matrix is assembled from the diagonal only.

Self-calibration

For self-calibration projects, the only modification of the classic Algorithm 1 is to have steps 1-2 include the IO blocks, i.e., to match

N =   Nii Nie Nio Nei Nee Neo Noi Noe Noo   = A B B T D . ( 22
)
No modification of the inverse Cholesky algorithm is necessary. Instead the L21 and L22 blocks are re-interpreted to correspond to all non-OP parameters.

Data sets

Ten data sets of varying sizes were prepared for the experiments:

Camcal A 21-image camera calibration data set from [START_REF] Börlin | Camera Calibration using the Damped Bundle Adjustment Toolbox[END_REF].

Vexcel A 41-image vertical aerial dataset flown at an altitude of 950 m above sea level (ASL) above the city of Strasbourg with an UltraCam Osprey Mark 3 Premium camera.

Roma A 60-image terrestrial close-range data set of the Arco di Constantino monument in Rome, Italy from [START_REF] References Börlin | Bundle Adjustment with and without Damping[END_REF].

StPierre A 239-image close-range data set of the St-Pierre-le-Jeune church in Strasbourg from [START_REF] Murtiyoso | Reprocessing Close Range Terrestrial and UAV Photogrammetric Projects with the DBAT Toolbox for Independent Verification and Quality Control[END_REF].

Sewu A 210-image data set acquired by the DJI Phantom 4 RTK drone flown at 60 m ASL over the Sewu temple in Central Java, Indonesia.

UmU Five data sets from several flights with a DJI Phantom 4 V2.0 drone over the Umeå University Campus, Umeå, Sweden.

UmU-3D A 272-image data set flown at 70 m ASL in "3D" cross-hatch mode. The data set includes both east-west and north-south strips. UmU-2D A 133-image subset of UmU-3D with the eastwest strips only. UmU-3D-R The UmU-3D data set, reduced to 210 images to match Sewu. The reduction was done by removing every other north-south strip. UmU-2-lyr UmU-3D extended with another "3D" flight at 120 m ASL over a larger area. UmU-3-lyr UmU-2-lyr extended with another "3D" flight at 40 m ASL over a smaller area.

All data sets were self-calibration data sets. The Camcal and Roma images were processed by PhotoModeler Scanner 20121 and imported to DBAT. The remaining data sets were processed by AgiSoft Metashape v1.62 . During import of the Metashape projects into DBAT, all object points with fewer than 3 rays and/or an estimated intersection angle below 5 degrees were removed before processing by DBAT. The statistics for the data sets are presented in Table 1.

Algorithm variants

The following algorithm variants were implemented and used to process normal matrices:

CC The classic Algorithm 1 of Section 3.1.

SI 1

The sparse inverse algorithm of Section 3.3 with the normal matrix in standard ordering.

SI 2 Same as SI 1 , but with the normal matrix in arrowhead ordering.

IC The inverse Cholesky Algorithm 2 of Section 3.2.

IC-S

The inverse Cholesky Algorithm 2 where the sparse-todense-conversion step 4 was omitted and all subsequent operations were performed using sparse linear algebra.

LD The legacy DBAT Algorithm 3 of Section 3.4.

Experiments

A total of four experiments were performed:

Self-calibration (SC) The unchanged normal matrices of each data set of Section 4.1 were given to each of the algorithms of Section 4.2.

No self-calibration (NSC) Same as the SC experiment, except the normal matrices were stripped of the IO blocks before processing.

Total variance only (TV) The SC and NSC experiments were repeated with the diagonal-only-variants of the CC and IC algorithm described in Section 3.5.

The effect of the number of cores The SC experiment was repeated with the CC, SI 1 , and IC algorithms where the number of available CPU cores were varied from one to six.

In all cases, the execution time for each algorithm on each data set was recorded. The normal matrix was assumed to have been computed by the last iteration of the preceding bundle adjustment. Thus, the time to compute the normal matrix was not included. Furthermore, the time to permute the normal matrix to the ordering preferred by each algorithm was excluded. Additionally, the density of the L22 and K21 blocks of the IC algorithm were recorded, as well as the density of the K22 block.

The timings were performed on an HP Z440 workstation from 2015 with a 6-core Intel Xeon E5-1650 v3 @ 3.50GHz and 64GB of RAM running Matlab 2019b Update 1 under Debian 10.

Results

The density numbers are given in Table 2. The L22 density varied between 30% and 50%. Except for the camera calibration data set, the K21 density varied between 53% and 66%. All K22 densities were very close to 50%.

The execution times for the CC, SI 1 , IC, and LD algorithms for the SC and NSC experiment are given in Table 3. The CC execution times were dominated (80-90%) by the computation of the 3-by-3 blocks in step 8. For the large problems (data sets 6-10), the IC execution times were dominated by (50-85%) by the K21 (Q) computation in step 5. The execution times for the SI 2 algorithm was within 15% of the SI 1 algorithm. The IC-S algorithm was several orders of magnitude slower than the IC algorithm. The results show that the IC algorithm is fastest on SC data sets 1-8 and on NSC data sets 1-7. The CC algorithm is fastest on the remaining data sets, with the SI algorithm somewhere in between. In both cases, the cross-over point correspond to a B density of about 2%. Except for the camera calibration data set, the IC algorithm is 20-120 times faster than the legacy LD algorithm. For the TV experiment, also presented in Table 3, the IC algorithm was fastest for SC data sets 1-7 and NSC data sets 1-6. In this case, the cross-over point corresponds to a B density of about 2.5%.

On average, the execution times for the SC data sets were 40-50% higher for the CC and SI algorithms compared to their NSC counterparts. In contrast, with the exception of the camera calibration data set, the SC and NSC times were within 1% of each other. In the diagonal-only TV experiment, the execution times decreased by about 40% for the CC algorithm and about 20% for the IC algorithm compare the full-block SC and NSC results.

The Sewu and UmU-3D-R data sets has the same number of images and object points. However, the number of image points are 38% higher in the Sewu data set. The same difference is seen in the number of rays per object point and the B-block density for the NSC data set, with a slightly lower difference for the SC data set. The increased point density translated to a 75-100% longer execution time for the CC algorithm. In contrast, the corresponding increase for the IC algorithm was about 6%.

The result of the speedup Experiment 4 are shown in Figure 5.

Only the IC algorithm benefited from an increased number of cores, but with a decreasing efficiency as the number of cores are increased. Using six cores the speedup was about 2.5. On a single core, the IC algorithm was faster the the CC algorithm for data sets 1-4, and within a factor of 2 for data sets 1-8.

DISCUSSION

In this paper, several algorithm for the computation of posterior covariance of object points were applied to self-calibration and non-self-calibration data sets of varying sizes, ranging from 21-944 images. On the smallest data sets, the fastest algorithm was the Inverse Cholesky (IC) algorithm. On the largest data sets, the classic (CC) algorithm was fastest. The performance of the Sparse Inverse (SI) algorithm was somewhere in between CC and IC. The relative results between SI and IC are consistent with the results by [START_REF] Kallin | Computation of posterior covariances of object points in bundle adjustment[END_REF], on similar problem sizes. Compared to the legacy algorithm in DBAT, the IC algorithm was 1-2 orders of magnitude faster. Several factors affect the computation time for a given algorithm, including the number of images, object points, and image points. A critical parameter appears to be the density of the mixed block of the normal matrix. As the image and point numbers increase, the density of the mixed block typically decreases. When the full 3-by-3 blocks of the posterior covariance matrix are computed, the density cross-over point for the fastest algorithm appear to be around 2%.

If the posterior correlations of the estimated object points are needed, the full 3-by-3 block must be computed. Otherwise, if only the total variance of each point is needed, only the diagonal elements of the posterior covariance have to be computed. In this case, a decrease in execution time of about 40% for the CC algorithm and 20% for the IC algorithm was observed, compared to the full block computations and the density cross-over point shifted to about 2.5%.

The IC algorithm is written completely in the Matlab language. As such, it benefits from an increased number of available cores due to the efficiency of the underlying numeric libraries used by Matlab. The performance increase was sub-linear. In contrast, the SI and CC algorithms were partially implemented in the C language and called from Matlab. The C code was not optimised for multiple cores. Further performance improvements are likely if the C code is parallellised, although the level is uncertain.

For a general toolbox as DBAT, a pure Matlab implementation is advantageous, as it removes any dependencies on, e.g., compilers. However, for a specific lab with mostly large projects, compiling and using the classic algorithm will give improved performance.

Compared to the previous version of DBAT, users can expect posterior covariance computations that are a few orders of magnitude faster. In summary, the computation time for the largest, 944-image dataset in this paper decreased by a factor of over 100. The actual computation time was around three minutes.

CONCLUSIONS

On medium-sized data sets of up to about 300 images, the Matlab-only IC algorithm is faster than the classic algorithm. On larger data sets, the classic algorithm is faster. However, on data sets up to about 900 images, and if a computation time of a few minutes at the end of a bundle adjustment processing is acceptable, the IC algorithm is still a useful high-level alternative.

We have shown throughout this paper that the novel IC algorithm implemented in DBAT managed to accelerate its computation of covariance values significantly. Experiments by comparing it to other methods, including the classical approach, showed that the results are comparable, albeit with a few caveats. For DBAT this is a significant addition as it enables a relatively rapid computation of more statistical metrics, further reinforcing its application for reprocessing bundle adjustment results of black-box solutions. Further investigations and experiments in this regard may also reveal other options for optimisation.

Figure 2 :

 2 Figure 2: If the object points are ordered first, the normal matrix forms an arrowhead matrix, where the Noo block forms the shaft of the arrow, and the remaining blocks form the tip at the lower and right borders of the matrix (left). The arrowhead shape is more clear for real projects, where the number of OP make the Noo block dominate the matrix. The right figure shows the normal matrix for the ROMA data set (see Section 4.1) with 60 images and about 25000 object points.

) 2

 2 Algorithm The inverse Cholesky algorithm for computing posterior covariance for object points. Require: N has arrowhead ordering 1: L ← chol(N ) Sparse Cholesky factorisation 2: P ← L -1 11 Sparse block inverse 3: S ← L21P Sparse multiplication 4: T ← full(L22) Convert to full 5: Q ← -T -1 S Dense triangular solve 6: d ← colsum(P P ) Sparse 7: d ← + colsum(Q Q) Dense 8: for i = 1, 2, 3 do 9: Pi ← extract(P, i) Sparse 10: Qi ← extract(Q, i) Dense 11: end for 12: for (i, j) = (2, 1), (3, 1), (3, 2) do 13: sij ← colsum(Pi Pj) Sparse 14: sij ← + colsum(Qi Qj) Dense 15: end for 16: C ← assemble(d, s21, s31, s32) Assemble output Algorithm 3 The legacy DBAT algorithm for computing posterior covariance for object points. Require: N has arrowhead ordering 1: L ← chol(N ) Sparse Cholesky factorisation 2: for all object points i in steps of 30 do 3: E ← columns from I.

Figure 5 :

 5 Figure5: Average speedup as a function of the number of cores. Only the IC algorithm benefits from an increased number of cores, but at a diminishing rate as the number of available cores increases.

Table 1 :

 1 Statistics for the data sets used in the paper.

	Ray	Point

Table 2 :

 2 Block non-zero density for the data sets. The B column correspond to the B block of the original problem in equations (22) and (5), respectively. The K21 (Q), and L22 columns are densities for corresponding blocks in the IC algorithm. All recorded K22 densities were very close to 50%.

		Self-calibration	No self-calibration
	Data Set Name	B	K21	L22	B	K21	L22
	1 Camcal	98.8 100.0 50.4 98.7 100.0 50.4
	2 Vexcel	10.9	62.8 44.0	8.0	61.6 43.5
	3 Roma	8.0	54.8 31.0	5.7	53.6 30.0
	4 StPierre	5.1	66.6 47.4	4.6	66.5 47.4
	5 UmU-2D	4.0	55.3 31.8	3.1	54.9 31.4
	6 Sewu	3.5	55.6 33.6	2.9	55.3 33.4
	7 UmU-3D-R	2.7	62.4 38.0	2.1	62.1 37.8
	8 UmU-3D	2.2	63.8 41.0	1.7	63.6 40.9
	9 UmU-2-LYR	1.0	58.5 43.4	0.8	58.4 43.3
	10 UmU-3-LYR	0.7	59.2 45.4	0.6	59.2 45.4

Table 3 :

 3 Execution times in seconds for the SC, NSC, and TV experiments. The fastest algorithm within each block is highlighted. Overall, the IC algorithm is the fastest for small data sets, whereas the CC algorithm is fastest for large data sets. The SI algorithm has intermediate performance. The cross-over point depends on the requested computation and whether the data set is a selfcalibration data set or not. The IC algorithm is faster than the legacy LD algorithm by 1-2 orders of magnitude.

	Full 3-by-3 blocks	Diagonal only

https://www.photomodeler.com/

https://www.agisoft.com/


François Schmitt

Olivier Piccin

Laurent Barbé

Bernard Bayle

Franc ¸ois Schmitt

An Origami-Inspired Flexible Pneumatic Actuator

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

I. INTRODUCTION

Soft materials offer new possibilities to design robot technologies well adapted to interactions with their environment, with new solutions to access confined spaces [START_REF] Calderón | Design, fabrication and control of a multi-material-multi-actuator soft robot inspired by burrowing worms[END_REF], grasp variable shape objects [START_REF] Correll | Soft Autonomous Materials-Using Active Elasticity and Embedded Distributed Computation[END_REF] and adapt to more situations [START_REF] Florez | Soft pneumatic actuators for legged locomotion[END_REF]. Soft robots are also known for their high resilience [START_REF] Tolley | An untethered jumping soft robot[END_REF], [START_REF] Terryn | Development of a self-healing soft pneumatic actuator: A first concept[END_REF] and good power to weight ratio [START_REF] Robertson | Soft Pneumatic Actuator Fascicles for High Force and Reliability[END_REF]. They may offer other interesting features, such as variable stiffness [START_REF] Stilli | A Novel Concept for Safe, Stiffness-Controllable Robot Links[END_REF], or zero backlash thanks to monolithic designs. However, their structural softness is not well adapted to actuate conventional rigid robotic structures.

In parallel to the development of soft robotics, new manufacturing methods have spread very fast. This is in particular the case of additive manufacturing which offers nearly unlimited possibilities to build shapes that were impossible to manufacture using conventional methods. More recently, Multi-Material Additive Manufacturing (MMAM) has enabled to combine heterogeneous materials. This is especially interesting for soft robots design, as it allows mixing soft and rigid properties in a monolithic part [START_REF] Hiller | Automatic Design and Manufacture of Soft Robots[END_REF], [START_REF] Bruyas | Toward unibody robotic structures with integrated functions using multimaterial additive manufacturing: Case study of an MRI-compatible interventional device[END_REF]. Compared to conventional plastics industry processes, it does not require complex multi-step processes such as overmolding or lamination.

Manufacturing soft robots using MMAM is not easy however. Soft robots often employ highly anisotropic materials, such as paper or fabric sheets used to reinforce the softer structure [START_REF] Sun | Characterization of silicone rubber based soft pneumatic actuators[END_REF], [START_REF] Mosadegh | Pneumatic Networks for Soft Robotics that Actuate Rapidly[END_REF]. The soft materials used are generally silicon polymers, exhibiting low stiffness and very high breaking elongation [START_REF] Ilievski | Soft Robotics for Chemists[END_REF]. When using MMAM however, the designer has a limited choice in the materials and manufacturing process. Even though commercially available MMAM materials exhibit elastic behavior comparable to other classical materials, they tend to break at much lower stress levels. This drawback becomes even more pronounced when trying to scale down the parts size, as commercially This work was supported by French state funds managed by the ANR within the Investissements d'Avenir program (Robotex ANR-10-EQPX-44, Labex CAMI -ANR-11-LABX-0004) and by the Région Grand Est. 1 F. Schmitt (f.schmitt@unistra.fr), L. Barbé In this paper, a new kind of soft fluidic actuator is proposed (see Fig. 1), inspired by recent developments in origami robots. It is based on a prismatic deformable pocket, built as a quasi-monolithic block using MMAM. This technology allows force generation for short stroke movements on a lightweight structure, with the motivation to develop building blocks for collaborative robotics where high power to weight ratio is more important than positioning accuracy. In this context, monolithic MMAM parts offer interesting opportunities to simplify the building process for such an actuator while keeping its footprint low.

In Section II, the geometry of the actuator as well as its parametrization are presented, followed in Section III by the study of the influence of several parameters on the theoretical performances of the actuator. The actuator implementation is then presented in Section IV, and the performances of the resulting prototypes are then assessed in Section V. Finally, in Section VI, the obtained results and the perspective of this work are discussed.

II. DESIGN PRINCIPLES

Origami designs have seen recent spikes in interest as they allow to build lightweight and stiff structures, while also being able to reconfigure themselves [START_REF] Miura | Synthesis of rigid-foldable cylindrical polyhedral[END_REF], [START_REF] Tachi | Rigid-foldable cylinders and cells[END_REF], [START_REF] Filipov | Origami tubes assembled into stiff, yet reconfigurable structures and metamaterials[END_REF], [START_REF] Liu | Three-Dimensional Printable Origami Twisted Tower: Design, Fabrication, and Robot Embodiment[END_REF]. Those properties are of particular interest in the design of fluidic actuators [START_REF] Li | Fluidic origami: A plant-inspired adaptive structure with shape morphing and stiffness tuning[END_REF], [START_REF] Kim | Origami Structured Compliant Actuator (OSCA)[END_REF], [START_REF] Li | Fluid-driven origamiinspired artificial muscles[END_REF] as they can be compact, lightweight, offer relatively large ranges of displacement and can also exhibit large off-axis stiffness [START_REF] Paez | Design and Analysis of a Soft Pneumatic Actuator with Origami Shell Reinforcement[END_REF]. These systems have inspired the system presented in the following.

A. Description and parametrization

The geometry of the proposed actuator can be modeled as a polyhedron and is parametrized as presented in Fig. 2. It is composed of three sub-volumes, the central one being a prismatic body obtained by the extrusion of a honeycomb pattern along the (O 0 , z) axis. The other two volumes are the end caps formed by a combination of both trapezoidal and triangular faces, in a pattern similar to the origami "waterbomb" [START_REF] Onal | An Origami-Inspired Approach to Worm Robots[END_REF], [START_REF] Qiu | Repelling-Screw Based Force Analysis of Origami Mechanisms[END_REF]. The obtained geometry exhibits three planar symmetries thus the actuator geometry can be defined by a set of five constant parameters (a, b, g, h, w) and the variable parameter e describing the actuator's extension. Moreover, two variable parameters θ and ϕ are introduced to describe the configuration of the prismatic shape and the two caps. When the angle θ > 0 (resp. < 0), the prismatic part profile is convex (resp. concave) and when ϕ > 0 (resp. < 0), the end caps are re-entrant (resp. salient) with respect to the prismatic volume. The design objective is to obtain an articulated polyhedron with rigid faces in which the edges act as hinges, allowing the faces to rotate with respect to their neighbors. Moreover, the polyhedron should provide a variable internal volume during its operation. However, following the bellows conjecture [START_REF] Connelly | The Bellows conjecture[END_REF], the internal volume of an articulated polyhedron remains constant under rigid deformations and thus it is mandatory that some faces or edges of the polyhedron change their geometry during operation. In order to better guide the motion, deformations should be localized in the end caps. Thus, the length f of the edge O 2 S 1 was chosen as a variable dimension to accommodate the volume variation.

a = O1O2 θ = x, # » O1O2 b = O1O6 ϕ = x, # » H1S1 g = S1S2 e = O1O3 h = H1S1 f = O2S1 w = O1O7

B. Geometric model and deployment limits

The deployment of the proposed polyhedron can be described using one parameter chosen among e, θ, ϕ and these variable parameters are related by the two equations e = 2a cos θ e = 2h cos ϕ

Considering that the faces have no thickness, the polyhedron deployment is limited by the following conditions:

• flat folded: when e = 0 (or θ = ϕ = ± π 2 ). • fully unfolded: either when h > a, the point O 2 is on

[O 1 O 3 ]
leading to e = 2a and θ = 0 (fully unfolded prismatic shape condition), or when h < a, the point

S 1 is on [H 1 H 2 ]
leading to e = 2h and ϕ = 0 (fully unfolded caps condition). When h = a, both conditions are met simultaneously.

Self collisions may only occur for concave configurations of the prismatic shape or with re-entrant caps, namely when b < 2a or when w < 2h. In this case, self collisions happen when θ = -arcsin b a or ϕ = arcsin w h . To summarize, the output limits [e min , e max ] of the deployable polyhedron can be calculated as:

e min = max      4a 2 -b 2 (if b < 2a) 4h 2 -w 2 (if w < 2h) 0 e max = min(2a, 2h) (2)

III. PARAMETERS SELECTION

In order to study and improve the performances of the actuator, the selection of the dimensional parameters is not limited to its mechanical limits. Indeed, several factor have to be accounted for, ranging from the base geometry to considerations related to the implementation of the system.

A. Internal volume

The analysis of the chamber internal volume gives two characteristics of particular interest in order to predict the behavior of the actuator: i) the fluid flow required in order to actuate the structure, and ii) the inflation (resp. deflation) of the structure reaches a stable configuration when the internal volume reaches a local maximum (resp. minimum). Using both polyhedron parametrization and equations ( 1), the internal volume can be computed as

Vact = V prism ew b + sign (θ) 4a 2 -e 2 -sign (ϕ) e 6 4h 2 -e 2 2b + g + sign (θ) 4a 2 -e 2 Vcaps ( 3
)
where V prism and V caps represent the volume of the prismatic central part of the actuator and the combined volume of both end caps (counted positively). Figure 3 plots the evolution of these volumes in function of the extension e, for the selected set of dimensions given in Table I. As expected, V prism has the most significant influence on the overall volume of the polyhedron. The figure also shows that in the case of a concave prismatic section (θ < 0, dashed lines), the volume increases monotonously with respect to the deployment whereas in the case of the convex section (θ > 0, solid lines) it reaches a maximum before the full extension (e = 2a), represented by the vertical dashed black line. One should note however that if θ is allowed to change sign, the actuator should naturally switch configuration from θ < 0 and reach the stable state corresponding to the maximum volume when θ > 0. Finally, concave configurations for either the prismatic section or the caps also reduce the footprint of the actuator, all other dimensions being equal.

B. Study of the edge of variable length f

As presented in section II-A, in order to allow the actuator's extension, the length f of the four edges O 2 S 1 is I.

assumed to be variable during operation and its expression can be written as

f = 1 2 4h 2 -e 2 + b -g + sign (θ) 4a 2 -e 2 2 . (4)
The figure 4a) shows the evolution of f in function of e for a convex (orange color) and a concave (blue color) shape of the actuator with the parameters of table I. One can see that the variation of f over the range of e is smaller for the convex (θ < 0) configuration. Furthermore, in the practical implementation of these extensible edges, this variation in length will generate traction/compression stress on the material, combined with a flexion due to the relative rotation of the adjacent faces can lead to the rupture of the internal pocket. Limiting this variation could improve the robustness of the design and therefore the concave (θ < 0) configuration is more desirable. The figure 4b) shows the evolution of f in the concave configuration for different values of h which drive the size of the caps. In order to define a design criterion, we form the quantity

f var = emax 0 (f -f min )de (5)
where f min is the minimum of f in the [0; e max ] range. Minimizing f var leads to a smaller amplitude of variation for f when e varies. The criterion f var depends on the construction parameters a, b, g and h. However, for a given bulkiness of the actuator which is mostly governed by the values of a and b, the remaining parameters are g and h and their influence can be explored to minimize f var . Figures 5a) and 5b) give the plots of f var as a function of h and g and provide the following design insights:

• higher values for h should be preferred as it also leads to smaller variations of f but this comes at the cost of increasing the footprint of the actuator as it increases the size of the end caps. • lower values of g leads also to smaller variations of f . However, when g = 0, S 1 and S 2 are coincident, this leads to the usual waterbomb pattern: a 6-th order vertex. Higher order vertices may have a detrimental effect on the resistance of the actuator as they will combine flexions coming from a higher number of edges and thus become a critical point in the structure. This problem is generally tackled by removing the vertex, but as the closure of the volume is mandatory for our application, we opted for low positive values of g, creating two 4-th order vertices.

C. Chosen parameters

Considering the elements presented in the previous sections, but also other practical limitations imposed by the manufacturing process, including but not limited to precision and material resistance, we have chosen to build the proposed prototype using the parameters indicated in Although not optimal for limiting the variations of f , h has been chosen small in order to favor compactness. In addition, a parameter named e con has been added, defining the value of e at construction time. This parameter influences the displacement range and also the strain on the hinges of the actuators. Two version are proposed with two different values for e con : i) Version A built with e con = 22 mm, and ii) Version B built with e con = 25 mm. Both versions have been prototyped and evaluated.

IV. PROTOTYPE IMPLEMENTATION

A. Manufacturing method

In the last section, the proposed actuator was presented as an ideal model in which most motions are rigid with the only exception of the extensible edges. In this model, the faces were also supposed infinitely stiff and thin. Obviously all of these assumptions cannot be reproduced on an actual prototype where material and manufacturing constraints have to be taken in account. Although several manufacturing methods exist to produce origami-like geometries [START_REF] Martinez | Elastomeric Origami: Programmable Paper-Elastomer Composites as Pneumatic Actuators[END_REF], [START_REF] Zhakypov | Design Methodology for Constructing Multimaterial Origami Robots and Machines[END_REF], we chose both for convenience and availability to use a multimaterial additive manufacturing machine in order to produce our prototype.

The machine used, a Connex 350 (Stratasys Ltd. [START_REF]Stratasys, Polyjet Process[END_REF]), allows to use a rigid polymer (VeroWhiteplus RGD 835 (VW)) and a "rubber-like" elastomer (TangoBlackplus FLX 980 (TB)), or combinations of both in order to build geometries with varying mechanical properties. This allows the design of prototypes with complex geometries that can be produced with a minimum amount of machining and assembly operations, but at the cost of material performances that are generally inferior to material used in typical molding processes.

B. Flexible origami-like compliant mechanism

Although the proposed actuator design may work with any kind of fluid, we chose to use compressed air as actuation fluid. Thus, when considering the system, several functions have to be included on the prototype in order to work as initially intended: i) "undeformable" faces should be as rigid as possible, ii) hinges located at the edges and deformable faces have to be as compliant as possible, iii) the internal pocket needs to be airtight and iv) the overall structure should be able to withstand both internal pressure and material deformation. To this end, the CAD model has been created from the geometrical concept by building first a rigid skeleton using 4 mm-thick VW plates and then adding flexures between the faces using 1.5 mm thick living hinges in TB to link all rigid faces with their relative neighbors. Deformable faces have also been built in TB, with a thickness of 2 mm as a trade off between compliance and resistance. The obtained CAD model is presented in Fig. 6.

From this base concept, both the base and end plates have then been cut open to provide an access to the internal room for manufacturing purpose. During assembly phase, two rigid closure caps and a polyurethane inlet tube are glued on the semi-flexible structure in order ensure the internal room tightness and allow exchanges of pressurized air with an external control system. In red are indicated the flexible regions in the section planes.

C. Deployment limit

The inflation of pressurized air in the internal room tends to push the faces of the prismatic section toward the exterior of the actuator, allowing them to fold from a concave to a convex configuration, as it would allow the structure to reach its maximum volume (as discussed in Section III-A). However, this behavior is not desired as: i) it reduces the extension of the actuator and ii) it may increase the stress on the living hinges of the mechanism. Therefore we added on the design stopping geometries (see Fig. 6) in order to limit the rotation of the legs without affecting the deployment of the structure along the (O 0 , x) axis.

D. Experimental testbed

During characterization experiments, the prototypes have been controlled using a computer controlled pneumatic testbed as shown in Fig. 7. The input pressure is regulated by a proportional valve (VPPM series, Festo) and a fast 3/2-way solenoid valve (MHE2 series, Festo) triggers the inflation or deflation of the actuator. The pressure of each line is also measured with a pressure sensor (SPTW series, Festo). Displacements of the actuator are then measured by a laser telemeter (optoNCDT 2300 series, Micro-Epsilon) and interaction forces are measured using a 100N force sensor (K1563, Scaime) as shown in Fig. 7. Custom software running under Linux Xenomai real-time operating system is used to synchronize all sensor acquisitions with the control of servo-valve and solenoid valves. V. EXPERIMENTAL CHARACTERIZATION Before starting the experimental characterization, a preliminary study not presented here has been carried out to assess the maximum pressure before failure of both versions of the proposed actuator. Both versions broke at similar pressure levels, and as such this gave us an estimated safe working pressure of 24.9 kPa that we used as the maximum pressure for all following tests. The minimum pressure level used was determined by the hardware as we could not generate pressure lower than 7.44 kPa. The first experiment carried out on the actuator was a displacement characterization. In this test, the actuator was allowed to move freely without load, and several steps of increasing pressures were applied. The results displayed in Fig. 8 show on the left the system transient response for the minimum and maximum pressures, where ε = e-econ econ represents the deformation relative to the initial size of the actuator (expressed in percentage). On the right is plotted the maximum relative deformation of the actuator at several pressure levels.

A. Displacement characterization

The observed deformation show a maximum of 51.3% (resp. 32.7%) for the version A (resp. B) at a pressure of 24.9 kPa, corresponding to a maximum displacement of about 11.3 mm (resp. 8.2 mm). In both cases, the maximum value for e is superior to its theoretical maximum of 32 mm (see equation ( 2)). This phenomenon is due to a deformation along the x axis of the living hinges that had not be taken in account in the geometric model.

B. Force characterization

The second experiment carried out on the actuator was a static force characterization. In this test, the extension of the actuator was blocked by a force sensor in order to measure the force generated by the system in its constructed position under several steps of pressure. Indeed as only a limited amount of energy is stored in material elasticity, the force generated at the constructed configuration should be maximal. The results displayed in Fig. 9 show on the left the transient response of the system for the minimum and maximum pressure. On the right is plotted the maximum force observed at each pressure levels. At low pressure, both systems generate similar levels of force, but as the pressure raises, the difference between both versions of the system increases. At the maximum pressure, version A is able to generate 33.2 N (resp. 40.0 N for version B). By comparison, a frictionless cylinder with a similar footprint (14 cm 2 for the base of the actuator) would generate a force of 34.9 N at 24.9 kPa.

C. Cyclic loading

Both of the previous experiments have been carried out by using only one pressure step for each pressure level. In order to characterize the repeatability and the resistance of the system, a cyclic test has been performed. Both displacement and force generation have been tested using a periodic step pressure load with a period of approximately 8s (due to software limitation, the triggering was performed manually). The selected pressure for this test was 17.4 kPa in order to limit the stress on the structure while generating significant displacements and forces. Both versions of the prototype have been tested and the displacement (or forces) were measured at each cycle. Table II shows the mean and standard deviation of the maximum forces and displacements measured during the tests. Although version A (resp. B) broke after 26 (resp. 30) cycles during the cyclic displacement test, both version were still in working condition after more than 60 cycles in the force tests. This difference is probably due to the lower level of material deformation during force tests with respect to displacements tests.

D. Experimental results discussion

In parallel to the quantitative results shown in this section, some observations can be made on the overall experimental assessment. No dynamic experiment have been carried out as they are influenced by both the actuator and the air supply performance. It has been observed however that at lower working pressures, the air pressure supplied by the testbed rises with an unrepeatable profile, leading some experimental results that were discarded. Changing the servo-valve for a better suited model should improve further experiments. Testing the off-axis stiffness and on-axis force displacement relationship would also be of interest. It would however require a more complex setup than currently available as forces need to be measured while the displacement are constrained (or reciprocally) along several directions. Finally, a limitation of the manufacturing process has also been observed during the tests. The MMAM process used creates ridges aligned with the horizontal plane of printing. These fracture initiations, combined with mixed tensile and flexural deformations along the living hinges, drastically limit the lifetime of the actuator.

VI. CONCLUSIONS

In this paper, a new robotic actuator was proposed, based on a monolithic design combining both stiff and soft elements. Two variants have been built and their performances have been assessed for both displacements and force generation. Forces in the ranges of 40N at a pressure of 24.9 kPa were measured as well as relative deformations without load in the order of 51.3% starting from the initial actuator size. Although both versions of the prototype exhibit a limited number of work cycles, especially at higher deformation ranges, they perform with a decent robustness and repeatability considering the limitations of the manufacturing process.

Further work will focus on improving the current design by trying to reduce extension at the living hinges in order 1) to increase their resistance; 2) to be closer to the theoretical model of the system, improving the system mechanical accuracy. In order to solve this issue we will adapt the present concept to other manufacturing methods such as overmolding on composite structure, in order to combine fully rigid parts, flexible fibrous materials and elastomers. Improvements of the testbed are also required in order to carry out more indepth assessments of the system.

The proposed prototype is also only capable of producing forces and displacements in one direction. We will work on implementing a bidirectional system by either using two actuators in antagonistic configurations, or by implementing a positive/negative pressure control with a single actuator in order to favor a more compact actuation solution.

Fig. 1 .

 1 Fig. 1. General view of the actuator prototype (Version B). The ruler is graduated in centimeters.

Fig. 2 .

 2 Fig. 2. Geometry and parametrization of the proposed actuator.

3 ]Fig. 3 .

 33 Fig. 3. Evolution of Vcaps and Vact with respect to e. The evaluation parameters chosen are those of Table I, with negative and positive values for θ and ϕ.

Fig. 4 .

 4 Fig. 4. Evolution of f with respect to e with a) different convexity for the prismatic part, and b) various values for the h parameter considering θ < 0. The other evaluation parameters chosen are those of TableI.

Fig. 5 .

 5 Fig. 5. Evolution of fvar with respect to e and with various values of a) h, and b) g. The other evaluation parameters used for each plot are those of Table I and θ < 0.

Fig. 6 .

 6 Fig. 6. CAD model of the proposed design (Version B) in full view and quarter section view with detailed views of the hinges implementation. Flexible regions (TB) are represented in black and rigid ones (WV) in white. In red are indicated the flexible regions in the section planes.

Fig. 7

 7 Fig. 7. a) Schematic layout of the pneumatic control testbed. Experimental setup for: b) Displacements measurements, and c) Forces measurements.

Fig. 8 .

 8 Fig. 8. a) Displacement of both versions of the unloaded actuator at minimum and maximum pressure with respect to time. b) Maximum displacement of the unloaded actuator with respect to different pressures.

FittingFig. 9 .

 9 Fig. 9. a) Static force response generated by both versions of the actuator at minimum and maximum pressure with respect to time. b) Maximum force generated by the actuator the actuator with respect to different pressures.

  and B. Bayle are with ICube Laboratory, University of Strasbourg-CNRS, France.

2

O. Piccin is with ICube Laboratory, and INSA Strasbourg, France available MMAM machines generally have a manufacturing accuracy of about a tenth of millimeter.

  Table. I.

	a [mm]	b [mm]	g [mm]	h [mm]	w [mm]	θ	ϕ
	16	35	5	16	40	< 0	< 0
			TABLE I
		SELECTED PROTOTYPE PARAMETERS.