Kaustuv Chaudhuri
email: kaustuv.chaudhuri@inria.fr

Undecidability of Multiplicative Subexponential Logic

Subexponential logic is a variant of linear logic with a family of exponential connectives-called subexponentialsthat are indexed and arranged in a pre-order. Each subexponential has or lacks associated structural properties of weakening and contraction. We show that classical propositional multiplicative linear logic extended with one unrestricted and two incomparable linear subexponentials can encode the halting problem for two register Minsky machines, and is hence undecidable.

⊢ a, ¬ a init ⊢ Γ, A ⊢ ∆, B ⊢ Γ, ∆, A ⊗

Introduction

The decision problem for classical propositional multiplicative exponential linear logic (MELL), consisting of formulas constructed from propositional atoms using the connectives {⊗ ⊗ ⊗, 1, , ⊥ ⊥ ⊥, !, ?}, is perhaps the longest standing open problem in linear logic. MELL is bounded below by the purely multiplicative fragment (MLL), which is decidable even in the presence of first-order quantification, and above by MELL with additive connectives (MAELL), which is undecidable even for the propositional fragment [START_REF] Lincoln | Decision problems for propositional linear logic[END_REF]. This paper tries to make the undecidable upper bound a bit tighter by considering the question of the decision problem for a family of propositional multiplicative subexponential logics (MSEL) [START_REF] Nigam | Exploiting non-canonicity in the sequent calculus[END_REF][START_REF] Nigam | Algorithmic specifications in linear logic with subexponentials[END_REF], each of which consists of formulas constructed from propositional atoms using the (potentially infinite) set of connectives {⊗ ⊗ ⊗, 1, , ⊥ ⊥ ⊥} ∪ u∈Σ ! u , ? u , where Σ is a pre-ordered set of subexponential labels, called a subexponential signature, that is a parameter of the family of logics. In particular, we show that a particular MSEL with a subexponential signature consisting of exactly three labels can encode a two register Minsky machine (2RM), which is Turing-equivalent. This is the same strategy used in [START_REF] Lincoln | Decision problems for propositional linear logic[END_REF] to show the undecidability of MAELL, but the encoding in MSEL is different-simpler-for the branching instructions, and shows that additive behaviour is not essential to implement branching. We use the classical dialect of linear logic to show these results. The intuitionistic dialect has the same decision problem because it is possible to faithfully encode (i.e., linearly simulate the sequent proofs of) the classical dialect in the intuitionistic dialect without changing the signature [START_REF] Chaudhuri | Classical and intuitionistic subexponential logics are equally expressive[END_REF].

This short note is organized as follows: in section 2 we sketch the one-sided sequent formulation of MSEL and recall the definition of a 2RM. In section 3 we encode the transition system of a 2RM in a MSEL with a particular signature. In section 4 we argue that the encoding is adequate, i.e., that the halting problem for a 2RM is reduced to the proof search problem for this MSEL-encoding, by appealing to a focused sequent calculus for MSEL. The final section 5 discusses some of the ramifications of this result.

Background

Propositional Subexponential Logic

Let us quickly recall propositional subexponential logic (SEL) and its associated sequent calculus proof system. This logic is sometimes called subexponential linear logic (SELL), but since it is possible for the subexponentials to have linear semantics it is redundant to include both adjectives. Formulas of SEL (A, B, . . . ) are built from atomic formulas (a, b, . . . ) according to the following grammar:

atomic multiplicative additive subexponential A, B, . . . ::= a | A ⊗ ⊗ ⊗ B | 1 | A ⊕ ⊕ ⊕ B | 0 | ! u A | ¬ a | A B | ⊥ ⊥ ⊥ | A & B | ⊤ ⊤ ⊤ | ? u A
Each column in the grammar above is a De Morgan dual pair. A positive formula (depicted with P or Q when relevant) is a formula belonging to the first line of the grammar, and a negative formula (depicted with N or M) is a formula belonging to the second line. The labels (u, v, . . . ) on the subexponential connectives ! u and ? u belong to a subexponential signature defined below. The additive fragment of this syntax is just used in this section for illustration; we will not be using the additives in our encodings. The fragment without the additives will be called multiplicative subexponential logic (MSEL).

(u ∈ U Σ ) ⊢ Σ Γ ⊢ Σ Γ, ? u A weak (u ∈ U Σ ) ⊢ Σ Γ, ? u A, ? u A ⊢ Σ Γ, ? u A contr
Definition 1. A subexponential signature Σ is a structure Λ, U, ≤ where:

• Λ is a countable set of labels;

• U ⊆ Λ, called the unbounded labels; and • ≤ ⊆ Λ × Λ is a pre-order on Λi.e., it is reflexive and transitive-and ≤-upwardly closed with respect to U, i.e., for any u, v ∈ Λ, if u ∈ U and u ≤ v, then v ∈ U.

We will assume an ambient signature Σ unless we need to disambiguate particular instances of MSEL, in which case we will use Σ in subscripts. For instance, MSEL Σ is a particular instance of MSEL for Σ.

The true formulas of MSEL are derived from a sequent calculus proof system consisting of sequents of the form ⊢ A 1 , . . . , A n (with n > 0) and abbreviated as ⊢ Γ. The contexts (Γ, ∆, . . .) are multi-sets of formulas of SEL, and Γ, ∆ and Γ, A stand as usual for the multi-set union of Γ with ∆ and {A}, respectively. The inference rules for SEL sequents are displayed in figure 1. Most of the rules are shared between SEL and linear logic and will not be elaborated upon here. The differences are with the subexponentials, for which we use the following definition. Definition 2. For any n ∈ N and lists u = [u 1 , . . . , u n ] and A = [A 1 , . . . , A n ], we write ? u A to stand for the context

? u 1 A 1 , . . . , ? u n A n . For v = [v 1 , . . . , v n ], we write u ≤ v to mean that u ≤ v 1 , . . . , and u ≤ v n .
The rule for !, sometimes called promotion, has a side condition that checks that the label of the principal formula is less than the labels of all the other formulas in the context. This rule cannot be used if there are non-?-formulas in the context, nor if the labels of some of the ?-formulas are strictly smaller or incomparable with that of the principal !-formula. Both these properties will be used in the encoding in the next section. The structural rules of weakening and contraction apply to those principal ?-formulas with unbounded labels.

Two Register Minsky Machines

Like Turing machines, Minsky register machines have a finite state diagram and transitions that can perform I/O on some unbounded storage device, in this case a bank of registers that can store arbitrary natural numbers. We shall limit ourselves to machines with two registers (2RM) a and b, which are sufficient to encode Turing machines.

Definition 3. A 2RM is a structure Q, * , C, -→ where: • Q is a non-empty finite set of states; • * ∈ Q is a distinguished halting state;
• C is a set of configurations, each of which is a structure of the form q, v , with q ∈ Q and v : {a, b} → N, that assigns values (natural numbers) to the registers a and b in state q; • -→ ⊆ C × I × C is a deterministic labelled transition relation between configurations where the label set I = {halt, incra, incrb, decra, decrb, isza, iszb} (called the instructions). By usual convention, we write -→ infix with the instruction atop the arrow. We require that every element of -→ fits one of the following schemas, where in each case q, r ∈ Q and q r:

q, v halt --------→ * , {a : 0, b : 0} (with q * ) q, {a : m, b : n} incra --------→ r, {a : m + 1, b : n} q, {a : m, b : n} incrb --------→ r, {a : m, b : n + 1} q, {a : m + 1, b : n} decra --------→ r, {a : m, b : n} q, {a : m, b : n + 1} decrb --------→ r, {a : m, b : n} q, {a : 0, b : n} isza --------→ r, {a : 0, b : n} q, {a : m, b : 0} iszb --------→ r, {a : m, b : 0} (1) 2
For a trace i = [i 1 , . . . , i n ], we write q 0 , v 0

i --→ q n , v n if q 0 , v 0 i 1 --→ • • • i n --→ q n , v n .
The 2RM halts from an initial configuration q 0 , v 0 if there is a trace i such that q 0 , v 0 i --→ * , {a : 0, b : 0} . (The configuration * , {a : 0, b : 0} will be called the halting configuration.) The halting problem for a 2RM is the decision problem of whether the machine halts from an initial configuration.

The requirement that -→ be deterministic amounts to: q, v i --→ q 1 , v 1 and q, v j --→ q 2 , v 2 imply that i = j, q 1 = q 2 , and v 1 = v 2 . Note that a trace that does not end with a halting configuration will not be considered to be halting, even if there is no possible successor configuration. It is an easy exercise to transform a given 2RM into one where every configuration has a successor except for the halting configuration.

