Bruno Mermet
email: bruno.mermet@univ-lehavre.fr

Gaële Simon

Side effects of agents are not just random

Side effects are an important characteristic of MAS, and proving them is an interesting issue. They often can be expressed as liveness properties. But there is no system dedicated to this kind of proof. The GDT4MAS framework allows to specify and prove the correctness of multiagent systems. This framework is mainly dedicated to prove safety properties about the system and to prove that agents achieve their goal(s). However, there is no proof principle to prove that agents satisfy liveness properties that are not part of their goal(s). In this article, we propose a proof mechanism that addresses this kind of problem: we show how we can add to GDT4MAS a proof mechanism adapted to prove leads-to properties, a subclass of liveness properties.

Introduction

During the execution of a MAS, unexpected system properties are often observed. These properties can either be useful (they can be for example called emergent properties) or harmful. In both cases, it may be interesting to prove that such properties will eventually happen in order to understand how they happen. However, to our knowledge, there is no system suitable to prove such properties.

Proving the correctness of multiagent systems is a hard problem that has been tackled for several years. Most of the works on the subject have established that a new formal specification and proof system, dedicated to multiagent systems, should be developed. Among them, GDT4MAS [START_REF] Mermet | GDT4MAS: an extension of the GDT model to specify and to verify MultiAgent Systems[END_REF] proposes interesting characteristics. Especially, the proof obligation generation process is fully automatisable and it can be applied to very large systems, essentially because it relies on theorem proving and first-order logic rather than on model-checking and propositional logic.

However, if this system is well-suited to prove invariant properties (also called safety properties) and to guarantee that agents satisfy their main goal, it does not propose any proof system to guarantee that an agent establishes liveness properties that are not part of its main goal, which is necessary when considering side effects.

So, in this article, we propose to add a new proof system to GDT4MAS dedicated to the verification of a well-known kind of liveness properties: leads-to properties.

Of course, there are techniques to verify leads-to properties in distributed systems [START_REF] Chandy | Parallel Program Design: A Foundation[END_REF][START_REF] Lamport | The syntax and semantics of tla + . Part 1: Definitions and Modules[END_REF][START_REF] Milner | A calculus of mobile processes[END_REF], but these techniques are dedicated to systems where all the processes are globally taken into account for each proof that must be performed.

Most of these techniques are dedicated to systems where processes work on independant variables, and synchronize occasionnaly to exchange information (This for instance the case of the π-calculus [START_REF] Milner | A calculus of mobile processes[END_REF]).

On the contrary, in a method such as TLA + [START_REF] Lamport | The syntax and semantics of tla + . Part 1: Definitions and Modules[END_REF], shared variables can be specified, but the proof process requires to consider, for each step of the system trace, all the actions that may be considered. This is not suitable for multiagents systems, because the number of possible actions is very large and thus, the property to prove would be to complicated to be performed by an automatic (or human) prover.

To perform efficient proofs on multi-agents systems, a compositional proof system is required. This is the case of the GDT4MAS proof system, and this is a property we require for a proof system dedicated to leads-to properties in multi-agents systems.

In the next section, we briefly introduce the notion of liveness property. In section 3, we recall the main concepts of the GDT4MAS framework. The new proof system we propose is described in section 4, and its application is examplified in section 5. Finally, a comparison with other works is proposed in section 6.

Invariant and liveness properties

When dealing with formal verification of software, many kinds of properties may be considered. In this section, we present two kinds of them: invariant properties and liveness properties.

Invariant properties

Invariant properties specify a set of states the system must satisfy at every moment. They are often presented as properties specifying that "nothing bad happens". Theses properties are mainly safety properties. Indeed, when specifying safety critical systems (as, for instance, a train control system), a first critial step is to verify that the system does not reach an unsafe state. In a formal verification system such as the B method, used by corporates developing safety critical systems, this kind of property is the only one that is formally proven.

Using temporal logic, an invariant property IP is specified as:

(IP ) However, a system doing nothing may trivially verify such invariant properties. Indeed, these properties do not specify anything about the task of the system.

Liveness properties

Contrary to invariant properties, liveness properties specify how the system should modify its state. They are often presented as properties specifying that "something good will happen". There are many kinds of liveness properties. Here are a small presentation of some of them. More details can be found for instance in [START_REF] Chandy | Parallel Program Design: A Foundation[END_REF][START_REF] Lamport | The syntax and semantics of tla + . Part 1: Definitions and Modules[END_REF].

-One-day: a given property OD will eventually become true:

♦(OD)
-Leads-to: a given property LT will eventually become true every time another property P is true:

(P → ♦LT )
-Until: a given property UP remains true until another property becomes true (and P will eventually become true):