Theorem 4 ( [START_REF] Minsky | Recursive unsolvability of Post's problem of 'tag' and other topics in the theory of Turing machines[END_REF]). The halting problem for 2RMs is recursively unsolvable.

The Encoding

For a given 2RM, which we fix in this section, we will encode its halting problem as the derivability of a particular MSEL sequent that encodes its labelled transition system and the initial configuration. We will use the following subexponential signature in the rest of this section. Definition 5. Let Ξ stand for the signature {∞, a, b} , {∞} , ≤ where ≤ is the reflexive-transitive closure of ≤ 0 defined by a ≤ 0 ∞ and b ≤ 0 ∞. Definition 6 (encoding configurations). For c = q, v , we write E(c) for the following MSEL Ξ context:

? a ¬ ra, ? a ¬ ra, . . . , ? a ¬ ra length = v(a) , ? b ¬ rb, ? b ¬ rb, . . . , ? b ¬ rb length = v(b)
, ¬ q Definition 7 (encoding transitions). The transitions [START_REF] Andreoli | Logic programming with focusing proofs in linear logic[END_REF] of the 2RM are encoded as a context Π with:

• to represent q, v halt -----→ * , {a : 0, b : 0} , the elements: q ⊗ ⊗ ⊗ ¬ h, h ⊗ ⊗ ⊗ ! a ra ⊗ ⊗ ⊗ ¬ h, h ⊗ ⊗ ⊗ ! b rb ⊗ ⊗ ⊗ ¬ h, h ⊗ ⊗ ⊗ ! ∞ 1 (for some h Q):
• to represent q, {a : m, b : n} incra -----→ r, {a : m + 1, b : n} , the element q ⊗ ⊗ ⊗ (¬ r ? a ¬ ra);

• to represent q, {a : m, b : n} incrb -----→ r, {a : m, b : n + 1} , the element: q ⊗ ⊗ ⊗ (¬ r ? b ¬ rb);

• to represent q, {a : m + 1, b : n} decra -----→ r, {a : m, b : n} , the element: q ⊗ ⊗ ⊗ ! a ra ⊗ ⊗ ⊗ ¬ r;

• to represent q, {a : m, b : n + 1} decrb -----→ r, {a : m, b : n} , the element: q ⊗ ⊗ ⊗ ! b rb ⊗ ⊗ ⊗ ¬ r;

• to represent q, {a : 0, b : n} isza -----→ r, {a : 0, b : n} , the element: q ⊗ ⊗ ⊗ ! b ¬ r; and

• to represent q, {a : m, b : 0} iszb -----→ r, {a : m, b : 0} , the element: q ⊗ ⊗ ⊗ ! a ¬ r.
Note that Π contains a finite number of elements. Definition 8 (encoding the halting problem). If Γ is A 1 , . . . , A n , then let ? u Γ stand for ? u A 1 , . . . , ? u A n . The encoding of the halting problem for the 2RM from the initial configuration c 0 = q 0 , v 0 is the MSEL Ξ sequent ⊢ ? ∞ Π, E(c 0 ).

Theorem 9. If the 2RM halts from c 0 , then ⊢ Ξ ? ∞ Π, E(c 0 ) is derivable. Proof. We will show that if c = q 1 , v 1 i --→ q 2 , v 2 = d (for some i), then the following MSEL Ξ rule is derivable: ⊢ ? ∞ Π, E(d) ⊢ ? ∞ Π, E(c)
This is largely immediate by inspection. Here are three representative cases.

• The case of i = incra: it must be that v 2 (a) = v 1 (a) + 1 and v 2 (b) = v 1 (b), so E(d) = E(c) \ {¬ q 1 } , ¬ q 2 , ? a ra. Moreover, q 1 ⊗ ⊗ ⊗ (¬ q 2 ? a ¬ ra) ∈ Π. So:

⊢ ¬ q 1 , q 1 init ⊢ ? ∞ Π, E(c) \ {¬ q 1 } , ¬ q 2 , ? a ¬ ra ⊢ ? ∞ Π, E(c) \ {¬ q 1 } , ¬ q 2 ? a ¬ ra ⊢ ? ∞ Π, E(c), q 1 ⊗ ⊗ ⊗ (¬ q 2 ? a ¬ ra) ⊗ ⊗ ⊗ ⊢ ? ∞ Π, E(c) contr, ?
The cases for incrb, decra, and decrb are similar. • The case of i = isza: it must be that v 2 (a

( u ∈ U Σ ) ⊢ Σ ? u A, ¬ a, [a] [init] ( u ∈ U Σ ) ⊢ Σ ? u A, Ω 1 , [B] ⊢ Σ ? u A, Ω 2 , [C] ⊢ Σ ? u A, Ω 1 , Ω 2 , [B ⊗ ⊗ ⊗ C] [⊗ ⊗ ⊗] ( u ∈ U Σ ) ⊢ Σ ? u A, [1] [1] ⊢ Ω, [A] ⊢ Ω, [A ⊕ ⊕ ⊕ B] ⊕ ⊕ ⊕ 1 ⊢ Ω, [B] ⊢ Ω, [A ⊕ ⊕ ⊕ B] ⊕ ⊕ ⊕ 2 no rule for 0 (u ≤ Σ v) ⊢ Σ ? v A, C ⊢ Σ ? v A, [! u C] [!] ⊢ Ω, N ⊢ Ω,
(u U Σ ) ⊢ Σ Ω, [A] ⊢ Σ Ω, ? u A ldecide (u ∈ U Σ ) ⊢ Σ Ω, ? u A, [A] ⊢ Σ Ω, ? u A udecide
) = v 1 (a) = 0 and v 2 (b) = v 1 (b). Hence, E(d) = E(c) \ {¬ q 1 } , ¬ q 2 and ? a ra E(c) ∪ E(d). Moreover, q 1 ⊗ ⊗ ⊗ ! b ¬ q 2 ∈ Π. So: ⊢ ¬ q 1 , q 1 init ⊢ ? ∞ Π, E(c) \ {¬ q 1 } , ¬ q 2 ⊢ ? ∞ Π, E(c) \ {¬ q 1 } , ! b ¬ q 2 ! ⊢ ? ∞ Π, E(c), q 1 ⊗ ⊗ ⊗ ! b ¬ q 2 ⊗ ⊗ ⊗ ⊢ ? ∞ Π, E(c) contr, ?
The instance of ! is justified because b ≤ ∞ and b ≤ b, and there are no ?-formulas labelled a or non-? formulas in the sequent. The case of iszb is similar.

• The case of i = halt. Here, we know that q 1 ⊗ ⊗ ⊗ ¬ h ∈ Π, so:

⊢ ¬ q 1 , q 1 init ⊢ ? ∞ Π, E(c) \ {¬ q 1 } , ¬ h ⊢ ? ∞ Π, E(c), q 1 ⊗ ⊗ ⊗ ¬ h ⊗ ⊗ ⊗ ⊢ ? ∞ Π, E(c) contr, ?
Now, as long as there are any occurrences of ? a ra or ? a rb in E(c), we can apply one of the decrementing rules

h ⊗ ⊗ ⊗ ! a ra ⊗ ⊗ ⊗ ¬ h or h ⊗ ⊗ ⊗ ! b rb ⊗ ⊗ ⊗ ¬ h ∈ Π.
The general case looks something like this, where ∆ ra = {¬ ra, . . . , ¬ ra} and ∆ rb = {¬ rb, . . . , ¬ rb}.

⊢ h, ¬ h init ⊢ ¬ ra, ra init ⊢ ? a ¬ ra, ra ? ⊢ ? a ¬ ra, ! a ra ! ⊢ ? ∞ Π, E(c) \ {¬ q 1 , ? a ∆ ra , ? b ∆ rb , ? a ¬ ra} , ¬ h ⊢ ? ∞ Π, E(c) \ {¬ q 1 , ? a ∆ ra , ? b ∆ rb , ? a ¬ ra} , ? a ¬ ra, ¬ h, h ⊗ ⊗ ⊗ ! a ra ⊗ ⊗ ⊗ ¬ h ⊗ ⊗ ⊗, ⊗ ⊗ ⊗ ⊢ ? ∞ Π, E(c) \ {¬ q 1 , ? a ∆ ra , ? b ∆ rb , ? a ¬ ra} , ? a ¬ ra, ¬ h contr, ?
There is a symmetric case for contracting the h ⊗ ⊗ ⊗ ! b rb ⊗ ⊗ ⊗ ¬ h. Eventually, the right branch just becomes ⊢ ? ∞ Π, ¬ h, at which point we have:

⊢ h, ¬ h init ⊢ 1 1 . . . . weak ⊢ ? ∞ Π, 1 ⊢ ? ∞ Π, ! ∞ 1 ! ⊢ ? ∞ Π, ¬ h, h ⊗ ⊗ ⊗ ! ∞ 1 ⊗ ⊗ ⊗ ⊢ ? ∞ Π, ¬ h contr, ?

Adequacy of the Encoding via Focusing

By the contrapositive of theorem 9, if the sequent ⊢ Ξ ? ∞ Π, E(c 0 ) is not derivable, then the 2RM does not halt from c 0 . This gives half of the reduction. For the converse of theorem 9, we need to show how to recover a halting trace by searching for proofs of a MSEL Ξ encoding of a halting problem. The best way to do this is to build a focused proof which will have the derived inference rules in the above proof as the only possible synthetic rules, in a sense made precise below. We will begin by sketching the focused proof system for SEL that is sound and complete for the unfocused system of figure 1, and then show how the synthetic rules for the encoding are in bijection for all instructions (with a small correction needed for halt).

Focusing is a general technique to restrict the non-determinism in a cut-free sequent proof system. Though originally defined for classical linear logic in [START_REF] Andreoli | Logic programming with focusing proofs in linear logic[END_REF], it is readily extended to many other logics [START_REF] Chaudhuri | A logical characterization of forward and backward chaining in the inverse method[END_REF][START_REF] Liang | Focusing and polarization in linear, intuitionistic, and classical logics[END_REF][START_REF] Nigam | Exploiting non-canonicity in the sequent calculus[END_REF]. This section sketches the basic design of a focused version of the rules of figure 1, and omits most of the meta-theoretic proofs of soundness and completeness, for which the general proof techniques are by now well known [START_REF] Chaudhuri | A logical characterization of forward and backward chaining in the inverse method[END_REF][START_REF] Miller | From proofs to focused proofs: a modular proof of focalization in linear logic[END_REF][START_REF] Simmons | Structural focalization[END_REF]. To keep things simple, we will define a focused calculus by adding to the unfocused system a new kind of focused sequent, ⊢ Ω, [A], where the formula A is under focus. Contexts written with Ω, which we call neutral contexts, can contain only positive formulas, atoms, negated atoms, and ?-formulas. The rules of the focused proof system for SEL are depicted in figure 2.

Focused sequents are created-reading from conclusion upwards to premises-from unfocused sequents with neutral contexts by means of the rules decide, ldecide, or udecide. In a focused sequent, only the formula under focus can be principal, and the focus persists on the immediate subformulas of this formula in the premises, with the exception of the rule [!]. In the base case, for [init], the focused atom must find its negation in the context, while all formulas in the context must be ?-formulas with unbounded labels. When the focused formula is negative, the focus is released with the [blur] rule, at which point any of the unfocused rules { , ⊥ ⊥ ⊥, &, ⊤ ⊤ ⊤} of figure 1 can be used to decompose the formula and its descendants further. Eventually, when there are no more negative descendants-i.e., the whole context has the form Ω-a new focused phase is launched again and the cycle repeats. Note that the structural rules contr and weak of the unfocused calculus are removed in the focused system. Instead, weakening is folded into [init] and [START_REF] Andreoli | Logic programming with focusing proofs in linear logic[END_REF], and contraction is folded into [⊗ ⊗ ⊗] and udecide. The rules contr and weak remain admissible for either sequent form in the focused calculus.

Theorem 10. The SEL sequent ⊢ Γ is provable in the unfocused system of figure 1 iff it is provable in the focused system of figure 2.

Sketch. Straightforward adaptation of existing proofs of the soundness and completeness of focusing, such as [START_REF] Chaudhuri | A logical characterization of forward and backward chaining in the inverse method[END_REF][START_REF] Miller | From proofs to focused proofs: a modular proof of focalization in linear logic[END_REF][START_REF] Simmons | Structural focalization[END_REF]]. An instance for SEL can be found in [8, chapter 5].

Theorem 11. The 2RM halts from c 0 if ⊢ Ξ ? ∞ Π, E(c 0 ) is derivable.
Proof. We will show instead that the 2RM halts from c 0 if the sequent ⊢ Ξ ? ∞ Π, E(c 0 ) is derivable in the focused calculus, and we will moreover extract the halting trace from such a focused proof. The required result will then follow immediately from theorem 10, since any provable SEL sequent has a focused proof.

Let a focused proof of ⊢ Ξ ? ∞ Π, E(c) (for c = q, v ) be given. We proceed by induction on the lowermost instance of udecide in this proof. Note that the MSEL Ξ context ? ∞ Π, E(c) is neutral; moreover, all the elements of E(c) are either negated atoms or ?-prefixed negated atoms with bounded labels. So, the only rules of the focusing system that apply to this sequent are ldecide or udecide. However, if we use ldecide, then the premise becomes unprovable, as there is no way to remove an occurrence of ¬ ra or ¬ rb from a context that also contains ¬ q. Thus, the only possible rule will be an instance of udecide, with the focused formula in the premise being one of the Π. First, consider the case where the focused formula does not contain h, i.e., it corresponds to one of the instructions in I \ {halt}. In each of these cases, the focused phase that immediately follows is deterministic. As a characteristic case, suppose the focused formula is q⊗ ⊗ ⊗! b ¬ r; then we have:

⊢ ¬ q, q init ⊢ ? ∞ Π, E(c) \ {¬ q} , ¬ r ⊢ ? ∞ Π, E(c) \ {¬ q} , [! b ¬ r] [!] ⊢ ? ∞ Π, E(c), [q ⊗ ⊗ ⊗ ! b ¬ r] [⊗ ⊗ ⊗] ⊢ ? ∞ Π, E(c) udecide
The right premise is now itself neutral and an encoding of a different configuration. We can appeal to the inductive hypothesis to find a halting trace for it, to which we can prepend the instruction isza to get the halting trace from c. A similar argument can be used for the other instructions in I \ {halt}.

This leaves just the formulas involving h for the lowermost udecide. We cannot select any formula but q⊗ ⊗ ⊗ ¬ h from Π, for the derivation would immediately fail because h Q and there is no ¬ h in E(c) to use with [init]. So, as the formula selected is q ⊗ ⊗ ⊗ ¬ h, we have:

⊢ ¬ q, q [init] ⊢ ? ∞ Π, E(c) \ {¬ q} , ¬ h ⊢ ? ∞ Π, E(c) \ {¬ q} , [¬ h] [blur] ⊢ ? ∞ Π, E(c), q ⊗ ⊗ ⊗ ¬ h ⊗ ⊗ ⊗ ⊢ ? ∞ Π, E(c) udecide
The context of the right premise is now neutral, so the only rule that applies to it is udecide. A simple nested induction will show that sequents of this form ⊢ ? ∞ Π, E(c) \ {¬ q} , ¬ h are always derivable in the focused calculus. Therefore, the trace that corresponds to the configuration c is just the singleton halt.

Corollary 12. The derivability of MSEL Ξ sequents is recursively unsolvable. Proof. Directly from theorems 4, 9, and 11.

Conclusion and Perspectives

We have given a fairly obvious encoding of a 2RM in a suitable instance of MSEL containing a three element subexponential signature. The encoding of the 2RM halting problem is very similar to that of [START_REF] Lincoln | Decision problems for propositional linear logic[END_REF] for MAELL; the main difference is in the encoding of the isz transitions where we can directly check for emptiness of the relevant zone instead of making an additive copy of the world and checking this property in the copy. Additives are therefore not necessary for undecidability.

Yet, this conclusion is not entirely satisfactory. If MSEL Ξ can simulate Turing machines, then it can obviously simulate a theorem prover that implements a complete search procedure for MAELL. Thus, in an indirect fashion, this paper establishes that additive behaviour can be encoded using subexponentials and multiplicatives alone. It would be interesting to build this encoding of additives more directly as an embedding of MAELL-or even just MALL-in MSEL.

This work leaves open the questions of decidability of an arbitrary MSEL with a two-element signature or a one-element signature; the latter is equivalent to the decidability of MELL itself. We also conjecture that the decision problem for an arbitrary MSEL with no unbounded subexponentials is PSPACE-hard, because it is very likely possible to polynomially and soundly encode a MALL sequent in such an MSEL.

Finally, this undecidability result should be taken as a word of caution for the increasingly popular uses of SEL as a logical framework for the encodings of other systems, such as [START_REF] Nigam | A general proof system for modalities in concurrent constraint programming[END_REF][START_REF] Nigam | Specifying proof systems in linear logic with subexponentials[END_REF]. If one is to avoid encoding a decidable problem in terms of an undecidable one, subexponentials must be used very carefully.

Figure 1 :

 1 Figure 1: Inference rules for a cut-free one-sided sequent calculus formulation of SEL. Only the rules on the last line are sensitive to the signature.

Figure 2 :

 2 Figure 2: Inference rules for a focused sequent calculus formulation of SEL.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

		[N]	[blur]
		Rules , ⊥ ⊥ ⊥, &, ⊤ ⊤ ⊤ shared with the unfocused system
	⊢ Ω, [P] ⊢ Ω, P	decide

.


Sajid Ullah Butt
email: sajid-ullah.butt@ensam.eu

Jean-François Antoine

Patrick Martin

An Analytical Model for Repositioning of 6 D.O.F Fixturing System

Keywords: fixturing system, fixture design, mechanical modelling, optimal balancing, positioning error, Lagrangian formulation

Dimensional errors of the parts from a part family cause the initial misplacement of the workpiece on the fixture affecting the final product quality. Even if the part is positioned correctly, the external machining forces and clamping load cause the part to deviate from its position. This deviation depends on the external load and the fixture stiffness.

In this article, a comprehensive analytical model of a 3-2-1 fixturing system is proposed, consisting of a kinematic and a mechanical part.

The kinematic model relocates the initially misplaced workpiece in the machine reference through the axial advancements of six locators taking all the fixturing elements to be rigid. The repositioned part then shifts again from the corrected position due to the deformation of fixturing elements under clamping and machining forces. The mechanical model calculates this displacement of the part considering the locators and clamps to be elastic. The rigid cuboid baseplate, used to precisely relocate the workpiece, is also considered elastic at the interface with the locators. Using small displacement hypothesis with zero friction at the contact points, Lagrangian formulation enables us to calculate the rigid body displacement of the workpiece, deformation of each locator, as well as the stiffness matrix and mechanical behavior of the fixturing system. This displacement of the workpiece is then finally compensated by the advancement of the six axial locators calculated through the kinematic model.