(U P → ( (U P ∨ P ) ∧ ♦P )
In the rest of this article, we will only consider leads-to properties. Indeed, a oneday property is a subtype of a leads-to property (with P true and LT OD), and an until property is also a special kind of leads-to property.

3 The GDT4MAS framework and the GDT model

Main concepts

In the GDT4MAS framework, the MAS is described by an environment, mainly specified by variables, an invariant (denoted i E in the sequel) and a population of agents evolving in this environment. Each agent is described as an instance of an agent type. As a consequence, in the following, after a short description of the notations we use, the notions of agent type and of agent behaviour are described.

Notation

Notation 1 (primed and unprimed variable) When the value of a variable v in two execution states is considered, the value of v in the first state, called the current state, is written v, and its value in the second state is written v ′ . For instance, the action consisting in increasing the value of v by 1 is specified by the postcondition v ′ = v + 1.

Agent Type Specification

Simplified Definition 1 (Agent Type) An agent type t is mainly described by a name (name t ), a set of variables (V arI t ), an invariant (i A ), and a behaviour (b t ) defined by a GDT. Definition 1 (Goal Decomposition Tree (GDT)). A Goal Decomposition Tree describes the behaviour of the agents of a given type. Each node of this tree is a GDT goal. The tree structure is defined thanks to the decomposition of each GDT goal into subgoals using decomposition operators. A predicate called Triggering Context (TC) is associated to each GDT: an agent begins the execution of its behaviour every time its TC is true.

Simplified Definition 2 (GDT goal) A GDT goal g is described by a name (name g ), a satisfaction condition (sc g ), a gpf (gpf g ), a decomposition or an action and an ns flag (ns g ). The satisfaction condition is a predicate specifying the property the goal must establish when it succeeds, whereas the gpf (Guaranted Property in case of Failure) is a predicate specifying the property the goal must establish when its execution fails. The ns flag specifies whether the goal always succeeds (Necessarily Satsfiable or NS) or not (NNS).

Please notice that when the execution of a node fails, the invariant must still remain true. The failure of a node represents the fact that, in a real world, an agent is not always guaranted to succeed in realizing a task dealing with the environment. For instance, a robot that must move its arm may be blocked by an object, or its arm may be rusty, reducing the amplitude of its move.

Simplified Definition 3 (Action) An action α is specified by a name (name α ), a precondition (pre α ), a postcondition (post α ), an ns flag (ns α ) and a gpf (gpf α ). The precondition is a predicate specifying when the action is enabled, the postcondition specifies what that action does (x ′ = x -1 for instance expresses that the action decreases the value of x by 1), the ns flag has the value NS if the action is guaranteed to always succeed, and NNS if the action may fail. The gpf is a predicate specifying what is however guaranteed to be true if the action fails. [START_REF] Mermet | Specifying, verifying and implementing a MAS: A case study[END_REF].

Definition 2 (Goal decomposition). A GDT goal is either a leaf goal or an intermediate goal. An action is attached to a leaf goal, whereas an intermediate goal is decomposed into several subgoals, thanks to a decomposition operator. A list of decomposition operators can be found in

Among others, we can introduce the following decomposition operators:

-SeqOr: Sequential Or. It decomposes the parent goal N into several subgoals N i . Subgoals are executed from the left to the right. If the considered subgoal succeeds, N is achieved and the execution of the decomposition ends. But if it fails, the next subgoal is considered. If the last subgoal is executed and fails, the satisfaction condition of N must be evaluated to know if N is however achieved or not.

-SeqAnd: Sequential And. It decomposes the parent goal N into several subgoals N i . Subgoals are executed from the left to the right. If the considered subgoal succeeds, the next one is executed. If the last subgoal is executed and succeeds, N is achieved. But if a subgoal fails, the satisfaction condition of N must be evaluated to determine whether N is achieved or not. -SyncSeqOr and SyncSeqAnd: These operators are similar to the SeqOr and SeqAnd operator, but a subset of environment variables can be locked during the whole execution of the parent goal decomposition.

An example of GDT is given in figure 1. In this figure, goals are described by their satisfaction condition. Moreover, N S goals are surrounded by a double ellipse. In satisfaction conditions, x and x ′ respectively represent the value of the variable x before and after the goal execution.

Proof system: general principles

The proof system fo GDT4MAS relies on Proof Schemas (PS). Applying a proof schema generates Proof Obligations (PO), that may be proven by an automatic prover, such as PVS [6]. At the moment, PS allow us to prove several kinds of properties. We first prove invariants at the agent-type level and at the systemlevel. Moreover, the proof system of the method verifies that goal decompositions are valid. Most PS rely on goal contexts. These contexts are computed automatically starting from the root goal. Intuitively, the context C G of a goal G is a predicate summarizing the state in which goal G will be executed.

Notations

In this section, we present two notations of GDT4MAS that will be used in the sequel.

Notation 2 (Priming) Let f be a predicate/expression. If f contains at least one primed variable, then pr(f ) = f . Otherwise, pr(f ) is the predicate/expression derived from f where each unsubscripted variable is primed.

Examples:

pr((x = x 0 )) ≡ (x ′ = x 0 ) and pr((x = x ′ )) ≡ (x = x ′ ).
Notation 3 (Invariant) Let A an agent situated in an environment E. We write:

i A the invariant regarding variables of the agent; -i E the invariant associated to the environment variables; -i EA the conjunction of i A and i E .

A new proof mechanism dedicated to leads-to properties

This section presents the proof mechanism we propose to verify that some leadsto properties are established by an agent. Here, we only consider leads-to properties that are associated to an agent; no other agent is required to establish this property.

In this section, we consider a leads-to property L defined so:

L ≡ (P L → ♦Q L )
A classical way to establish that a leads-to property is verified by a specification consists in associating a variant and a witness to this property [START_REF] Chandy | Parallel Program Design: A Foundation[END_REF].

Informally, a variant expresses the progress towards the establishment of Q L . If it is proven that an agent makes a variant decrease and that when this variant reaches its lower bound, Q L is established, then the leads-to property L is proven. A witness is a property that represents the fact that P L has been true, and thus, that Q L must be established. In this article we propose to adapt this mechanism to verify leads-to properties of agents.

Definitions and notations

We begin by a formal definition of a variant: Definition 3 (Variant). A variant is a decreasing sequence defined in a wellfounded structure.

Of course, this definition requires to define what a well-founded stucture is. Definition 4 (Well-founded Structure). A well-founded structure (S, <) is a set S with an order relation < such that every decreasing sequence in S has a lower bound. For instance, (N, <) is a well-founded structure.

Corollary 1. A variant has a lower bound. We write V 0 L the lower bound of a variant V L .

In this article, we will only consider variants defined on (N, <). This property must be added to the invariant of the agent. Definition 5 (Witness). Let L ≡ (P L → ♦Q L ) a leads-to property. A witness is a property that must be true when P L is true, and remains true until Q L is true. Notation 4 (Variant and Witness) Let L a leads-to property associated to an agent A. We write:

-V L the variant we associate to L to prove it; -V 0 L the lower bound of the variant V L ; -W L the witness we associate to L to prove it.

Sketch of the proof process

Thanks to the variant and witness we associate to a leads-to property, proving that an agent establishes a leads-to property L consists in proving that:

1. The chosen variant is a variant:

when it has reached its lower bound, Q L is established; -once P L has been true and until Q L becomes true, the agent must execute its gdt. -there is no other agent that increases the variant; 2. The chosen witness is a witness; -it is true when P L is true; -when W L is true, it remains true until Q L becomes true. 3. The agent progresses: when the agent executes its gdt, it makes the variant decrease or it establishes Q L .

In the next parts of this section, we detail each of these steps.

The chosen variant is... a variant!

To prove that V L is a variant, we have to prove that, when it has reached its lower bound, the desired property is satisfied. So, we have to add the following proof obligation:

i EA ∧ (V L = V 0 L ) → Q L (1)
Moreover, we also have to prove that once P L has been true, and until Q L becomes true, the agent is activated, and thus executes its GDT. This is established by proving the following property, where T C A is the triggering context of the agent:

i EA ∧ W L ∧ ¬Q L → T C A (2)
Finally, we also have to prove that no other agent makes the variant increase once P L has been established until Q L is established. So, for each other agent A in the system, we have to check for every action α used in a leaf goal G (we recall that post and gpf of actions contain primed variables):

i EA ∧ C G ∧ W L ∧ (post α ∨ gpf α ) → pr(V L ) ≤ V L (3)
4.4 The witness property... is a witness!

As explained before, we associate to our leads-to property L a witness property W L that verify both following properties:

-Initialisation : W L must be true when P L is true;

The property that must be verified is the following:

i EA ∧ P L → W L (4)
-Finalization : W L remains true until Q L becomes true.

For each agent, we have to establish that, when it modifies the environment (that is to say, when it performs an action, whether it succeeds or not), if the witness is true before the action, then it is still true after the action has been performed, unless Q L has become true. So, for each action α associated to a leaf goal G of each agent A, we have to verify:

i EA ∧ C G ∧ W L ∧ (post α ∨ gpf α ) → pr(W L ∨ Q L ) (5)

The agent progresses

In order to prove that each execution of the GDT of an agent defines a progress towards the establishment of property Q L , we have to prove that the execution of the main goal performs such a progress, that is to say, the main goal of the agent is a progress goal.

Definition 6 (Progress goal (pg)). We call Progress Goal a goal that either makes the variant decrease or establishes property Q L . For a leads-to property L, we associate to each goal G a boolean pg L G that is true if and only if G is a progress goal.

Determining that a goal is a progress goal can be done by inference rules relying on the structure of the gdt, once we know which leaf goals make progress. Moreover, as the gdt execution depends on the success status of goals, we must determine, for each goal, if it is a success progress goal and if it is a failure progress goal.

Definition 7 (Success Progress Goal (spg)). We call Success Progress Goal a goal that either makes the variant decrease or establishes property Q L when it is executed and succeeds. For a leads-to property L, we associate to each goal G a boolean spg L G that is true if and only if G is a success progress goal. Definition 8 (Failure Progress Goal (fpg)). We call Failure Progress Goal a goal that either makes the variant decrease or establishes property Q L when it is executed and fails. For a leads-to property L, we associate to each goal G a boolean f pg L G that is true if and only if G is a failure progress goal.

Corollary 2. A goal is a progress goal if and only if it is a success progress goal and a failure progress goal. So, for every goal G, we have pg

L G = spg L G ∧ f pg L G .
In the following paragraphs, we first present how we determine spg and fpg leaf goals, and then, we show how we infer these properties for non-leaf goals. Finally, we give proof schemas that we have to associate to non-spg and non-fpg leaf goals.

Determining the set of spg and fpg leaf goals To determine if a goal is a spg goal, we have to check that when this goal succeeds (and so, establishes its satisfaction condition), it either makes the variant decrease or establishes Q L . Of course, we must only consider executions of this goal performed when P L has been true, which is specified by the fact that W L is true. Hence the following property that must be established by each non lazy1 spg leaf goal G:

(i EA ∧ W L ∧ C G ∧ pr(sc G )) → ((pr(V L ) < V L ∨ pr(Q L )) (6)
In the same way, a goal G is a fpg leaf goal if and only if it verifies the following property:

(i EA ∧ W L ∧ C G ∧ pr(gpf G )) → ((pr(V L ) < V L ∨ pr(Q L )) (7)
Please notice that, the gpf of an NS goal being f alse, such goals are fpg goals.

Inference of spg and fpg properties A first way to ensure that a non-leaf goal is spg or fpg consists in demonstrating that it is a consequence of the decomposition. In this article, we only detail this process for the SeqAnd/SyncSeqAnd and SeqOr/SyncSeqOr operators.

SeqAnd and SyncSeqAnd

: Let G a goal decomposed into G 1 SeqAnd G 2 .
G is a spg goal, if, in all the cases where G may succeed, the variant decreases. Goal G may succeed in three cases, detailed below:

-Of course, G succeeds when G 1 then G 2 succeed. In this case, if either G 1 or G 2 are spg goals, goal G makes the variant decrease. -Because of side effects, G may also succeed even if G 1 has failed. Then, G 1 must be fpg. -Finally, G may also succeed when goal 1 has succeeded, leading to the execution of G 2 , which has failed. In this case, if G 1 is spg or G 2 is fpg, then the variant decreases.

So, we are guaranted that goal G is spg if:

   spg L G 1 ∧ spg L G 2 f pg L G 1 spg L G 1 ∧ f pg L G 2
As a consequence, here is a sufficient condition to determine that a goal is a spg goal:

f pg L G 1 ∧ (spg L G 1 ∨ pg LG2 ) → spg L G (8)
Now, to determine if G is a fpg goal, we consider both cases where it can fail, that is to say when its first subgoal fails or when its second subgoal fails after the first one has succeded. Hence:

f pg L G 1 ∧ (spg L G 1 ∨ f pg L G 2 ) → f pg L G (9)
SeqOr and SyncSeqOr : Let G a goal decomposed into

G 1 SeqOr G 2 .
Goal G may succeed in the three following cases:

goal G 1 succeeds; -goal G 1 fails, and then, goal G 2 succeeds.

goal G 1 fails, and then, goal G 2 fails but, because of side effects, goal G succeeds anyway.

So, we have:

spg L G 1 ∧ (f pg L G 1 ∨ pg L G 2 ) → spg L G (10)
The only case where Goal G may fail is when G 1 and G 2 fail. So, the fact that one of theses goals is spg ensure that G is spg. Hence:

f pg L G 1 ∨ f pg L G 2 → f pg L G ( 11
)
Using satisfaction conditions to determine spg goals Inference rules 8 and 10 to determine if a goal is spg give sufficient properties, but theses properties are not always necessary. A typical example is when the satisfaction condition of a non-leaf goal directly establishes either Q L or makes the variant decrease. So, for every non leaf goal G that has not been characterized as a spg goal by inference rules described above, we will also verify if property 6 is true. If this is the case, goal G can be identified as an spg goal.

Non-fpg and non-spg leaf goals When a leaf goal G is not a spg goal, we however must prove that this goal does not make the variant increase when it succeeds. So, for each non-spg goal, we have to prove the following formula:

i EA ∧ W L ∧ C G ∧ pr(sc G ) → pr(V L ) ≤ V L (12)
In the same way, for each goal G that is not a fpg goal, we must prove:

i EA ∧ W L ∧ C G ∧ gpf G → pr(V L ) ≤ V L (13)
Indeed, this is necessary to guarantee that between two steps during which the agent makes the variant decrease, it is not increased in another way.

We choose here of course a very simple example, in order to be able to present all the principles of the proof. We consider a "multiagent system" with only one agent modifying the variant.

Please notice that the system may contain several other agents. In this case, as explained in section 4.3, it has to be proven that their actions do not increase the variant. Taking into account the dynamicity of the environment relies on the same principle, because, as explained in previous articles, the dynamicity of the environment can be modeled by an agent modifiying the state of the environment.

The environment contains two variables, x and d, and is specified by the following invariant:

i E =    x ∈ N d ∈ B d ↔ (x > 0 ∧ x ≤ 10) (14)
Our agent has a behaviour described by the GDT given in figure 2. In this figure, goals names (from A to E) and their simplified satisfaction conditions are given. By simplified SC, we mean that we did not write the part specifying that the value of other variables are not modified. For instance, the full SC of node D is

y ′ = 2 ∧ x ′ = x ∧ d ′ = d.
Informally, the goal of this agent is to decrease the value of the environment variable x, by 2 if possible, and otherwise by 1.

Moreover, the triggering context of the agent, its invariant and the gpf of node E are defined so:

T C a d ( 15
)
I a (y ∈ N) (16
)

gpf E x ′ = x ∧ d ′ = d (17) gpf B x ′ = x (18)
We want to prove that this agent establishes the following leads-to property:

(x = 10 → ♦x = 0) ( 19
)
We will use x as the variant and d as the witness. To conform to the notation used in the previous section, we have: This article beeing focused on the proof of liveness properties, we do not present other proofs that must be performed to guarantee the correctness of this specification.

P L (x = 10) (20) Q L (x = 0) (21) V L (x) (22) V 0 L (0) (23) W L (d) (24)
Moreover, in order to give readable formulae, we do not give full contexts of nodes and thus, hypotheses in theorems to prove are simplified.

Determining leaf progress goals

goal D As goal D is a NS goal, it is a fpg goal.

To determine if it is spg, we must establish property 6 for this goal. Thus, we have:

WL d CD d pr(scD) (y ′ = 2)
Of course, the conjunction of these properties with the invariant does not imply x ′ < x or x ′ = 0. So, D is not an spg goal. So:

spgD = f alse (25)
f pgD = true (26)
goal E When goal E is considered, we have:

WL d CE d-2 ∧ y-1 = 2 ∧ x-1 = x-2 d-1 = d-2 ∧ y = y-1 ∧ x = x-1 pr(scE) x ′ = x -y ∧ (d ′ ↔ x ′ = 0)
The context of goal E given above is calculated by the context inference rules of the GDT4MAS method. Itexpresses the fact that goal E is considered only after goal D has succeeded when it has been executed in its context.

To establish that E is a spg goal, according to 6, we must demonstrate that the conjunction of these properties imply that the variant decreases, that is to say x ′ < x. This is obvious because, from C E , we can deduce y = 2 and from pr(sc E ), we can deduce x ′ = xy. So, E is a spg goal.

We also have to determing if E is a fpg goal, thanks to rule 7. Among the hypotheses of this rule, we have gpf E (which implies x ′ = x, see 18) and requires as conclusion either x ′ < x (which cannot be true!) or x ′ = 0 which cannot be guaranted because the context does not provide any knowledge about the value of x. So, goal E is not a fpg goal.

So, we have:

spgE = true (27)
f pgE = f alse (28)
goal C About goal C, we have the following properties:

WL d CC (d-2 ∧ x-1 = x-2 ∧ x = x-1) pr(scC ) (x ′ = x -1 ∧ (d ′ ↔ x ′ = 0))
To establish that goal C is a spg goal, we must try to establish rule 6. This rule requires to prove, from the conjunction of the above properties, that the variant decreases (x ′ < x) or that property Q L is true. This is obvious because, from sc C , we deduce that x ′ = x -1, which implies x ′ < x. So, goal C is a spg goal. Moreover, as this goal is a NS goal, this is also a fpg goal. So we have:

spgC = true (29)
f pgC = true (30)
Conclusion As a conclusion, we know that no leaf goal make the variant increase. Moreover, spg goals and fpg goals are respectively the following:

SP G = {C, E} (31)
F P G = {C, D} (32)
P G = {C} (33)

Inference of the progress property

Goal B To determine if goal B is a spg goal, we apply rule 8 that provides the following sufficient condition to guarantee that goal B is spg:

f pgD ∧ (spgD ∨ pgE)
However, D is not a spg goal and E is not a pg goal. Thus, with this rule, we cannot determine that goal B is spg. So, we try to apply rule 6. Considering goal B, we have:

WL d CB d pr(scB) x ′ = x -2
And we have to establish that the conjunction of these formulae implies either x ′ < x of x ′ = 0. As sc B implies x ′ = x -2, we obviously have x ′ < x. So, B is a spg goal.

We now have to determine if goal B is a fpg goal, applying rule 9:

f pgD ∧ (spgD ∨ f pgE)
As goal E is not a fpg goal and D is not a spg goal, we can deduce that goal B is not a fpg goal. So we have:

spgB = true (34)
f pgB = f alse (35)
Goal A To determine if goal A is a spg goal, we apply rule 10, which gives:

spgB ∧ (f pgB ∨ pgC )) → spgA
As we have established before that goal B is spg (34) and that goal C is pg (33), we can establish that goal A is a spg goal.

Conclusion Goal A being a N S goal and a spg goal, we now know that each execution of the GDT of the agent makes the variant decrease.

The chosen variant is a variant

Correctness According to equation 1, to prove that the chosen variant is effectively a variant, we have to prove:

iEA ∧ x = 0 → x = 0
This is obviously true ! Activation According to equation 2, we have to prove:

iEA ∧ d ∧ ¬(x = 0) → (x = 10 ∨ d)
Once again, this formula is obviously true.

the witness is a witness

Initialisation From formula 4, weh have to verify:

x = 10 ∧ (d ↔ (x > 0 ∧ x ≤ 10)) → d
This is still an obviously true formula.

Finalization We have to apply proof schema 5 for every leaf goal (and we recall here that, according to GDT4MAS principles, the gpf of an NS action is f alse).

Goal D The NS action δ associated to goal D is defined by:

post δ y ′ = 2 ∧ d ′ = d gpf δ f alse
So, with the context of goal D given above, we must establish:

iEA ∧ d ∧ d ∧ ((d ′ = d ∧ y ′ = 2) ∨ f alse) → pr(d ∨ x = 0)
That can be simplified into:

iEA ∧ d ∧ d ′ = d ∧ y ′ = 2 → d ′ ∨ x ′ = 0
This property is obviously true (as d and d ′ = d can be found among the hypotheses).

Goal E The action η associated to goal E is defined by:

postη x ′ = x -y ∧ (d ′ ↔ x ′ = 0) gpfη x ′ = x ∧ d ′ = d
Using C E given above, applying proof schema 5, we obtain the following proof obligation:

iEA ∧ d-2 ∧ y-1 = 2 ∧ x-1 = x-2 d-1 = d-2 ∧ y = y-1 ∧ x = x-1 ∧ d ((x ′ = x -y ∧ (d ′ ↔ x ′ = 0)) ∨ (x ′ = x ∧ d ′ = d)) → pr(d ∨ x = 0)
In order to simplify the explanation of the demonstration (that can be however easily performed by an automatic prover), we remove useless hypotheses. So, we have to prove:

(d ∧ x ′ = x -y ∧ (d ′ ↔ x ′ = 0)) ∨ (d ∧ x ′ = x ∧ d ′ = d) → d ′ ∨ x ′ = 0
The structure of this formula being a∨b → c, we will successively demonstrate a → c and b → c.

-(d ∧ x ′ = x -y ∧ (d ′ ↔ x ′ = 0)) → d ′ ∨ x ′ = 0
We use a proof-by-case on the value of d ′ . Either d ′ is true, and so, the goal is true, or d ′ is false. In the latter case, according hypothesis 2, x ′ = 0, and so the goal is true. QED.

-(d ∧ x ′ = x ∧ d ′ = d) → d ′ ∨ x ′ = 0
As d and d ′ = d are hypotheses, we obviously deduce d ′ . QED.

So the proof obligation generated by applying proof schema 5 to goal E is true.

Goal C The action associated to goal C is defined by:

postγ x ′ = x -1 ∧ (d ′ ↔ x ′ = 0) gpfγ f alse
Using C C and applying proof schema 5 to goal C, we have to prove:

iEA ∧ d-2 ∧ x-1 = x-2 ∧ x = x-1 ∧ d ((x ′ = x -1 ∧ (d ′ ↔ x ′ = 0)) ∨ f alse) → pr(d ∨ x = 0)
In order to simplify the explanation of the demonstration (that can be however easily performed by an automatic prover), we remove useless hypotheses. So, we have to prove:

(d ′ ↔ x ′ = 0) → (d ′ ∨ x ′ = 0)
The proof is obvious: either d ′ is true, and so, the goal is true, or d ′ is false and so, from hypotheses, x ′ = 0 is false, and so, x ′ = 0. QED.

Conclusion

Following the proof system described in section 4, we have been able to establish that an agent whose behaviour is described by the gdt given in figure 2 satisfies a liveness property that is not a part of its main goal.

Several formal specification languages dedicated to multiagent systems exist. However, they are often not dedicated to the proof. This is for instance the case of 2apl [START_REF] Dastani | 2APL: a practical agent programming language[END_REF], that is finally more a programming language than a specification language suited to proof. MetateM [START_REF] Fisher | A survey of concurrent METATEM -the language and its applications[END_REF] gives the developer a way to specify properties, and the system controls that the execution does not violate these properties. However, this is a proof-by-construct process; this means that the proof is performed only at the execution time, and if the initial conditions change, a new proof (consisting in an execution of this new initial state) must be performed.

Finally, most works dealing with the verification of multiagent systems rely on model-checking principles. One of the most recent work in this area is the definition of AJPF [START_REF] Dennis | Model checking agent programming languages[END_REF], a model-checker relying on JPF [START_REF]Java Path Finder[END_REF] and the Agent Infrastructure Layer AIL. This is, as far as we know, the only system that proposes a way to verify leads-to properties on multi-agent systems. However, a first drawback of the method is the time taken by the system to establish the property (several hours for a very simple system). Of course, a more optimized modelchecker such as spin [START_REF] Holzmann | The Model Checker SPIN[END_REF], may greatly reduce the time required. However, such systems remain dedicated to small-size systems. Moreover, such systems have a more serious drawback: also they can be used to prove a property such as the property we have proven in section 5: (x = 10 → ♦x = 0), they cannot be applied when the left-hand side property (here, x = 10) characterize an infinite number of states. For instance, if we would be interested in proving the following leads-to property: (x ≥ 10 → ♦x = 0), a model-checking-based method would fail, whereas the process we propose would be as efficient as it is in the given example.

The same problem can be found with MCMAS [START_REF] Lomuscio | Mcmas: A model checker for the verification of multi-agent systems[END_REF], which moreover does not provide a way to verify leads-to properties. This model-checking technique tries to verify formulae specified in propositionnal logic, as AJPF. The main disadvantage of this technique is that, relying on propositionnal logic, proofs cannot be generalized on systems of any size. For instance, in the cited article, it is shown that the verification of the dinning cryptographers must be performed for each number of cryptographers we are interested in. Moreover, even if the time taken for 10 cryptographers is quite good, performances decrease dramatically when the number of cryptographers increase. Finally, with such a technique, to prove that the MAS work with any number of cryprgrapher, an infinite number of verifications must be performed, requiring, of course, an infinite time.

Indeed, as model checking techniques may be applided on systems with several millions of states, their complexity is a critical aspect that must be taken into consideration. But with theorem proving techniques, this criterion is quite less important. Indeed, each proof requires a very short time, and the number of proofs is very low, compared to the number of states generated in model checking techniques (for instance, even on a very large industrial system, less that 50,000 proofs had to be verified [START_REF] Abrial | Formal methods in industry: achievements, problems, future[END_REF]). For instance, with the GDT4MAS model, if we call n(t) the number of nodes of an agent type t and T the set of agent types, the number of proofs to perform is approximately 2Σ t∈T n(t).

In this article, we have shown that the GDT4MAS model, that was mainly dedicated to the proof of invariant properties, can be extended to prove liveness properties such as lead-to properties. As other proof obligations of the GDT4MAS framework, the new proof obligations generated are easily proven by an automatic theorem prover such as PVS.

This kind of proof can help in analyzing the behaviour of a MAS. In the work presented here, we have only considered liveness properties associated to a single agent. Of course, more general liveness properties at the system level will have to be considered, especially properties that are established not only by a single agent, but by a subset of the agents in the system. This is a shortterm perspective. Moreover, at it is classically performed in standard verification systems, our proof system can only prove leads-to properties P leadsto Q for which there is a continuous progress to Q once P has been true. In a multiagent system where agents are fully autonomous, we also have to consider properties for which this progress is not continuous. This is a long-term perspective for us.

Fig. 1 .

 1 Fig. 1. Simple GDT

Fig. 2 .

 2 Fig. 2. GDT of the example

In this article, we only focus on non lazy goals, that is to say goals that are always executed even if their satisfaction condition is already true when the goal is considered.


Olivier Goury

Pierre Kerfriden

Stéphane P A Bordas

Bridging analytical and computational homogenisation for nonlinear multiscale problems: a reduced order modelling approach for a damage problem

Keywords: model order reduction, proper orthogonal decomposition (POD), nonlinear fracture mechanics, system approximation, computational homogenisation

published or not. The documents may come

Automatised selection of load paths to construct reduced-order models in computational damage micromechanics: from dissipation-driven random selection to Bayesian optimization

Introduction

Many natural or engineered materials are heterogeneous. More precisely, they are heterogeneous at a microscopic length scale, but homogeneous at a macroscopic length scale. For a better understanding of the material behaviour, it is necessary to modelise the macroscale while taking the microscale into account. At the microscale, the material can be made of inclusions, voids, fibers, etc.

Using supercomputers, a direct numerical simulation of the material might be at reach but at tremendous computational cost, since the microscale features need to be meshed within the coarse scale material. There is hence a need for various homogenisation techniques. Those are based on the definition of a representative volume element(RVE) [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF], also called unit cell in the literature. That element has to be small enough to be considered as a point from the macroscale point of view, but large enough so that it statistically contains enough heterogeneities to be representative of the microstructure. From there, approaches can essentially be hierarchical or concurrent (or semi-concurrent), analytical of computational.

Analytical homogenisation (See for example [START_REF] Zohdi | A description of macroscopic damage through microstructural relaxation[END_REF]) assumes a certain constitutive law of the macrostructure that depends on a certain number of parameters that can be found by studying the behaviour of the RVE under certain loading conditions. This approach is limited to simple inclusion geometries, as it involves solving the RVE problem analytically, or at least guessing the form of its solution. This approach is not always applicable, especially in the case of materials with complicated geometries and behaving nonlinearly. The advantage is that obtaining the material parameters is then very fast. The other way is computational homogenisation. Unlike an analytical approach, it does not require any a priori guess of the solution of the RVE problem since it involves solving it anytime the material properties are needed [START_REF] Feyel | A multilevel finite element (fe2) to describe the response of highly non-linear structures using generalized continua[END_REF]. It is computationally more expensive since the RVE boundary value problem has to be solved repeatedly at any point of the macro-mesh where material properties are needed. 1 This paper will aim at reducing the cost of solving the RVE boundary value problem with the use of a model order reduction method. This can be seen as a bridge between analytical and computational homogenisation. Model order reduction provides a surrogate of a full order model that retains its behaviour but can be evaluated at reduced cost. Several authors have looked at ways to alleviate the cost of solving the RVE problem. This includes [START_REF] Oskay | Eigendeformation-based reduced order homogenization for failure analysis of heterogeneous materials[END_REF][START_REF] Fish | Hybrid impotent-incompatible eigenstrain based homogenization[END_REF] where low order approximations of eigenstrains were used or [START_REF] Ghosh | Micromechanical analysis and multi-scale modeling using the Voronoi cell finite element method[END_REF] with the use of the Voronoi cell finite element method.

In this paper, a projection-based method will be used. The solution will be searched for in a space of small dimension spanned by a few basis vectors. This basis will be obtained by using the popular proper orthogonal decomposition (POD) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF]. Solutions of the RVE boundary value problem when loaded with various boundary conditions (tensile, mixed mode,...) will be stored in an "offline" stage. This provides the so-called snapshot. With the use of a POD transform, few relevant modes capturing the necessary information will be selected. In the "online" stage, the solution of the RVE will be obtained from solving the RVE reduced problem. Various applications of POD can be found in [START_REF] Legresley | Investigation of non-linear projection for pod based reduced order models for aerodynamics[END_REF][START_REF] Meyer | Efficient model reduction in non-linear dynamics using the Karhunen-Loeve expansion and dual-weighted-residual methods[END_REF][START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF][START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF]. This method hence combines both features of the analytical and the computational homogenisation, making it a accurate and versatile method. Note that the concept of applying POD-based model order reduction to computational homogenisation was already proposed in [START_REF] Yvonnet | The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains[END_REF] for hyperelastic solid or in [START_REF] Monteiro | Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction[END_REF] for a conduction problem, though no system approximation was applied.

Within the POD-System approximation framework, the computations are now of the order of the reduced model, instead of the one of the full scale model. To solve the RVE boundary value problem, all that is needed now is to store the gappy operator. The approach is still computational, since the RVE boundary value problem still needs to be solved, but it is solved at very reduced costs. The "offline" stage can be thought of as a "learning" stage, where some information about the behaviour of the RVE under different scenarios of loading is stored and treated with the output being the most "informative " degrees of freedom and the corresponding gappy operator that enables a reconstruction of the missing information. In the "online" stage, the form of the solution of the RVE boundary value problem is partly restricted by the fact it needs to fit in the POD bases. This can be seen as a partly analytical method as well.

The paper is organised as follows. In section 2, we define the type of problem we want to solve and define the representative volume element equations. In section 3, we define the model order reduction approach and incorporate it in the RVE boundary value problem formulation. Section 4 presents the system approximation technique that allows a fast integration of the material properties. Section 5 presents an example and the numerical results associated. Conclusion are drawn in section 6.

Problem statement

We are considering a heterogeneous material with a damage constitutive law at the microscale. We assume that to correctly model the full strucuture, it is necessary to take the microstructure into account. Hence a multiscale approach will be used, namely a semi-concurrent approach (FE 2 type, see [START_REF] Feyel | A multilevel finite element (fe2) to describe the response of highly non-linear structures using generalized continua[END_REF]). A representative volume element (RVE) will be used to numerically evaluate the constitutive law of the macrostructure. In this paper, the RVE will be modelled with a lattice structure, made of damageable beams with various Young's moduli materialising heterogeneities (inclusions in this case). This model is chosen for its simplicity and may not be very realistic, but this is just a tool for testing the reduced order model and is thus not the focus of this paper. An illustration of the RVE is displayed in Figure 1.

RVE boundary value problem

Over the RVE, the following balance equation is verified:

Z Ω σ m (δu):✏ m (δu)dΩ = Z ∂Ω t.δudΓ , (1)

Matrix

Inclusions where σ m is the Cauchy stress at the micro-level and ✏ m the micro-strain. Unlike the macro-problem, the constitutive law is know. In this paper, we consider a model based on standard damage mechanics applied to a lattice structure made of beams. The model is described in [START_REF] Kerfriden | Statistical extraction of process zones and representative subspaces in fracture of random composites[END_REF].

Scale coupling

The Hill-Mandel principle states that the energy should be conserved from the micro to the macro scale:

σ M : ✏ M = Z Ω σ m : ✏ m dΩ (2)
The strain at the macro-level is taken as the average of micro-strain over the RVE:

✏ M = 1 |Ω| Z Ω ✏ m dΩ (3)
Equations 2 and 3 imply the following relationship between microstress and macrostress:

σ M = 1 |Ω| Z Ω σ m dΩ (4)

Boundary conditions

The micro problem is defined with a damageable lattice model. The displacement over the RVE is represented as a sum of a coarse contribution ū given by the macro-model and a fluctuation ũ throughout the RVE at the micro-level:

u(x; µ)=ū (x; µ)+ũ (x; µ) (5)
The fluctuation is assumed to be null over the RVE domain on average, which guarantees the fact no factitious energy is created:

Z Ω ũ (x; µ)dx = 0 (6)
The coarse contribution is fixed by the macro-strain tensor: ū (t)=✏ M (t)(xx ), where x is the barycentre of the RVE. In this framework, the true unknown of the RVE problem is the fluctuation ũ .

Effective damage tensor

Effective strain Parameter space Within our model order reduction framework, the parameter will be the effective macro-strain

Macrostructure

✏ M =  ✏ xx (t) ✏ xy (t) ✏ xy (t) ✏ yy (t)
, or more precisely, all its values on the discretised time line [0,T]:

✏ M (0), ✏ M (t 1 ),...,✏ M (T ). This space of a dimension theoretically infinite, but in practice its dimension depends on the time discretisation. In any case, it is of high dimension and building a reduced model that can handle any parameter in that high-dimensional space presents a challenge.

Newton-Raphson formulation

Inserting equation (5) into equation [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF], the problem can be expressed in a more general form with the unknown being the fluctuation, and the effective macro-strain being a parameter:

F int (ũ ; ✏ M )+F ext (✏ M )=0. ( 7
)
It is a non-linear system that is solved using a Newton-Raphson algorithm. At each iteration, one has to solve the following system until convergence:

K i ∆ũ i+1 = -r i , (8)
Figure 3: Representation of the parameter space of the reduced model. It is essentially any 3dimensional curves parametrised by the time. This space is theoretically of infinite dimension. Some sample loading paths are represented.

where

K i = ∂Fint ∂ũ (ũ i ) is the tangent stiffness matrix, r i = F int (ũ i
)+F ext is the residual and ∆ũ i+1 = ũ i+1ũ i is the variation in the fluctuation.

Model Order Reduction

Our goal is to solve the balance equations for any macro-strain ✏ M (t) t2[0,T ] at a cheap cost while keeping the accuracy of a quantity of interest to be defined. At the RVE level, we assume that for any boundary conditions applied on the RVE, the displacement field lies in a subspace of the admissible displacement space that is of small dimension, and can be span by a few basis vectors.

Projection-based model order reduction

The displacement fluctuation on the RVE is searched for in a vector space of small dimension U MOR = span ((φ i ) i=1,N ):

u(x; ✏ M )=ū (x; ✏ M )+ũ (x; ✏ M )=ū (x; ✏ M )+ N X 1 φ i α i = ū (x; ✏ M )+Φ ↵ (9)
The unknown of the RVE boundary value problem is now the state variable α rather than the displacement fluctuation ũ . Plugging this into the i th Newton iteration (8) leads to:

K i Φ∆↵ i+1 = -r i , (10)
The system is now under-determined, since the fluctuation and hence its variation is now restricted to a reduced space. Instead, at each step of a Newton algorithm, the state variable variation ∆↵ has to minimise the expression: kK i Φ∆↵ i+1 + r i k. In a Galerkin framework, this becomes:

Φ T K i Φ∆↵ i+1 = -Φ T r i ( 11
)
This system is of much smaller dimension than the original equation [START_REF] Legresley | Investigation of non-linear projection for pod based reduced order models for aerodynamics[END_REF].

Proper Orthogonal Decomposition

Principle

The proper orthogonal decomposition (POD) is a popular method for generating a set of orthogonal vectors representing well a high-dimensional space. It is used to define a hierarchical basis that best represent a set of data for a given amount of modes. In our case, it provides an orthonormal set of basis vectors [φ 1 , φ 2 ,...]=Φ that minimises the following cost function for any number n POD of modes:

J h.,.i (Φ)= Z µ X t ku ex (t; µ) - nPOD X k φ k . hφ k , u ex (t; µ)ik 2 , ( 12
)
where h., .i is a scalar product that remains to be defined and k.k = p h., .i is the induced norm. The parameter µ is should be written µ ({⌧ |⌧  t}) as it varies with t but we denote it µ for simplifying the notations. The cost function evaluates the sum of the differences between the exact solution u ex (t; µ) and its projection onto the basis Φ. The basis Φ is hierarchical in the sense that the vector φ 1 gives the best1 approximation of the solution of the RVE problem with 1 vector, [φ 1 , φ 2 ] the best approximation of the solution of the RVE problem with 2 vectors and so on.

This problem requires the exact solution over the entire parameter domain and can't hence be solved. However, it can be approached by calculating a number of solutions for certain values of the parameter. This is presented in the following section.

Snapshot POD

Let us now assume that, for a discrete set of values P s of the parameter µ at some specific timesteps, the exact solution is available. The set of those solutions is called the snapshot. We define the snapshot matrix S = ⇥ s 1 ,...,s nµ ⇤ =(u ex (t; µ)) t2[0,T ],µ2P s . The discrete version of the cost function ( 12) defined over P s reads:

J s h.i (Φ)= X i ks i - nPOD X k φ k . hφ k , s i ik 2 . ( 13
)
The minimisation problem associated to this cost function is:

8 < : min φ1,...,φl J s h.i (φ 1 ,...,φ l ) hφ i , φ j i = δ ij . ( 14
)
This problem can be solved analytically by defining the Lagrangian:

L(φ 1 ,...,φ l , µ)=J s h.i K 0 (φ 1 ,...,φ l )+ l X i,j λ ij ⌦ φ i , φ j ↵ -δ ij . ( 15
)

Canonical scalar product

The next step is to define a scalar product h.i. One common choice is to chose the canonical scalar product, that is: hx, yi = x T y.

In this case, the necessary optimality conditions read:

∂L ∂µ ij (Φ, µ)=φ i T φ j -δ ij =02 R, 8 i, j < l (16) ∂L ∂φ k (Φ, µ)=r φk L(Φ, µ)=-2 X i s i ⇣ φ k T s i ⌘ + l X i=1 (λ ik + λ ki )φ i , 8 i<l. ( 17
)
One can show that this is equivalent to:

SS T φ k = λ kk φ k and λ ij =0 8i 6 = j, (18)
which is a simple symmetric eigenvalue problem. It can be solve using a singular value decomposition:

S = UΣW T , (19)
where U and V are unitary matrices and Σ a rectangular matrix containing the (hierarchically ordered) singular values of S on its diagonal. An approximation of S of rank i th order is given by U b ΣW T , where b Σ is given by the truncation Σ at order k ( only the first k singular values of S remains on the diagonal). In this decomposition, the first k columns of U forms the optimal basis of order k,i . e .

8 k, φ k = u k .
Using the orthonormality of the basis Φ, one can know determine a simple error formula the POD basis of order l:

⇣ J s h.i (φ 1 ,...,φ l ) ⌘ 2 = X i ks i - l X k=1 φ k hφ k , s i ik 2 = X i n X k=l+1 kφ k hφ k , s i ik 2 = n X k=l+1 λ kk (20)
When normalised, we obtain the following error that represents how well the POD basis of order l approximates the snapshot S:

ν(φ 1 ,...,φ l )= s P n k=l+1 λ kk P n k=1 λ kk . ( 21
)
3.2.3 POD using a structure-dependent scalar product

Using a canonical scalar product means that what is being minimised is the average error of the displacement over the entire RVE domain (this is according to the definition of the cost function J s h.i ). This may not be the best choice. In the multiscale framework of computational homogenisation, the reduced model should be able to predict the average stress on the RVE, since it is the output quantity of the RVE simulation that is then used in the macro-simulation. The energy E(t) over the RVE is hence a quantity of interest. It is evaluated as:

E(t)=u(t; µ) T K(t; µ)u(t; µ), (22)
with u(t; µ) the displacement and K(t; µ)t h es t i ffness matrix. This motivates the use of the following scalar product:

hx, yi K(t;µ) def = x T K(t; µ)y. ( 23
)
With that definition, the basis vectors φ i are now defined so that they minimise the reduced-order model error in terms of energy over the RVE, which is meaningful. However, the material non-linear behaviour and history-dependence mean that the stiffness matrix evolves during time and changes with µ. This scalar product would hence vary too, which is impractical since we are looking for a basis that is fixed in both the time and parameter space. An alternative is to pick a simplified version:

hx, yi K0 def = x T K 0 y, (24)
where K 0 is the stiffness of the structure at the initial time (which does not depend on the parameter µ). This scalar product, though symmetric since K 0 is, is not well-defined: K 0 is singular since rigid body motions create no forces. However, the basis Φ is a representation of the displacement fluctuation only (which has zero displacement on the RVE boundary), it has hence no rigid body modes contribution. Therefore, on the displacement fluctuation space, the operator h., .i K0 is indeed a well-defined scalar product. In this case, the new minimisation problem to solve is:

8 < : min φ1,...,φl J s h.i K 0 (φ 1 ,...,φ l ) hφ i , φ j i K0 = φ i T K 0 φ j = δ ij . ( 25
)
The Lagrangian becomes:

L(φ 1 ,...,φ l , µ)=J s h.i K 0 (φ 1 ,...,φ l )+ l X i,j λ ij ⇣ φ i T K 0 φ j -δ ij ⌘ . (26)
The necessary optimality conditions read:

∂L ∂µ ij (Φ, µ)=φ i T K 0 φ j -δ ij =02 R, 8 i, j < l (27) ∂L ∂φ k (Φ, µ)=r φk L(Φ, µ)=-2 X i K 0 s i ⇣ φ k T K 0 s i ⌘ + l X i=1 (λ ik + λ ki )K 0 φ i , 8 i<l. ( 28
)
One can show that this is equivalent to:

K 0 SS T K 0 φ i = λ ii K 0 φ i . ( 29
)
Taking into account the fact that φ i is a basis of the fluctuation on the RVE and hence has no component of a rigid body mode (φ i 2 Ker(K 0 ) ? ), equation (29) reduces to the eigenvalue problem:

SS T K 0 φ i = λ ii φ i . (30)
Solving that equation will provide us a set of K 0 -orthogonal vectors that best represent our snapshot space in terms of energy. It is not symmetric anymore, so a singular value decomposition can not be used to solve that problem, and we have to use a more costly eigenvalue decomposition. This is not a big concern since all those computations are done in the offline stage and are done once and for all. One can prove that the error formula (21) obtained with the canonical scalar product remains true with this special scalar product.

System Approximation

Constraining the displacement in a low-dimensional space does not provide a significant computational gain, even if the systems to be solved are of smaller dimension. This is because the material of study is nonlinear and history-dependent, and its stiffness varies not only in different areas of the material but also with time. This requires to evaluate the stiffness everywhere in the material and this at each time step of the simulation. This means that the numerical complexity remains despite the simplification on the displacement. Hence, to decrease the numerical complexity, the domain itself need to be approximated. Several authors have looked into that. Notable contributions include the Hyperreduction method [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF], the missing point estimation [START_REF] Astrid | Missing point estimation in models described by proper orthogonal decomposition[END_REF], or system approximation [START_REF] Carlberg | Efficient non-linear model reduction via a leastsquares petrov-galerkin projection and compressive tensor approximations[END_REF]. Those methods share the idea that the material properties will be evaluated only at a small set of points or elements within the material domain. They differ in the way of selecting those points and in the treatment of that reduced information. In this paper, we will use the "gappy" method, very much like in [START_REF] Carlberg | Efficient non-linear model reduction via a leastsquares petrov-galerkin projection and compressive tensor approximations[END_REF][START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF].

Gappy Method

The internal forces generated by the reduced displacement F int (Φ α) will be evaluated only in a small subset of the degrees of freedom I of the domain Ω. All the elements in contact with those degrees of freedom have to be considered. We refer to those as the controlled elements. The internal forces will then be reconstructed by writing the internal forces as a linear combination of a few basis vectors themselves (just like it was made for the displacement).

F int (Φ α) ≈ ngap X 1 ψ i β i = Ψβ, (31)
where

h ψ 1 , ••• , ψ ngap i = Ψ
is the forces basis of size n gap and β the associated scalar coefficients. The coefficients β of the expansion are found so that to minimise the norm of the difference between the linear expansion and the nonlinear term over the subset I:

argmin β ? kF int (Φ α) -Ψβ ? k P , (32)
with P being a matrix so that Pij = ( 1i fi 2 I and i = j 0 otherwise and kxk P = kP T xPk 2 . P can be written EE T with E being an extractor matrix so that E T x is the restriction of x to the set I.I f the number of points in I is identical to the number of basis vectors (ψ i ) i=1,ngap , β ? can be found by solving the equation:

E T Ψβ = E T F int (Φ α) (33)
which implies:

β =(E T Ψ) -1 E T F int (Φ α) (34)
At a newton iteration of our POD-Galerkin framework, this reduces equation [START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF] to:

Φ T Ψ(E T Ψ) -1 E T KΦ∆α = -Φ T r. ( 35
)
This can be rewritten in the form:

GE T KΦ∆α = - ⇣ GE T F int (Φ α)+Φ T F ext ⌘ , (36)
where we define the gappy operator G = Φ T Ψ(EΨ) -1 .

Remark: Note that once the "offline" stage operations are done, all that the bases Φ and Ψ are calculated, the set of control points I is selected and the gappy operator is evaluated. In the "online" stage, all that remains to do is build a system of small dimension and solve it which is computationally much cheaper.

Selection of the controlled elements

The selection of the control elements will be done using the discrete empirical interpolation method (DEIM) [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF]. This method finds a set of degrees of freedom I in a greedy manner from the internal forces basis Ψ.W eq u i c k l yd e s c r i b eh e r et h em e t h o d . At iteration j of the greedy algorithm, j -1 points have been already selected. We define the extractor E j that extracts those j selected degrees of freedom (i.e. for any vector v, E j v is a smaller vector containing only the j entries of v corresponding to the selected degrees of fredom). The residual r gap = [1,j] β j -j+1 is evaluated, where [1,j] is the matrix containing the first j vectors of the basis Ψ and j+1 is the j +1 th vectors in that basis. β is the solution of the minimisation problem

β = argmin β ? E j [1,j] β ? -E j j+1 2 . ( 37
)
The solution is easily found:

β = ⇣ E j [1,j] ⌘ -1 E j j+1 .
The greedy procedure then selects the index of the highest entry in r gap as the j +1 th control degree of freedom. This procedure essentially selects the set of degrees of freedom that maximises the conditioning of the system (33). At the end of the greedy algorithm, the number of control degrees of freedom chosen equals the number of basis vectors ( i ) ngap which makes system (33) well defined.

Example and Numerical Results

We consider a sample RVE made of damageable lattice bars. The RVE is composed of a matrix with inclusions. Bars with materialising inclusions are given a high Young's modulus compared to the matrix bars. This model is simple but the focus of this study is not the multiscale model itself but the application of model order reduction.

Snapshot selection

The parameter space, made of all possible values of the macro-strain ✏ M (t) for t ∈ [0,T], is of very high dimension.

Arbitrary regular sampling

We simplify the snapshot selection to monotonic increase of the strain with time in the following form:

✏ M (t)= t N .  ✏ xx ✏ xy ✏ xy ✏ yy . ( 38
)
The snapshots are then chosen so that the macro-strains have an equivalent "norm":

✏ 2 xx +✏ 2 yy +✏ 2 xy = C
, where C is a constant. They are equally distributed on the 3-dimensional sphere parametrised by ✏ xx , ✏ yy and ✏ xy . 100 snapshots are taken on that sphere. This is schematically represented in Figure 6.

Greedy enrichment of the snapshot space

The accuracy of the reduced model greatly depends on the snapshot space and how well it samples the parameter space. Here, the parameter space contains any load path (based on the macro-strain ✏ M ) over a certain period of time until fracture is reached. After time discretisation, the parameter space is of dimension 3 × n t , since in 2 dimensions the load can be uniaxial in the x or y direction or in shear. n t stands for the number of time steps required to reach fracture. The parameter space dimension grows with the number of timesteps and can hence be quite large. To ensure an exhaustive Figure 6: Sample loading paths taken used to generate the snapshot space. The loading paths are linear and monotonic. sampling, it is necessary to enrich the arbitrary sampling defined in the previous section. We will do so following a greedy iterative procedure.

The idea is to enrich the snapshot space with the solution u enrich obtained from a certain discrete loading path ✏ M max (t 0 ,t 1 ,...) that leads to the worst approximation by the current POD basis. This "path of worst approximation" is built up in a greedy manner, incrementally at each timestep. This is the loading path that will intuitively provide the most information to the POD basis.

The greedy procedure is as follows: assuming ✏ M max (t 0 ,t 1 ,...,t i ) is the "path of worst approximation" up to timestep t i , the procedure looks for the load increment ∆✏ i max that maximises the error between the full order model and the reduced model at time t i+1 having followed the loading path ✏ M max (t 0 ,t 1 ,...,t i )u pt ot i m et i . The loading path is then extended to time t i+1 as:

✏ M max (t i+1 )=✏ M max (t i )+∆✏ i max . (39)
The procedure iterates until reaching fracture. The maximisation process to find the the load increment ∆✏ i max is described in the following paragraphs. Check Figure 7 to see load paths samples generated by this algorithm.

A sequence of maximisation problems

We define a sequence of error functions parametrised by the current timestep t i and the value of the load up to that time ✏ M ⌧ <ti :

f err ⇣ t i , ✏ M ⌧ ≤ti ⌘ = ku ex ⇣ t i , ✏ M ⌧ ≤ti ⌘ -u r ⇣ t i , ✏ M ⌧ ≤ti ⌘ k K0 . (40)
This function simply evaluates the error between the exact solution and the reduced order solution at timestep t i . Now, at step i + 1 of the greedy procedure, we are looking for the load increment that maximises the reduced order model error. In this purpose, we solve the following maximisation problem:

:

max ∆✏ M ? f err ⇣ t i+1 , ✏ M ⌧ ≤ti , ✏ M ti + ∆✏ M ? ⌘ ∆✏ M ? = δ step , (41)
where δ step is a predefined load increment value that we keep constant during the whole greedy procedure.

∆✏ M ? =  ∆✏ xx ∆✏ xy ∆✏ xy ∆✏ yy = 0 @ ∆✏ xx ∆✏ yy ∆✏ xy 1
A is the load increment from timestep t i through to t i+1 .

Problem (41) can be rewritten to become unconstrained: max

∆✏ ? xx ,∆✏ ? yy f err ✓ t i+1 , ✏ M ⌧ ≤ti , ✏ M ti + [ ∆✏ M (∆✏ ? xx , ∆✏ ? yy ) ◆ = d f err ⇣ ∆✏ ? xx , ∆✏ ? yy ; t i+1 , ✏ M ⌧ ≤ti ⌘ , (42)
where

[ ∆✏ M (∆✏ ? xx , ∆✏ ? yy )= 0 B @ ∆✏ ? xx ∆✏ ? yy q δ 2 step -∆✏ ? xx 2 -∆✏ ? yy 2 1 C
A which guarantees [ ∆✏ M = δ step .I n t h i s case we assume ∆✏ xy ≥ 0 so we keep a positive sign in front of the square root.

Gradient descent algorithm

To solve problem (42), we use a gradient descent method (also famously known as steepest descent). This is a first order optimisation method. From an initial guess

∆✏ 0 = ✓ ∆✏ 0 xx ∆✏ 0 yy ◆
,as e q u e n c e of iterates is found using the following update:

∆✏ k+1 = ∆✏ k + γ @ d f err @ (∆✏ xx , ∆✏ yy ) (∆✏ k ; t i+1 , ✏ M ⌧ ≤ti ), (43)
where γ is a scalar whose value is found through a linesearch using a bisection algorithm. The derivative @ d ferr @(∆✏xx,∆✏yy) (∆✏ k ; t i+1 , ✏ M ⌧ ≤ti ) is evaluated numerically by taking small variations (of size ⌫) around the value ∆✏ k :

@ d f err @ (∆✏ xx , ∆✏ yy ) (∆✏ k ; t i+1 , ✏ M ⌧ ≤ti ) ≈ 1 ⌫ 2 6 6 4 d f err ✓✓ ∆✏ k xx + ⌫ ∆✏ k yy ◆ ; t i+1 , ✏ M ⌧ ≤ti ◆ -d f err ⇣ ∆✏ k ; t i+1 , ✏ M ⌧ ≤ti ⌘ d f err ✓✓ ∆✏ k xx ∆✏ k yy + ⌫ ◆ ; t i+1 , ✏ M ⌧ ≤ti ◆ -d f err ⇣ ∆✏ k ; t i+1 , ✏ M ⌧ ≤ti ⌘ 3 7 7 5 , (44)
In theory, the gradient descent algorithm would run until reaching some convergence criterion. In practice, only a few iterations are performed to reach a decent optimum, as each step involves solving the full order model several times and is hence costly.

POD basis

The basis Φ is selected using the POD method explained in section 3.2. The first few modes are displayed in Figure 8. We will test the results using various number of basis vectors. Procedures to select the optimal number of modes according to robust cross-validation procedures can be found in [START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF].

System approximation

We follow the procedure described in 4. The basis Ψ is extracted from the same snapshot space as used for the displacement basis Φ. The set of controlled elements is selected using the DEIM [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF]. The amount of vectors in the basis Ψ is chosen so that the system approximation does not increase the global error of the reduced order model. The error ν tot between the exact solution and the reduced model solution with system approximation can be decomposed in the following way (with u ex (t) the exact solution, u r (t; Φ), the reduced order solution without the system approximation using the dynamic basis Φ, and u r,sa (t; Φ, Ψ) the complete reduced order model with system approximation using the dynamic basis Φ and the static basis Ψ):

ν tot (t) 2 = ku ex (t) -u r,sa (t; Φ, Ψ)k 2 2 = k(u ex (t) -u r (t; Φ)) + (u r (t; Φ) -u r,sa (t; Φ, Ψ))k 2 2
(45)

k(u ex (t) -u r (t; Φ))k 2 2 + k(u r (t; Φ) -u r,sa (t; Φ, Ψ))k 2 2 . ( 46
)
Taking this in consideration, the basis Ψ is chosen to be the smallest that verifies the inequality:

k(u r (t; Φ) -u r,sa (t; Φ, Ψ))k 2  1 10 k(u ex (t) -u r (t; Φ))k 2 . ( 47
)
This guarantees that error generated by the system approximation remains insignificant compared to the error generated by approximating the displacement. More details can be found in [START_REF] Kerfriden | A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics[END_REF]. The location of controlled elements (which are all the elements in contact with the control degrees of freedom) is shown in Figure 9 for various basis sizes.

(a) Controlled elements with 2 "dynamic" basis vectors.

16 "static" basis vectors are needed to achieve the minimal accuracy (b) Controlled elements with 4 "dynamic" basis vectors. 28 "static" basis vectors are needed to achieve the minimal accuracy (c) Controlled elements with 10 "dynamic" basis vectors. 34 "static" basis vectors are needed to achieve the minimal accuracy (d) Controlled elements with 15 "dynamic" basis vectors. 48 "static" basis vectors are needed to achieve the minimal accuracy

Figure 9: Controlled elements selected using various basis sizes. The bigger the basis, the more controlled elements are needed. The elements tend to gather around the regions where the variation of displacement is the highest, hence where the variation of the internal forces will be high.

Numerical results

In this section, we will test the performance of the method by comparing the relative error between the "truth" solution of the RVE problem, which is the solution obtained when using the full order model, and the reduced model. We will focus on the simulation of the RVE problem only, and not consider the complete computational homogenisation framework (the macroscale problem will not be solved).

For the first test, the load path is set using the following effective strain: ✏ M (t)= t T .

 11 11 . Note that this case is not in the snapshot. First, the relative error is plotted with various numbers of POD basis vectors. See Figure 11. Figure 10: Regions of interest selected by the system approximation procedure. Those regions (circled in the Figure) are matching the areas of higher displacement found in the POD bases. This is intuitively good, since those elements have to give enough information to be able to reconstruct the internal forces over the entire domain. Those are the elements whose behaviour vary the most when changing the loading path (which is the parameter of the reduced model), hence containing the core information necessary to build up an accurate reconstruction. Several remarks can be made:

• As expected, the error decreases when the number of basis vectors increases.

• However, after about 4 vectors in the basis, the error does not decrease very much and reaches a plateau. This means that no matter how many vectors in the basis, a maximum accuracy is achieved. This can be explained by the fact that the loading path tested is not part of the Figure 11: Relative error when using the reduced model without system approximation. As the number of POD basis vectors increases, the error decreases. Using more than 4 basis vectors does not have a strong effect. The error reaches a threshold. That threshold corresponds to the projection error of the exact solution onto the snapshot space.

snapshot. The only way to decrease this residual error is to enrich the snapshot space. We define u ex,snap (t) as the projection of the exact solution onto the snapshot space. Using the same principle than equation (46), we can decompose the error further (dropping parameters for clarity):

ν tot (t) 2 = k(u ex (t)u ex,snap (t)) + (u ex,snap (t)u r (t)) + (u r (t)u r,sa (t))k 2

= k(u ex (t)u ex,snap (t))k 2 + k(u ex,snapu r )+(u ru r,sa )k 2 , (49) k(u ex,snapu r )+(u ru r,sa )k 2 can be made as small as desired by taking high dimensional bases Φ and Ψ. The residual error that remains is k(u ex (t)u ex,snap (t))k 2 , which entirely depends on the richness of the snapshot space.

The error with respect to the speedup is displayed in Figure 13. What we call speedup here is the ratio of the elapsed time of the full order simulation over the elapsed time of the reduced model. It represents how many times faster is the reduced order model compared to the full order model. It can be seen that there is a proportional relation between speedup and error: as the number of basis functions increases, the speedup and the error decrease. The user can reduce the error with the price of having a slower simulation. What makes the reduced model faster is purely the bypassing of most of the elements when computing the internal forces or the tangent stiffness. Note that the speedup is not purely equal to the ratio between controlled elements and total number of elements since the Newton-Raphson procedure requires more steps to converge in the reduced order model scheme than in the full order model. However that speedup is still significant.

Speedup Relative error

Number of basis functions increases

Conclusion and perspectives

This paper aimed at reducing the computational burden of solving the representative volume element (RVE) boundary problems arising in semi-concurrent computational homogenisation the context of damage. A POD-based reduced order modelling method was applied, combined with a system approximation technique. One big challenge in trying to apply this method is to select an appropriate snapshot space. The snapshot space was selected using various monotonic loadings paths on the RVE. The example shows a good speedup and a reasonable accuracy. The controlled elements (which are the key elements that monitor information used in the reduced order model) were selected using the discrete empirical interpolation method (DEIM). Those elements gather in the area of high damage and of high variation of displacement, intuitively being the most informative areas. This approach can be though as a bridge between analytical and computational homogenisation: the reduced bases are pseudo-analytical solution of the RVE problem that is still computationally solved at reduced cost. An interesting development to the method would be to have adaptive reduced bases changing depending on the loading applied and the state of the material.

Figure 1 :

 1 Figure 1: Lattice model of the RVE. Bars have different properties depending if they represent the matrix or the inclusions

RVEFigure 2 :

 2 Figure 2: Schematic representation of semi-concurrent computational homogenisation. The constitutive law of the macro-structure is unknown. An effective strain is applied as boundary condition to the RVE boundary value problem. Once solved, the local material properties are obtained.

Figure 4 :

 4 Figure 4: Representation of a RVE boundary value problem solution in the reduced order model setting: it is a sum of the coarse scale contribution and a fluctuation that is represented as a linear combination of basis vectors

  (a) Original structure (b) Example of a surrogate structure

Figure 5 :

 5 Figure 5: Example of a surrogate structure. The stiffness of the structure is evaluated on controlled elements only, while the other ones are just like ghosts

Figure 7 :

 7 Figure 7: Sample loading paths generated by the algorithm

Figure 8 :

 8 Figure 8: First 3 modes obtained through the POD. The damage localises between pairs of inclusions. 14

1

  Number of displacement basis vectors Number of static basis vectorsRelative error in energy norm

Figure 12 :

 12 Figure 12: Evolution of the error varying the number of displacement and static basis vectors.

Figure 13 :

 13 Figure 13: Evolution of the error with respect to the speedup while increasing the number of basis vectors.

In the sense of the k.k norm

Acknowledgements

The authors thank the financial support of EPSRC High End Computing Studentship for Mr. Olivier Goury as well as the support of Cardiff and Glasgow Universitys Schools of Engineering. Pierre Kerfriden and Stéphane Bordas would like to also acknowledge the financial support of the Royal Academy of Engineering and of the Leverhulme Trust for S. Bordas' Senior Research Fellowship entitled "Towards the next generation surgical simulator" which funded Pierre Kerfriden's post in 2009-2010 and allowed him to lay the foundations for the work presented in this paper.