Introduction

Fixtures are used to support, locate and hold a workpiece at a desired orientation in machine space during manufacturing [START_REF] Ryll | Towards an intelligent fixturing system with rapid reconfiguration and part positioning[END_REF][START_REF] Kang | Computer-Aided fixture planning: A review[END_REF][START_REF] Cecil | Computer-Aided fixture design -a review and future trends[END_REF]. The quality of a part is influenced by the capability of a fixture to secure and locate it on the machine considering different functional conditions during fabrication. The design of fixtures is important to precisely hold the workpiece in place and compensate the errors during machining or assembling operation to ensure high product quality.

The need for high quality production, at low cost, has accelerated the research efforts in fixture design aiming at producing cost effective products without compromising on quality. To cope with current market demand, Ryll et al. [START_REF] Ryll | Towards an intelligent fixturing system with rapid reconfiguration and part positioning[END_REF] emphasize on the need for "intelligent" fixtures which should be capable of self-configuring; reducing and compensating dimensional errors; providing stability and adapting clamping forces to guarantee optimum performance.

During production, two rough parts from the same part family can have small dimensional variations before machining operation. The rough workpiece, placed on the locators, may not be completely included in the required position due to these geometrical variations, which cause the workpiece to be wasted. To avoid the time and the material loss, it is necessary to place each new part precisely. This placement necessitates the mobilisation of machine tool which is only possible on the machines having large number of DOF and high number of geometric transformations. If a processing operation (machining or assembling) only requires a low number of axes or only small displacements, the choice of a 5-axis machine is not an economically feasible one.

Generally, the fixture design involves four main steps [START_REF] Kang | Geometric and kinetic model based Computer-Aided fixture design verification[END_REF][START_REF] Wang | Computer aided fixture design: Recent research and trends[END_REF][START_REF] Boyle | A review and analysis of current computer-aided fixture design approaches[END_REF] which are; setup planning, fixture planning, unit design and verification. This work focuses on the fixture planning which ensures the precise placement of the workpiece with respect to machine tool reference. To do so, all possible errors should be eliminated or compensated. The main causes of machining errors, which in turn cause misalignment of the workpiece, are the following: Generally, the fixture design involves four main steps [START_REF] Kang | Geometric and kinetic model based Computer-Aided fixture design verification[END_REF][START_REF] Wang | Computer aided fixture design: Recent research and trends[END_REF][START_REF] Boyle | A review and analysis of current computer-aided fixture design approaches[END_REF] which are; setup planning, fixture planning, unit design and verification. This work focuses on the fixture planning which assures the precise placement of the workpiece with respect to machine tool reference. To do so, all possible errors should be eliminated or compensated. The main causes of machining errors, which in turn cause misalignment of the workpiece, are the following: The first stage of the fixture design is the choice of fixturing configuration. Somashekar [START_REF] Somashekar | Fixturing features selection in feature-based systems[END_REF] proposed a model to select the primary, secondary and tertiary planes of the workpiece, and the number of fixturing elements based on the moment and force acting on the workpiece. Menassa & Devries [START_REF] Menassa | Optimization methods applied to selecting support positions in fixture design[END_REF] determined the secondary and tertiary planes with their respective locator positions considering the primary locating details to be known. Roy & Liao [START_REF] Roy | Fixturing analysis for stability consideration in an automated fixture design system[END_REF] relocated the supports on a 3-2-1 fixturing system considering the stability of the workpiece by applying a virtual wrench. System stability is enhanced by increasing the area of the triangle formed by three supports on the primary plane. Li & Melkote [START_REF] Li | Improved workpiece location accuracy through fixture layout optimization[END_REF] presented a model to improve workpiece location accuracy by fixture layout optimization taking fixture-workpiece contact deformation to be elastic.

Local geometrical defects of rough workpiece also cause its dislocation from the desired position. A bed of nails type fixture [START_REF] Aoyama | Development of fixture devices for thin and compliant workpieces[END_REF] can be used to hold thin or soft workpieces with reduced deformation. Small Displacement Torsor (SDT) [START_REF] Clement | A study of Optimal-Criteria identification based on the Small-Displacement screw model[END_REF][START_REF] Bourdet | Logiciels des machines à mesurer tridimensionnelles[END_REF] approach is used to find the localization error of the workpiece.

Using the same approach, the geometrical deviation of machining surface relative to its nominal position is presented in [START_REF] Villeneuve | Tolerancing for manufacturing: a three-dimensional model[END_REF]. Asante [START_REF] Asante | A small displacement torsor model for tolerance analysis in a workpiece-fixture assembly[END_REF] proposed a model using HTM and SD to calculate the positioning errors of the workpiece on a 3-2-1 fixture. Overall positioning error is taken as the sum of workpiece geometric error, locator geometric error and clamping error.

The workpiece shifts from its original position under clamping and machining forces due to the deformation of elastic fixturing elements. Clamping and machining forces are applied once the workpiece is initially located in the fixture. The effect of clamping on the workpiece displacement and the optimized or minimum clamping forces models are presented in [START_REF] Li | Analysis of reactions and minimum clamping force for machining fixtures with large contact areas[END_REF][START_REF] Li | Fixture clamping force optimisation and its impact on workpiece location accuracy[END_REF][START_REF] Deng | Determination of minimum clamping forces for dynamically stable fixturing[END_REF]. Jayaram et al. [START_REF] Jayaram | A fast analytical method to compute optimum stiffness of fixturing locators[END_REF] calculated the minimum stiffness of each locator required to withstand the applied load with acceptable workpiece displacement considering the workpiece to be rigid. Raghu & Melkote [START_REF] Raghu | Modeling of workpiece location error due to fixture geometric error and fixture-workpiece compliance[END_REF] predicted the final position and orientation of the workpiece due to fixture geometric errors using the part loading chart [START_REF] Raghu | Analysis of the effects of fixture clamping sequence on part location errors[END_REF] and fixture workpiece compliance [START_REF] Li | Fixture clamping force optimisation and its impact on workpiece location accuracy[END_REF].

Significant work has been performed on the determination of the machine tool error [START_REF] Lin | A Generic Kinematic Error Model for Machine Tools[END_REF][START_REF] Jha | Analysis of geometric errors associated with five-axis machining centre in improving the quality of cam profile[END_REF][START_REF] Wan | A unified framework of error evaluation and adjustment in machining[END_REF][START_REF] Martin | Virtual manufacturing: prediction of work piece geometric quality by considering machine and set-up accuracy[END_REF][START_REF] Zhu | Integrated geometric error modeling, identification and compensation of CNC machine tools[END_REF]. It is impossible to eliminate machine tool error due to structural and production limitations, therefore it is necessary to compensate them. The compensation can be achieved either by changing tool path, moving the cutting tool, or moving the workpiece in the machine coordinate system. The easiest error compensation method is changing the cutting tool path using the NC part program [START_REF] Ramesh | Error compensation in machine tools -a review: Part i: geometric, cutting-force induced and fixture-dependent errors[END_REF][START_REF] Jha | Analysis of geometric errors associated with five-axis machining centre in improving the quality of cam profile[END_REF][START_REF] Raksiri | Geometric and force errors compensation in a 3-axis CNC milling machine[END_REF][START_REF] Zhu | Integrated geometric error modeling, identification and compensation of CNC machine tools[END_REF] as shown in figure 1, but this compensation requires 4 or 5-axis machine tools to perform the necessary transformation. In this article, a fixturing system is proposed which can hold complex workpieces and which is capable of performing a 6-DOF repositioning in the machine coordinates without the use of a 4 or 5-axis machine. The proposed repositioning principle is shown in figure 2, where a high quality baseplate is introduced in between the locators and the workpiece. The workpiece is assumed to be fixed rigidly on the baseplate and any modification of its position can be performed through repositioning of the baseplate on locators. Baseplate repositioning will eliminate the position uncertainty due to local geometric errors and avoid the use of high DOF machines at each machine center; instead, a single fixture can perform the workpiece repositioning at each machine center. The proposed fixturing system is capable of; 1. 1. Determining the relative positioning error (Due to geometrical defects or due to deformation under load) between the workpiece and the tool before and during machining or between the two parts during assembling 2. 2. Ensuring the axial displacement of 6-locators in order to reorient the workpiece at an optimal position

The proposed system aims to perform automatic on-line or off-line workpiece repositioning operation. It ensures the pre-positioning of complex parts for precise machining operations. The system can also be used on the automatic production lines where the number of axis is limited for each station.

The proposed system allows better positioning of the workpiece on the fixture, hence limiting allowances. In this article, the kinematic and mechanical model of the proposed fixturing system is presented.

The article is composed as follows; section 2 presents the kinematic model of the reconfigurable fixturing system, by considering large displacement transformation for the repositioning of the workpiece. The kinematic model is illustrated through CATIA simulation. Section 3 discusses the mechanical model of the fixturing system assuming the locators to be non rigid elements. The assembly is clamped while machining forces are introduced. A case study demonstrates the mechanical model using the Lagrangian formulation and small displacement transformation (SDT). This case study identifies the vibrational attributes of the fixturing system, while taking into consideration the stiffness of each locator and clamp, as well as the mass of the baseplateworkpiece assembly. Assuming that unknown initial position could imply large displacements (LD) of workpiece during correction phase, the kinematic model is built using homogeneous transformation matrices (HTM) and LD formulation. The geometrical properties of the baseplate enable us to calculate the position of the baseplate from the positions of six locators. The position of the workpiece is obtained as per the hypothesis of the rigid workpiece-baseplate contact. This allows the kinematic model to be more efficient and repeatable. In addition, it would be easier to manoeuvre the workpiece-baseplate assembly through the locators. In the proposed fixturing system, the workpiece is located and clamped on the baseplate forming a single rigid assembly.

Kinematic model of the fixturing system

Formalization

It is assumed that the positioning error of the baseplate is negligible as compared to the positioning error of the workpiece. For the workpiece repositioning, initial position of the workpiece can be measured through CMM and the whole workpiece-baseplate assembly can then be placed on the fabrication machine. The initial position of the workpiece is compared with its required position. Repositioning is necessary if the difference between initial and final position is beyond limits. There needs to be a proper mathematical formulation to perform the transformation of the workpiece from its initial to final position with the help of six locators which are able to move only axially.

Positioning transformation scheme of the proposed fixturing system is shown in figure 4. Here, X i represents the position vector of reference i while [P ij ]

[P bP ] = [P Ob ] -1 [P OP ] (1)
The final position of the workpiece with respect to machine coordinate (X F )

is the required position of the workpiece and is known through part program.

The reorientation of the baseplate ([P bb ′ ]) enables us to reorient the workpiece ([P P F ]) as the contact between the baseplate and workpiece is supposed to be unchanged by the load modifications ([P bP ] ≡ [P b ′ F ]). As the initial and final positions (X P and X F ) of the workpiece are known, we get the HTM equations as,

[P OF ] = [P Ob ′ ][P b ′ F ] = [P Ob ′ ][P bP ] (2)
[

P Ob ′ ] = [P OF ][P -1 bP ]
Equations 1 and 2 give the final HTM of the baseplate in the machine coordinate system which can be written in terms of locators' positions in HTM form as,

[P Ob ′ ] =          a ′ 3 a ′ 2 a ′ 1 x ′ b b ′ 3 b ′ 2 b ′ 1 y ′ b c ′ 3 c ′ 2 c ′ 1 z ′ b 0 0 0 1          = [P OF ].([P Ob ] -1 .[P OP ]) -1 (3)
where, a ′ i , b ′ i and c ′ i are the components of unit normal vectors of each plane and [P Ob ′ ] is the absolute HTM of the baseplate with respect to machine coordinate system needed to reorient the workpiece at the position required by the part program. Due to no friction hypothesis, the final calculated position of each locator will not be in line with its axis, so that repositioning is not possible. To solve this issue, the cuboid formed by the actual locators positions is centred with the cuboid formed by the final calculated positions. For this purpose, the contacting points' positions along the axis of each locator are calculated on respective surface by using the equation of each surface. For example z position is fixed for locator 1. The final plane equation of primary surface can be calculated as,

a ′ 1 x ′ 1 + b ′ 1 y ′ 1 + c ′ 1 z ′ 1 = D 1 (4)
where, D 1 is the vertical distance of the plane from machine origin. The point z * on the plane, which corresponds to the fixed lateral position of locator 1(x 1 , y 1 ), can be found by,

z * 1 = D 1 -a ′ 1 x 1 -b ′ 1 y 1 c ′ 1 (5)
Similarly, the axial displacement of all the locators can be calculated to center the cuboid formed by the locators with the one required to reposition the workpiece.

Case Study

In order to validate the kinematic model, a case study is performed on a hip prosthesis repositioning through CATIA simulation. A CPT 12/14 Hip Prosthesis by Zimmer [START_REF] Zimmer | zimmer hip prosthesis, CPT 12/14 cemented stems[END_REF] is chosen as a demonstrative workpiece. This workpiece is fabricated in a single unit, by employing a precise cost effective machining process.

Prosthesis replacement of large human joints is one of the most promising methods in treating post-traumatic and degenerative dystrophic joint diseases.

In 2006, the annual number of prosthesis replacements of hip joints was 300K in the USA, 60K in Germany and 20K in Russia [START_REF] Rosenberg | Trends and developments in the manufacturing of hip joints: an overview[END_REF]. Due to the nature of custom design process, the CNC milling fabrication process is a suitable choice using a milling operation [START_REF] Dietrich | Designing and manufac turing of customized human bone endoprostheses[END_REF]. Here, it is assumed that the milling operation is performed initially on one half of the workpiece and then the workpiece is inverted to perform the machining in the remaining half part.

The hip prosthesis is created in CATIA and its original dimensions are slightly increased and supports are added to obtain a rough workpiece before machining. It is supposed that this workpiece is clamped rigidly on the baseplate which is further located through six rigid locators. An inverse impression of the workpiece (like a half die) is created with the original hip prosthesis dimensions and is placed on a fixed position with reference to the machine origin. This position represents the tool path on the machine as the tool moves with reference to machine and not with reference to workpiece. A boolean operation is performed to simulate the machining operation by subtracting the common material from the workpiece. Two slots are made during machining of the first half part which will help to place the workpiece on two well positioned blocks after inverting.

Data input

The analytical model is implemented in a worksheet directly linked to the CATIA model which furnishes the initial position of roughly placed workpiece(P OP of equation 1). This position should be obtained by CMM in real environment as presented in table 1 along with the initial positions of all the locators. The initial position of the baseplate (P Ob ) is a function of locators' positions. The machining performed on this initially roughly placed workpiece is shown in figure 5. The workpiece should be repositioned at the required position to perform a precise machining operation. This final position (P OF of equation 2) is known by the part program and is shown in the table 2.

Results

The algorithm calculates the final locators' positions (P Ob ′ ) to reorient the workpiece at therequired position. This final calculated position of each locator is shown in table 3(a) by rounding off the values to two decimals representing the locators' advancement precision of 1µm. These calculated values are again introduced as input to CATIA model to check the final attained po- The workpiece is then inverted and clamped on the baseplate. Two slots, created during the machining of the first side, help to place the workpiece on two precisely placed blocks. This helps to place the workpiece precisely with respect to the first machined surface. Same procedure is performed for the repositioning of the second side of the prosthesis on the baseplate. The final assembly and end-product after the repositioning and machining simulation on the second face of the prosthesis is shown in figure 7, and the error with respect to required position of the prosthesis is shown in table 4.

Figure 7: Machining performed on the workpiece after repositioning the second side and the final product The above example validates the algorithm. For the second face repositioning, angular error is almost zero because of input data being the output of the first correction, which was in the precision range. The procedure does not change with the complexity of the workpiece, therefore, the same repositioning procedure can be used for more complex parts.

Robustness of the model

The workpiece position uncertainty can be calculated from the Plücker coordinates [START_REF] Halevi | Principles of Process Planning: a Logical approach[END_REF] as the function of locators' advancements precision. In our case, using the locators' input positions (table 1), the uncertainty at reference point P (table 2) is a function of six advancements,

                               δx P δy P δz P δα δβ δγ                                =                               
                               (6)
where, dz 1 , dz 2 , dz 3 , dy 4 , dy 5 and dx 6 are uncertainties of the locators' advancements, in our case, assumed to be 0.01 mm maximum.

Mechanical model

In the previous section, a kinematic model is presented which is capable of correcting the positioning error if the initial and required positions of the workpiece are known. In reality, locators can deform under the weight of the baseplate, as well as under static and dynamic forces acting on the workpiecebaseplate assembly. The baseplate, repositioned by kinematic model, is clamped to the pallet to remain at the corrected position and undergo the machining forces. In fact, these mechanical actions will imply displacements on the baseplate-workpiece assembly, so that the previously corrected position gets modified which could result in the part being dislocated and wrongly processed.

In this section, a mechanical model of the fixturing system is formulated to calculate the displacement of the workpiece caused by the deformation of elastic fixturing elements under load. For this purpose, the overall stiffness and mass of the whole fixturing system are calculated considering the locators as being the elastic elements with negligible masses and the workpiece-baseplate assembly as a rigid mass element. The analytical formulation is used because it calculates the mass and stiffness matrices quickly for any configuration and advancements of locators. This analytical mechanical model is formulated and solved in a computational system (Mathematica ) to obtain the result. As the deformation must remain small for the system reliability and lifetime, small displacement theory may be applied. Also, it is assumed that the baseplateworkpiece contact stays rigid and unaffected under the machining forces.

For the computation, the mechanical model needs the initial values such as the positions and orientations of the locators and clamps, their stiffness matrices and mass matrix of the baseplate, magnitudes and location of external forces and moments etc., as input data. Lagrangian formulation (eq. 7) is used to calculate the mechanical behavior of the fixturing system considering small displacements. The reference transformation of the mechanical model is shown in figure 8 where, P F F * = P b * b ′ is the displacement of the workpiece-baseplate assembly under load. Once this displacement is known, The workpiece can be reoriented as detailed in the previous section.

∂ ∂t ∂(T -U ) ∂ qi - ∂(T -U ) ∂q i = ∂W ∂q i (7)
Figure 8: Transformation of reference axes for mechanical model

Energy calculation

Total energy of a system is composed of its potential and kinetic energies.

Kinetic energy is negligible for the locators due to their negligible mass while potential energy of all locators can be calculated using equation 8.

U = 1 2 {∆X} T i [K] i {∆X} i (8)
where, [K] i is the stiffness matrix of i th locator and ∆X i is the relative displacement vector of the contact point of that locator relative its other extremity. Total kinetic energy of workpiece-baseplate assembly, consisting of translational and rotational kinetic energy [START_REF] Lalanne | Mécanique des vibrations linéaires (avec exercices corrigés et programmes de calcul)[END_REF], can be calculated by equation 9.

T = 1 2 { V } T [M ]{ V } + 1 2 { Ω} T [I]{ Ω} (9)
where, { V } and { Ω} are the translational and rotational velocities respectively, while [M ] and [I] are mass and inertia matrices. Kinetic energy terms can be neglected for this specific case.

Clamping forces

Clamp is the means of tightening the workpiece on the fixture once it is located. As the locators are assumed to be elastic, the clamping force causes the locators to compress. Here, it is to be noted that if the initial position of the workpiece is measured through CMM, the deformation of the locators due to the load of the baseplate-workpiece assembly is already taken into account, so only the external forces will cause further deformation. Normally, the clamp is brought in contact with the baseplate and then it is tightened. Instead of measuring the magnitude of clamping force exerted by each clamp, the amount of distance at which the external end of the clamp is moved, is used to model the clamp. The clamp is taken as elastic element; a part of the external displacement will be transferred to the other end of the clamp (baseplate-clamp contact) while the rest will induce potential energy in the clamp. Here, transformed displacement vectors ({∆X} C,i ) are calculated for each clamp as the functions of workpiece displacement vector {∆X P , ∆Y P , ∆Z P , ∆β, ∆γ, ∆α} T .

The total potential energy contained in all the clamps can be calculated using equation 10.

U C = 1 2 {∆X} T C,i [K] C,i {∆X} C,i (10)

Machining forces

Clamps behave like static forces while the machining force changes continuously with time which makes it very difficult to compensate. At the moment, the proposed mechanical model cannot compensate the machining forces in real time but in processes such as drilling, boring, assembling etc., it is relatively practical. In drilling, for example, a force acts at the center of the hole along with a moment couple. A generalized equation is proposed for total work done caused by external forces and moments in equation 11, where, {F } is the applied force vector, {T} is the external torque vector, {∆X P } is the displacement of point P of the baseplate under applied load, while {∆Θ} is the angular displacement due to applied torque. If an external force is acting away from the center of gravity, it will cause linear displacement as well as angular displacement of the workpiece; in that case, a HTM is used with the vector {∆X P } which adds the angular displacements to the work done by the force. The generalized equation for the work done can be written as,

W = {F } i .{∆X P } + {T} i .{∆Θ} (11)
Here, it is important to note that clamping force is the potential energy contained in clamps, so they have to be used on the right hand side of Lagrangian equation (eq. 7) with potential energy, while machining forces and load are external forces, so they have to be used on the left hand side of the Lagrangian equation (eq. 7).

Case study

A case study is performed on the fixturing system to explain the working and findings of the proposed mechanical model. A fixturing system, having six locators and two clamps, considered as three dimensional springs, is shown in figure 9. {F } is taken as machining force at any point on the workpiece, {T} is the moment of cutting tool, [K E ] 1 and [K E ] 2 are the stiffness matrices 5 where, X p , Y p and Z p are the positions of point "P " while X i , Y i and Z i are the positions of point "i" in x, y and z directions respectively with i=1,2,...6. Stiffness matrix for first locator, as well as mass and inertia matrices for the baseplate are taken as,

[K] i =         k r 0 0 0 k r 0 0 0 k a         , [M ] =         m 0 0 0 m 0 0 0 m         [I] =         65 48 mL 2 0 0 0 65 48 mL 2 0 0 0 8 3 mL 2        
where, k r is the radial stiffness of the locator which is the sum of share and bending stiffness placed in series, while, k a is the axial stiffness of the locator.

Table 5: Positions of each locator with respect to the center point P (mm)

Locator X p -X i Y p -Y i Z p -Z i 1 L -L L/4 2 -L -L L/4 3 0 L L/4 4 L 2L 0 5 -L 2L 0 6 2L 0 0
A screw-nut system is shown in figure 11 which has been considered as the

-X C Y p -Y C Z p -Z C 1 0 -2L -L/4 2 -2L 0 -L/4
locator model for this case study. Stiffness matrix of each locator is transformed depending upon their orientations with respect to the first locator. k a remains unvaried with the advancement of locator while k r will change by the length of the locator projected out of the support (fig. 11). In our case, the projected length of locator, out of the support, is taken 15mm for all locators. Here, the workpiece is assumed to be small so that its mass can be neglected as compared to that of the baseplate. The baseplate is chosen to be made of steel with a density of 7850kg/m 3 , the volume can be derived from figure 10 to calculate mass of the baseplate. The values for each variable are,

{X} E1 =                0 -1.77 -1.77                × 10 -5 m {X} E2 =                -1.77 0 -1.77                × 10 -5 m
The algorithm calculates the energy contained by clamps from the applied external displacement. For this specific example, a drilling operation is chosen for which F is 74N vertically downward and moment couple is -0.8 N-m for a hole diameter of 5 mm [START_REF] Onwubolu | Response surface methodology-based approach to CNC drilling operations[END_REF]. The force is acting at x = -0.1m and y = -0.1m from point P . The weight of the baseplate is also assumed to be acting at point P .

Results

Total potential energy is the sum of potential energies of all the locators and clamps. The mass and stiffness matrices of the system are obtained from calculation carried out in Mathematica . The final stiffness matrix of the whole system is,

[K] =            101
           × 10 7 (12)
The above calculation is carried out considering only locators' body stiffness and it gives overestimated values. In fact, contact stiffness between spherical locator surface and baseplate surface is smaller than that of locator bodies.

As this stiffness is in series with the body stiffness, total system stiffness will be reduced. An iterative process is carried out to calculate the non-linear contact stiffness of each locator and their deformations. The overall stiffness of the system, calculated as the result of locators' body and contact stiffness, is shown in equation 13.

[

K] =           
           × 10 7 (13)
Here, it is to be noted that during the calculations, the baseplate rotates due to compliance of the springs, so the terms with angles appear in mass matrix calculation. These terms can be neglected as the displacements are assumed to be small. Lagrangian formulation calculates the workpiece displacement vector. The calculated final (linear and angular) displacements are shown in

M ] -1 [K] -ω 2 [I] = 0 (14)
The vibrational modes can be found by solving the equation 14 as in [START_REF] Kelly | Schaum's outline of Theory and Problems of Mechanical Vibrations[END_REF]. The calculated natural frequencies for this example, taking contact stiffness into account are,

ω =                                        6803
                                      
rad/sec

Results analysis

From the above results, it is observed that the system is most rigid in z direction (3 locators & 2 clamps) and least rigid in x direction (1 locator & 1 clamp), also angle α is negligibly small as only external couple of -0.8 Nm affecting it. The result validates the mechanical model for the example.

The mechanical model can be applied to more complex problems with multiple loads, multiple directional stiffness and different locations of clamps and locators. Also a huge difference of the stiffness of the fixturing system, by inclusion of locators' contact stiffness, is noticed. Therefore, the effect of contact stiffness on the precision of the workpiece cannot be neglected.

Natural frequencies of the system are calculated for static mean load. The speed of the spindle has to steer clear from these resonance frequencies for safe machining operation. External displacement of the clamps is increased from 0.025mm to 0.05mm and it is noted that the first resonance is shifted from 1698 rad/sec to 2725 rad/sec. This proves that the machine's safe operating frequency can be changed by tightening or loosing the clamps.

Conclusion and Perspectives

In this article, an analytical model for the machining fixture is demonstrated. A baseplate is placed in between the machine table and the workpiece to avoid the irregular contact of the workpiece and the locators. Stepwise Then in the mechanical model, the locators are assumed to be elastic and the baseplate is taken as the rigid mass element. The mechanical behavior of the fixturing system is calculated considering the body and contact stiffness of all the locators, stiffness and external displacements of clamps and machining forces. A case study is performed on a fixturing system in which a workpiece displacement vector is obtained due to the deformation of the elastic elements under external load. The advancement of each locator is also calculated to compensate this displacement. This model is also valid for more complex systems with multiple clamps, multi-directional stiffness and multiple external forces and moments.

The repositioning system has not yet been designed to compensate very fluctuating forces due to machining process: chosen actuators are not fast enough. To achieve this function in order to adjust the workpiece position following high accuracy, piezoactuators could be inserted between workpiece and baseplate for active compensation.

1 .

 1 Error due to placement of locators 2. Error due to geometrical/form defects of workpiece 3. Deformation under clamping and machining forces 4. Machine tool error or machine kinematic error, thermal error and mechanical error Here, we shall state a brief literature review of the above mentioned causes of errors.

Figure 1 :

 1 Figure 1: Compensation through NC part program

Figure 2 :

 2 Figure 2: Compensation through workpiece repositioning

Figure 3 :

 3 Figure 3: Proposed fixturing system

Figure 4 :

 4 Figure 4: Transformation of reference axes for kinematic model transformation

Figure 5 :

 5 Figure 5: Machining performed on the workpiece placed at the initial position

Figure 6 :

 6 Figure 6: Machining performed on the workpiece after repositioning

Figure 9 :

 9 Figure 9: Representation of the fixturing system as spring mass system

Figure 10 : 3 . 4 . 1

 10341 Figure 10: Positions of the locators, clamps and force

Figure 11 :

 11 Figure 11: Screw-nut controlled wedge-slope Locator system

k r = 7 .

 7 44 × 10 8 N/m , k a = 9.24 × 10 8 N/m L = 60mm , m = 13.56kgFor the calculation of clamps' potential energies, the positions of clamps are shown in table 6. Clamps are supposed to be having unidirectional stiffness (k E1 = k E2 = 9 × 10 7 N/m) and displacement of both clamps are taken as X E1 = X E2 = -0.025mm. Both the clamps are supposed to be inclined at 45 degree angles. The external displacement vectors for both clamps become,

  development is performed starting from the kinematic model in which the workpiece is rigidly fixed on the baseplate. The baseplate is located through a 3-2-1 locating configuration and all the elements of the fixture are considered to be rigid. The kinematic model calculated the locators' advancements which enabled us to relocate the workpiece indirectly by baseplate relocation.

Table 1 :

 1 Input data from the initial positions

	(a) Initial locators' positions
	(Axial positions are highlighted)
	Locator no x (mm) y (mm) z (mm)
	1	70	100	15.00
	2	180	100	15.00
	3	120	40	16.00
	4	70	10.00	40
	5	180	11.00	40
	6	10.00	60	40
	(b) Initial workpiece position
	Plane Angle Degree Point P	mm
	α i	0.52	x P	103.5
	β i	-0.96	y P	60.57
	γ i	0.01	z P	70.67

Table 2 :

 2 Chosen final position of the workpiece

	Plane Angle Degree Point P mm
	α f	0	x F	100
	β f	0	y F	60
	γ f	0	z F	70

Table 3 :

 3 Result obtained through the algorithm

	(a) Calculated Locators' positions
	(Axial positions are highlighted)
	Locator no x (mm) y (mm) z (mm)
	1	70	100	15.00
	2	180	100	15.00
	3	120	40	14.99
	4	70	9.92	40
	5	180	9.92	40
	6	6.50	60	40
		(b) Error after correction
	Plane Angle Degree Point P	mm
	α i	0.00	x P	0.000
	β i	0.01	y P	0.005
	γ i	0.00	z P	-0.003

Table 4 :

 4 Error of the second face after correction

	Plane Angle Degree Point P mm
	α i	0.00	x P	0.002
	β i	0.00	y P	0.00
	γ i	0.00	z P	0.00

Table 6 :

 6 Positions of each clamp with respect to the center point P(mm) Clamp no. X p

table 7 (

 7 a) and the advancements of each locator, required to compensate the displacement of workpiece-baseplate assembly, are shown in table 7(b). These advancements are calculated through the proposed kinematic model.

Table 7 :

 7 Displacements under static mean load

	(a) Displacement of baseplate
	Parameter Displacement
	∆x P	-10.47 µm
	∆y P	-7.98 µm
	∆z p	-5.29 µm
	∆β	-47.36 µrad
	∆γ	61.98 µrad
	∆α	0.67 µrad
	(b) Advancements of all the locator
	Locator Advancement (mm)
	1	-0.01
	2	-0.01
	3	0.00
	4	0.01
	5	0.02
	6	0.00
	[

Acknowledgements

The authors gratefully acknowledge the financial support from Higher Education Commission (HEC), Pakistan.


Binbin Xu

Stéphane Binczak

Sabir Jacquir

Oriol Pont

Hussein Yahia

Parameters Analysis of FitzHugh-Nagumo Model for a Reliable Simulation

Derived from the pioneer ionic Hodgkin-Huxley model and due to its simplicity and richness from a point view of nonlinear dynamics, the FitzHugh-Nagumo model has been one of the most successful neuron / cardiac cell model. It exists many variations of the original FHN model. Though these FHN type models help to enrich the dynamics of the FHN model. The parameters used in these models are often in biased conditions. The related results would be questionable. So, in this study, the aim is to find the parameter thresholds for one of the commonly used FHN model in order to pride a better simulation environment.

The results showed at first that inappropriate time step and integration tolerance in numerical solution of FHN model can give some biased results which would make some publications questionable. Then the thresholds of parameters α, γ and ε are presented. α controls the global dynamics of FHN. α > 0, the cell is in refractory mode; α < 0, the cell is excitable. ε controls the main morphology of the action potential generated and has a relation with the period (P = 3.065 × ε -0.8275 + 4.397). To show oscillations of relaxation with FHN, ε should be smaller than 0.0085. γ influences barely AP, it showed linear relationship with the period and duration of action potential. Globally, when |α| ≥ 0.1, ε < 0.0085, there is no definite threshold for γ, but smaller values are recommended.

Introduction

FitzHugh-Nagumo (FHN) model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF] has been since long time one of the basic models to study the cardiac / neuron dynamics, due to its simplicity (2 variables, 3 parameters) and relation to ionic models. It is also the second most cited cardiac / neuron model. The original model respects basically the conditions derived from the Hodgkin-Huxley model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. So, it can be considered as physiologically correct in certain A C C E P T E D conditions. However, there is nowadays important parameters' diversity in studies with FHN model, which are often considered them as the default / original ones. Unfortunately, this is not always true.

History became legend. Legend became myth. It has been more than 60 years since the publication of FHN model. Today's FHN model is rather than generalized, somehow lost the original conditions of its parameters. The derived models would help to understand FHN model from different aspects. But they act indeed differently. The obtained result could be sometime questionable and even misleading, when inappropriate conditions are applied.

So the aim of this study is to find the parameter thresholds for one of the commonly used FHN model. We track back to the origin of the FHN model in Section 2. In Section 3, we present at first that a problem from numerical solution of FHN model can indeed give some biased results which would make some publications questionable. Different parameter threshold and their relationship with the period of action potential are then presented.

Models 2.1 A little history of FitzHugh-Nagumo model

1928, Van der Pol model

The first model to representing the heart's dynamics is the Van der Pol model (BVP) [START_REF] Van Der Pol | The heartbeat considered as a relaxation oscillation, and an electrical model of the heart[END_REF]. In 1928, Van der Pol proposed a heartbeat equation based on the concept "relaxation-oscillations" system :

v -α 1 -v 2 v + ω 2 v = 0, (1)
where α determines the main dynamics : smaller α makes the system output be more-likely sinusoidal; larger values of α can produce oscillations of relaxation. This equation is then commonly studied as

v -α 1 -v 2 v + v = 0 by setting ω = 1.

1952, Hodgkin-Huxley model

In 1952, Hodgkin and Huxley presented in their pioneer work a four-dimensional model (V, m, h, n) of ionic mechanisms underlying current conduction and excitation of action potential (AP) in nerve [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]

= ẋ/c + x 3 /3 -x [5]
, he proposed a two-dimensional model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] of the electrical activities of nerve membrane (x corresponding to v in Eq. ( 1)) :

ẋ = c y + x -x 3 /3 + z , ẏ = -(x -a + by) , (2)
where 1

-2b/3 < a < 1, 0 < b < 1, b < c 2 .
1962, Nagumo's electrical circuit and Nagumo's equation

In 1962, Nagumo et al. confirmed this model (Eq. ( 2)) with experimental evidence from an electrical circuit model [START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF]. In this work, they proposed also a mathematical models of the nerve axon (z : normalized potential) :

∂ 3 z ∂t∂x 2 = ∂ 2 z ∂t 2 + µ(1 -z + ǫz 2 ) ∂z ∂t + z, (3)
where µ > 0, 3 16 > ǫ > 0. In fact, this transmission line simulator of nerve axon can be found in another work of theirs in 1961 [START_REF] Nagumo | Self-oscillation in a transmission line with a tunnel diode[END_REF], same time as FitzHugh's work. Later in 1965, Nagumo et al. [START_REF] Nagumo | Bistable transmission lines[END_REF] proposed a more general active transmission equation in neuron :

∂ 2 u ∂x 2 = ∂u ∂t + (u + 1)(u -m)(u -1), (4)
where 0 m > -1. It is known as Nagumo's equation. Due to their equal contribution, the Eq. ( 2) is thereby named as FitzHugh-Nagumo model (FHN).

A C C E P T E D

Today's FitzHugh-Nagumo model

The success of FHN model is not only due to its mathematical simplicity and its richness from a point of view of system dynamics, but also because of its correlation to HH model. In fact, the four-dimensional HH model could divided into two subsystems: (V, m) corresponds to v in FHN, fast dynamics representing excitability; (h, n) corresponds to w, slow dynamics representing accommodation and refractoriness. In consequence, the result obtained with FHN model can be qualitatively interpreted as other ionic models.

For these reasons, FHN model becomes thus one of the most studied models in analytical neuroscience / electrocardiology.

Many variations of FitzHugh-Nagumo model have been derived from the original one. One of the most used is :

v = v(v -α)(1 -v) -w, ẇ = ε(v -γw), (5)
where α is considered playing a dominant role in the fast dynamics (sodium currant) of the FHN model.

The condition of α is commonly chosen as α > 0 (either 0 < α < 0.5 or 0 < α < 1) to keep its qualitative electro-physiological meaning. However it is in contradiction to the one in Nagumo's equation.

So

, what is the real condition for α ?

Eq. ( 5) is in fact a variation of the combination of FitzHugh's model (Eq. ( 2)) and Nagumo's equation (Eq. ( 4)). This model takes advantages of both models : global dynamics from FitzHugh's model and more generalized cubic function from Nagumo's equation. The cubic functions from Eq. ( 2) and Eq. ( 5) can be represented respectively :

f (x) = x - x 3 3 = - x 3 x + √ 3 x - √ 3 , (6a)
f (v) = v(v -α)(1 -v) = -v (v -α) (v -1) . (6b)
Eq. ( 6) implies in a intuitive way that, to ensure that the model could produce electrical wave pulse (in other wards, excitable), α in Eq. ( 5) should be negative. However, this does not mean that positive α is wrong. In fact, α > 0 in Eq. ( 5) means that the cell simulated with FHN model is in refractory mode.

In this period, external stimulation cannot provoke action potential.

A C C E P T E D

3 Results

Integration tolerance and time step

Many studies showed that FHN model can be used to study chaotic phenomena in neuron system.

However, they are mostly based on discrete Euler method. As a first-order method, the local error of Euler method is proportional to the square of the step size. Its global error is proportional to the step size. To obtain relatively reliable results with Euler method for nonlinear stiff ODE systems like FHN, very small integration step size is only the basic condition. Even this condition is satisfied, the stability cannot be fully guaranteed. For these reasons, Euler method is not recommended for practical use [START_REF] Press | Integration of ordinary differential equations[END_REF]. Taking the same parameters and initial conditions, the first simulation is performed with larger tolerance (RelTol: 1e-3 and AbsTol: 1e-6), we obtained a arrhythmia-type v (Fig. 1, upper) which is misleading. Since the true result should be regular (Fig. 1, lower. smaller tolerance : RelTol: 1e-7 and AbsTol: 1e-10). It showed that if the stability condition is not satisfied, small numerical integration tolerance could "distort" the final results (even with high order integration methods). Especially in case of reaction-diffusion equations, the false new result could serve as stimulation to the system. If this stimulation happened in excitable mode, an action potential can be provoked. Since the new stimulations are not constant; the amplitude of newly provoked action potentials would be therefore irregular. In this way, a cardiac arrhythmia with FHN model can be generated. This phenomenon is in fact a course of

A C C E P T E D

different stimulation to the cell, instead of one real cardiac arrhythmia. The results are thus questionable.

In consequence, all other simulations presented in this study are performed with smaller error tolerance.

3.2 α > 0 vs α < 0

If α = 0, the cubic function in Eq. ( 5) will be f (v) = v 2 (1 -v). This is a special case of FHN model and lost the generality of the model. So, only α = 0 is studied. Two cases are considered : single cell and tissue (cells grid).

Single cell

All other parameters in this section are the same except for α. As shown in Fig. 2, when α > 0, since the membrane is in refractory period, external stimulation (initial condition serves as a single stimulation)

cannot induce an AP in the membrane. So the potential v decreases exponentially to resting state.

However, if α < 0, the membrane is in excitable period. Even a small initial stimulation (v ini = 0.1) can provoke a pulse train with full amplitude.

Tissue

two-dimensional cells grid. Eq. ( 5) becomes in this case :

v = D∆v + v(v -α)(1 -v) -w, ẇ = ε(v -γw),

A C C E P T E D

where D is the diffusion parameter and ∆ is the Laplace operator (4 points scheme in this study). The wave source. These two results implied same conclusion as in Fig. 2 : α > 0 in refractory period.

However, the reason of the plane wave generation in last figure when α > 0 is that border initiation is a special case of two-dimensional simulation. It is in fact due to the boundary condition of discrete ∆.

Applying a center initiation to the tissue (Fig. 4), α > 0 cannot indeed provoke a wave as in Fig. 3. For negative α, multiple circular waves are generated, which confirmed once more previous conclusion. To generate wave in case of α > 0, continuous stimulation is needed:

v = D∆v + v(v -α)(1 -v) -w + I stim
As shown in Fig. 5, a single circular wave is provoked in the tissue. The tissue is nevertheless in a refractory period. For these reasons, no other waves are generated.

This result is also important in clinical applications, for example the defibrillation by sub-threshold stimulation. To terminate fibrillation using multiple sites simulation of small amplitude, it is necessary to reduce firstly tissue's excitability, otherwise the defibrillation cannot be achieved [9]. Results in this section confirmed that if the tissue is in excitable mode (α < 0), external stimulation will become an ectopic focus. In realistic situation, collision of new waves will generate spiral wave which makes the defibrillation be unsuccessful or even worse.

A C C E P T E D

Parameters' thresholds

In this section, the aim is to determine the influence of the three parameters α, γ and ε on the action potential's cycle length (period) and duration APD90 (Action Potential Duration at 90%). So some reasonable thresholds could be obtained. After the preliminary study, the final testing ranges for γ and ε are 0.001 : 0.0025 : 0.25, α is in -0.25 : 0.0025 : -0.001. In total, 10 6 simulations have been performed.

As shown in Fig. 6, ε influences dramatically the period of AP. Since ε controls the recovery current in FHN model, so larger ε means larger w, so v will thus be reduced. In fact, ε has a power function

A C C E P T E D

which is shown in the first sub-figure of Fig. 6. Fixing ε and α, γ is linear to the period. This relationship reflects the property of first order of w from Eq. ( 5). As discussed previously, α should be negative in order to produce an AP train. However, this is just one of the conditions. Even α < 0, FHN cannot always produce valid action potential (3 rd sub-figure in Fig. 6). It depends also on ε and γ. When α > -0.1, there exists a region where no valid AP is provoked. If α < -0.1, α and period of AP have practically a linear relationship. The relationships between the duration of AP and these three parameters are almost the same as those for period. To determine their thresholds, we introduce the ratio R DP between the AP period and its duration (Fig. 7). We found if R DP > 0.33, only sinusoidal signal are observed. In this case, the FHN lost its essential oscillations of relaxation property that it is designed for. To respect this property, it should be R DP ≤ 0.33. In consequence, the strict condition for ε should be ε ≤ 0.0085 to generate oscillations of relaxation. Under this condition of ε, then α and γ can be chosen properly. It is advisable to avoid the interval -0.1 < α < 0 due to the constraint shown in Fig. 6. As for γ, no absolute threshold exists. So small values are recommended, as conventionally did.

A C C E P T E D

Conclusion

In this article, we studied the parameters ε, γ and α in FitzHugh-Nagumo model. It is clear now that α controls the global dynamics of FHN. α > 0, the cell is in refractory mode and does not respond to external stimulation; if α < 0, the cell is excitable. ε controls the main morphology of the action potential

Figure 1 .

 1 Figure 1. Different v due to numerical integration tolerance error, fourth-order Runge-Kutta method. Upper figure with RelTol: 1e-3 and AbsTol: 1e-6; lower figure with RelTol: 1e-7 and AbsTol: 1e-10. Parameters : α = 0.008 (upper figure), α = -0.008 (lower figure), γ = 0.008, ε = 0.01; initial conditions : v ini = 0.1, w ini = 0; adaptive time step.

Figure 2 .

 2 Figure 2. Comparaison of α = 0.1 (upper) and α = -0.1 (lower) in FHN model (Eq. (5)). Parameters : γ = 0.008, ε = 0.01; initial conditions : v ini = 0.1, w ini = 0; adaptive time step.

Figure 3 .

 3 Figure 3. α = 0.1 (upper) and α = -0.1 (lower), border initiation : v ini (1, 1 : 200) = 1, zeros for other cells; w ini (1 : 200, 1 : 200) = 0.

Figure 4 .

 4 Figure 4. α = 0.1 (upper) and α = -0.1 (lower), center initiation : v ini (100, 100) = 1, zeros for other cells; w ini (1 : 200, 1 : 200) = 0.

Figure 5 .

 5 Figure 5. α = -0.1 (upper) : center stimulation, I stim (100, 100) = 1

Figure 6 .

 6 Figure 6. Parameters' influence on period of AP in FHN. Red curve in first sub-figure is the fitted curve : Period = 3.065 × ε α,γ -0.8275 + 4.397. (discontinuity in third sub-figure means no valid AP generated.)

Figure 7 .

 7 Figure 7. Ratio of Period/Duration. Red line is the ratio threshold at 0.33

A C C E P T E D

generated and has a relation with the period (P = 3.065 × ε α,γ -0.8275 + 4.397). When ε > 0.0085, the obtained curves are more likely to be sinusoidal. To show oscillations of relaxation with FHN, ε should be smaller than 0.0085. γ influences barely AP, it showed linear relationship with the period and duration of AP. When α ≤ 0.1 and ε < 0.0085, γ can be freely chosen for excitable cell, but smaller values are recommended. Another issue in numerical studies of FHN model is the error tolerance and the order of integration method. It is better to avoid using low order integration methods; otherwise, the integration stability error would alter the results and lead to questionable conclusion.

A C C E P T E D