Lucas Boczkowski

Emanuele Natale

Ofer Feinerman

Amos Korman

(Supplementary Information) Limits on Reliable Information Flows through Stochastic Populations

The models

We consider a population of n agents that interact stochastically and aim to converge on a particular opinion held by few knowledgable individuals. For simplicity, we assume that the set of opinions contain two opinions only, namely, 0 and 1.

As detailed in this section, we shall assume that agents have access to significant amount of resources, often exceeding more reasonable and realistic assumptions. Since we are concerned with lower bounds, we do not lose generality from such permissive assumptions. These liberal assumptions will actually simplify our proofs. One of these assumptions is the assumption that each agent is equipped with a unique identity id(v) in the range {1, 2, . . . , n} (see more details in Section Liberal assumptions).

Initial configuration

The initial configuration is described in several layers. First, the neutral initial configuration corresponds to the initial states of the agents, before the sources and the desired opinion to converge to are set. (The term neutral is motivated by a physical analogy, as opposed to a charged initial configuration.)

A random initialization is then applied to the given neutral initial configuration, which determines the set of sources and the opinion that agents need to converge to. This will result in what we call the charged initial configuration. It can represent, for example, an external event that was identified by few agents which now need to deliver their knowledge to the rest of the population.

Neutral initial configuration x (0)

Each agent v starts the execution with an input that contains, in addition to its identity:

• an initial state taken from some discrete set of states, and

• a binary opinion variable λ v ∈ {0, 1}.

(The opinion of an agent could have been considered as part of the state of the agent. We separate these two notions merely for the presentation purposes.) The neutral initial configuration x (0) is the vector whose i'th index, x (0) i for i ∈ {1, 2, . . . , n}, is the input of the agent with identity i.

Charged initial configuration and correct opinion

The charged initial configuration is determined in three stages. The first corresponds to the random selection of sources, the second to the selection of the correct opinion, and the third to a possible update of states of sources, as a result of being selected as sources with a particular opinion.

1st stage -Random selection of sources. Given an integer s ≤ n, a set S of size s is chosen uniformly at random (u.a.r) among the agents. The agents in S are called sources. Note that any agent has equal probability of being a source. We assume that each source knows it is a source, and conversely, each non-source knows it is not a source. 2nd stage -Random selection of correct opinion. In the main model we consider, after sources have been determined in the first stage, the sources are randomly initialized with an opinion, called the correct opinion. That is, a fair coin is flipped to determine an opinion in {0, 1} and all sources are assigned with this opinion. 3rd stage -Update of initial states of sources. To capture a change in behavior as a result of being selected as a source with a particular opinion, we assume that once the opinion of a source u has been determined, the initial state of u may change according to some distribution f source-state that depends on (1) its identity, (2) its opinion, and (3) the neutral configuration. Each source samples its new state independently.

Alphabet and messages

Agents communicate by observing each other according to some random pattern (for details see the Section Random interaction patterns). To improve communication agents may choose which content, called message, they wish to reveal to other agents that observe them. Importantly, however, such messages are subject to noise.

More specifically, at any given time, each agent v (including sources) displays a message m ∈ Σ, where Σ is some finite alphabet. The alphabet Σ agents use to communicate may be richer than the actual information content they seek to disseminate, namely, their opinions. This, for instance, gives them the possibility to express several levels of certainty [START_REF] Korman | Confidence sharing: An economic strategy for efficient information flows in animal groups[END_REF]. We can safely assume that the size of Σ is at least two, and that Σ includes both symbols 0 and 1. We are mostly concerned with the case where Σ is of constant size (i.e., independent of the number of agents), but note that our results hold for any size of the alphabet Σ, as long as the noise criterion is satisfied (see below).

δ-uniform noise

When an agent u observes some agent v, it receives a sample of the message currently held by v. The noise in the sample is characterized by a noise parameter 0 < δ ≤ 1/2. One of the important aspects in our theorems is that they are general enough to hold assuming any distribution governing the noise, as long as it satisfies the following noise criterion.

Definition 1 (The noise criterion with parameter δ). Any time some agent u observes an agent v holding some message m ∈ Σ, the probability that u actually receives a message m is at least δ, for any m ∈ Σ. We assume that all noisy samples are independent.

Observe that the aforementioned criterion implies that δ ≤ 1/|Σ|, and that the case δ = 1/|Σ| corresponds to messages being completely random, and the rumor spreading problem is thus unsolvable.

We next define a weaker criterion, that is particularly meaningful in cases in which sources are more restricted in their message repertoire than general agents. This may be the case, for example, if sources always choose to display their opinion as their message (possibly together with some extra symbol indicating that they are sources). Formally, we define Σ ⊆ Σ as the set of possible messages that a source can hold together with the set of messages that can be observed when viewing a source (i.e., after noise is applied). Our theorems actually apply to the following criterion, that requires that only messages in Σ are attained due to noise with some sufficient probability.

Definition 2 (Σ -relaxed noise ellipticity parameter δ). We say that the noise has Σ -relaxed ellipticity δ if P m,m ≥ δ for any m ∈ Σ and m ∈ Σ .

Random interaction patterns

We consider several basic interaction patterns. Our main model is the parallel-PULL model. In this model, time is divided into rounds, where at each round i ∈ N + , each agent u independently selects an agent v (possibly u = v) u.a.r from the population and then u observes the message held by v. The parallel-PU LL model should be contrasted with the parallel-PU SH model, in which u can choose between sending a message to the selected node v or doing nothing. We shall also consider the following variants of PU LL model.

• parallel-PULL(k). Generalizing parallel-PU LL for an integer 1 ≤ k ≤ n, the parallel-PU LL(k) model allows agents to observe k other agents in each round. That is, at each round i ∈ N + , each agent independently selects a set of k agents (possibly including itself) u.a.r from the population and observes each of them.

• sequential-PU LL. In each time step t ∈ N + , two agents u and v are selected uniformly at random (u.a.r) among the population, and agent u observes v.

• broadcast-PU LL. At each time step t ∈ N + one agent is chosen u.a.r. from the population and all agents observe it, receiving the same noisy sample of its message.

The broadcast-PU LL model is mainly used for technical considerations. We use it in our proofs as it simplifies our arguments while not harming their generality. Nevertheless, this broadcast model can also capture some situations in which agents can be seen simultaneously by many other agents, where the fact that all agents observe the same sample can be viewed as noise being originated by the observed agent.

Regarding the difference in time units between the models, since interactions occur in parallel in the parallel-PU LL model, one round in that model should informally be thought of as roughly n time steps in the sequential-PU LL or broadcast-PU LL model.

Liberal assumptions

As mentioned, we shall assume that agents have abilities that surpass their realistic ones. This does not only increases the generality of our lower bounds, but also simplifies their proofs. Specifically, the following liberal assumptions are considered.

• Unique identities. Each agent is equipped with a unique identity id(v) ∈ {1, 2, . . . , n}, that is, for every two agents u and v, we have id(u) = id(v). Moreover, whenever an agent u observes some agent v, we assume that u can infer the identity of v. In other words, we provide agents with the ability to reliably distinguish between different agents at no cost.

• Unlimited internal computational power. We allow agents to have unlimited computational abilities including infinite memory capacity. Therefore, agents can potentially perform arbitrarily complex computations based on their knowledge (and their id).

• Complete knowledge of the system. Informally, we assume that agents have access to the complete description of the system except for who are the sources and what is their opinion. More formally, we assume that each agent has access to:

the neutral initial configuration x (0) , all the systems parameters, including the number of agents n, the noise parameter δ, the number of sources s, and the distribution f source-state governing the update the states of sources in the third stage of the charged initial configuration.

• Full synchronization. We assume that all agents are equipped with clocks that can count time steps (in sequential-PU LL or broadcast-PU LL) or rounds (in parallel-PULL(k)). The clocks are synchronized, ticking at the same pace, and initialized to 0 at the beginning of the execution. This means, in particular, that if they wish, the agents can actually share a notion of time that is incremented at each time step.

• Shared randomness. We assume that algorithms can be randomized. That is, to determine the next action, agents can internally toss coins and base their decision on the outcome of these coin tosses. Being liberal, we shall assume that randomness is shared in the following sense. At the outset, an arbitrarily long sequence r of random bits is generated and the very same sequence r is written in each agent's memory before the protocol execution starts. Each agent can then deterministically choose (depending on its state) which random bits in r to use as the outcome of its own random bits. In particular, since agents are allowed to have distinct initial states (e.g. by having unique identity labels), they can choose to make use of disjoint sets of random bits, thus making use of independent random variables. On the other hand, the shared randomness also implies that, for example, two agents can possibly make use of the very same random bits or merely observe the outcome of the random bits used by the other agents. Furthermore, the above implies that, conditioning on an agent u being a non-source agent, all the random bits used by u during the execution are accessible to all other agents.

• Coordinated sources. Even though non-source agents do not know who the sources are, we assume that sources do know who are the other sources. This means, in particular, that the sources can coordinate their actions.

Algorithm

Upon observation, each agent can alter its internal state (and in particular, its message to be seen by others) as well as its opinion. In reality, the updates of these variables may follow different constraints.

In the case of ants for example, it may take a long time to change their message even if their internal state changes. As part our liberal approach, we allow agents to change any part of their internal state instantaneously.

The strategy in which agents update these variables is called "algorithm". As mentioned, algorithms can be randomized, that is, to determine the next action, agents can use the outcome of coin tosses in the sequence r (see the shared randomness assumption in Liberal assumptions). Overall, the action of an agent u at time t depends on:

1. the initial state of u in the charged initial configuration (including, in particular, the identity of u and whether or not it is a source), 2. the initial knowledge of u (including the system's parameters and the neutral configuration), 3. the time step t, and the list of its observations (history) up to time t -1, denoted x

(<t) u , 4
. the sequence of random bits r.

Convergence and time complexity

At any time, the opinion of an agent can be viewed as a binary guess function that is used to express its most knowledgeable guess of the correct opinion. The agents aim to minimize the probability that they fail to guess this opinion. In this context, it can be shown that the optimal guessing function is deterministic (see A remark about random guess functions, in the Appendix).

Definition 3. We say that convergence has been achieved if one can specify a particular non-source agent v, for which it is guaranteed that its opinion is the correct opinion with probability at least 2/3. The time complexity is the number of time steps (respectively, rounds) required to achieve convergence.

We remark that the latter definition encompasses all three models considered.

Remark 1 (Different sampling rates of sources). We consider sources as agents in the population but remark that they can also be thought of as representing the environment. In this case, one may consider a different rate for sampling a source (environment) vs. sampling a typical agent. For example, the probability to observe any given source (or environment) may be x times more than the probability to observe any given non-source agent. This scenario can also be captured by a slight adaptation of our analysis. When x is an integer, we can alternatively obtain such a generalization by considering additional artificial sources in the system. Specifically, we replace each source u i with a set of sources U i consisting of x sources that coordinate their actions and behave identically (recall that we assume that sources know who are the other sources and can coordinate their actions), simulating the original behavior of u i . Since the number of sources increases by a multiplicative factor of x, our lower bounds (see Theorem 4 and Corollary 14.1) decrease by a multiplicative factor of x 2 .

Related works in computer science

In Rumor Spreading problems (also referred to as Broadcast) a piece of information typically held by a single designated agent is to be disseminated to the rest of the population. It is the subject of a vast literature in theoretical computer science, and more specifically in the distributed computing community, see, e.g., [START_REF] Boczkowski | Minimizing message size in stochastic communication patterns: Fast self-stabilizing protocols with 3 bits[END_REF][START_REF] Censor-Hillel | Global computation in a poorly connected world: fast rumor spreading with no dependence on conductance[END_REF][START_REF] Demers | Epidemic algorithms for replicated database maintenance[END_REF][START_REF] Doerr | Asymptotically optimal randomized rumor spreading[END_REF][START_REF] Doerr | Stabilizing consensus with the power of two choices[END_REF][START_REF] Feinerman | Breathe before speaking: efficient information dissemination despite noisy, limited and anonymous communication[END_REF][START_REF] Goyal | Lower bounds for the noisy broadcast problem[END_REF][START_REF] Karp | Randomized rumor spreading[END_REF][START_REF] Pittel | On spreading a rumor[END_REF]. While some works assume a fixed topology, the canonical setting does not assume a network. Instead agents communicate through uniform PU SH/PULL based interactions (including the phone call model), in which agents interact in pairs with other agents independently chosen at each time step uniformly at random from all agents in the population. The success of such protocols is largely due to their inherent simplicity and fault-tolerant resilience [START_REF] Elsässer | On the runtime and robustness of randomized broadcasting[END_REF][START_REF] Karp | Randomized rumor spreading[END_REF]. In particular, it has been shown that under the PU SH model, there exists an efficient rumor spreading protocol that uses a single bit per message and can overcome flips in messages (noise) [START_REF] Feinerman | Breathe before speaking: efficient information dissemination despite noisy, limited and anonymous communication[END_REF].

The line of research initiated by El-Gamal [START_REF] El-Gamal | [END_REF], also studies a broadcast problem with noisy interactions. The regime however is rather different from ours: all n agents hold a bit they wish to transmit to a single receiver. This line of research culminated in the Ω(n log log n) lower bound on the number of messages shown in [START_REF] Goyal | Lower bounds for the noisy broadcast problem[END_REF], matching the upper bound shown many years earlier in [START_REF] Gallager | Finding parity in a simple broadcast network[END_REF].

Several works have investigated algorithmic properties of networks with unstable topological structure, such as ephemeral networks, evolving graphs and edge-Markovian evolving graphs [START_REF] Avin | How to Explore a Fast-Changing World (Cover Time of a Simple Random Walk on Evolving Graphs)[END_REF][START_REF] Clementi | Information Spreading in Stationary Markovian Evolving Graphs[END_REF][START_REF] Akrida | Ephemeral networks with random availability of links: The case of fast networks[END_REF]. Such works prove analytical results assuming that the evolution of the topology satisfies certain constraints, and did not consider the case of noise affecting communication in conjunction with the dynamicity of the topology.

The lower bounds

Throughout this section we consider δ < 1/|Σ|, such that (1-δ|Σ|) δsn ≤ 1 10 . Our goal in this section is to prove the following result.

Theorem 4. Assume that the relaxed δ-uniform noise criterion is satisfied.

• Let k be an integer. Any rumor spreading protocol on the parallel-PU LL(k) model cannot converge in fewer rounds than

Ω nδ ks 2 (1 -δ|Σ|) 2 .
• Consider either the sequential-PU LL or the broadcast-PU LL model. Any rumor spreading protocol cannot converges in fewer time steps than

Ω n 2 δ s 2 (1 -δ|Σ|) 2 .
To prove the theorem, we first prove (in Reducing to the broadcast-PU LL Model) that an efficient rumor spreading algorithm in either the noisy sequential-PULL model or the parallel-PULL(k) model can be used to construct an efficient algorithm in the broadcast-PU LL model. The resulting algorithm has the same time complexity as the original one in the context of sequential-PULL and adds a multiplicative factor of kn in the context of parallel-PU LL(k).

We then show how to relate the rumor spreading problem in broadcast-PULL to a statistical inference test (Rumor Spreading and hypothesis testing). A lower bound on the latter setting is then achieved by adapting techniques from mathematical statistics (Proof of Theorem 7).

Remark 2. The lower bound of Theorem 4 loses relevance when s is of order greater than √ n. Indeed, the following simple protocol in the broadcast-PU LL model turns out to match the lower bound. For simplicity's sake, let us consider the case of a binary alphabet Σ = {0, 1}, and assume without loss of generality that the sources' opinion is 1. Each non-source agent, at each time step, chooses a random message u.a.r. in Σ, while each source agent always displays the correct message. After n time steps, the agents have collected n observations. If s √ 10n, a straightforward application of the Chernoff bound shows that with high probability at least s/2 of the n observations come from source agents. Thus, at most n -s/2 of the observations have distribution Bernoulli( 12 ), while at least s/2 of them are identically 1 before the effect of noise is taken into account. Since s/2 is of the same order of the standard deviation of the non-source messages, the agents have a good probability to correctly infer the correct opinion by choosing the most frequent message among the n observations.

Reducing to the broadcast-PULL Model

The following lemma establishes a formal relation between the convergence times of the models we consider. We assume all models are subject to the same noise distribution. Lemma 5. Any protocol operating in sequential-PU LL can be simulated by a protocol operating in broadcast-PU LL with the same time complexity. Moreover, for any integer 1 ≤ k ≤ n, any protocol P operating in parallel-PULL(k) can be simulated by a protocol operating in broadcast-PU LL with a time complexity that is kn times that of P in parallel-PULL(k).

Proof. Let us first show how to simulate a time step of sequential-PU LL in the broadcast-PULL model. Recall that in broadcast-PULL, in each time step, all agents receive the same observation sampled u.a.r from the population. Upon drawing such an observation, all agents use their shared randomness to generate a (shared) uniform random integer X between 1 and n. Then, the agent whose unique identity corresponds to X is the one processing the observation, while all other agents ignore it. This reduces the situation to a scenario in sequential-PULL, and the agents can safely execute the original algorithm designed for that model.

As for simulating a time step of parallel-PU LL(k) in the broadcast-PULL model, agents divide time steps in the latter model into rounds, each composing of precisely kn time steps. Recall that the model assumes that agents share clocks that start when the execution starts and tick at each time step. This implies that the agents can agree on the division of time into rounds, and can further agree on the round number. For an integer i, where 1 ≤ i ≤ kn, during the i-th step of each round, only the agent whose identity is (i mod n)+1 receives the observation, while all other agents ignore it. Observe that receiving the observation doesn't imply that the agent processes this observation. In fact, it will store it in its memory until the round is completed, and process it only then. The aforementioned rule for receiving a message, ensures that when a round is completed in the broadcast-PU LL model, each agent receives precisely k independent uniform samples as it would in a round of parallel-PULL(k). Therefore, at the end of each round j ∈ N + in the broadcast-PULL model, all agents can safely execute their actions in the j'th round of the original protocol designed for parallel-PU LL(k). This draws a precise bijection from rounds in parallel-PU LL(k) and rounds in broadcast-PU LL. The multiplicative overhead of kn simply follows from the fact that each round in broadcast-PU LL consists of kn time steps.

Thanks to Lemma 5, Theorem 4 directly follows from the next theorem. Theorem 6. Consider the broadcast-PU LL model and assume that the relaxed δ-uniform noise criterion is satisfied. Any rumor spreading protocol cannot converge in fewer time steps than

Ω n 2 δ s 2 (1 -δ|Σ|) 2 .
The remaining of the section is dedicated to proving Theorem 6. Towards achieving this, we view the task of guessing the correct opinion in the broadcast-PU LL model, given access to noisy samples, within the more general framework of distinguishing between two types of stochastic processes which obey some specific assumptions.

Rumor Spreading and hypothesis testing

Consider the following class of problems.

Adaptive Coin Distinguishing Task (ACDT).

A distinguisher is presented with a sequence of observations taken from a coin of type η where η ∈ {0, 1}. We can think of the type η as initially set to 0 or 1 with probability 1/2 (independently of everything else). The goal of the distinguisher is to determine the type η, based on the observations. More specifically, for a given time step t, denote the sequence of previous observations (up to, and including, time t -1) by

x (<t) = (x (1) , . . . , x (t-1) ).

For each time t, we denote the probability distribution of the observation

X (t)
η ∈ Σ that the distinguisher receives, given the type η ∈ {0, 1} and the history of previous observations x (<t) , as

P (X (t) η = m | x (<t) ). (1)
We follow the common practice to use uppercase letters to denote random variables and lowercase letter to denote a particular realisation, e. g. X (≤t) for the sequence of observations up to time t, and x (≤t) for a corresponding realization.

Remark 3. We note that ACDT generalizes in several respects the canonical problem of distinguishing between an unbiased coin and a coin with fixed bias ε (see, e.g., Chapter 5 in [START_REF] Rigollet | High dimensional statistics[END_REF]). It is more general because 1. the probabilities of observations may vary adaptively as a function of the outcome of the previous samples, since the coins and p

(t) 1 (m) P (X (t) 0 = m | x (<t)
) and P (X

(t) 1 = m | x (<t)
) in (1) actually depend on x (<t) , the history of observations up to time t -1, and 2. instead of binary random variables we consider Σ-valued random variables.

We next introduce, for each m ∈ Σ, the parameter ε(m, x (<t) ) = P (X

(t) 1 = m | x (<t) ) -P (X (t) 0 = m | x (<t) ).
Since, at all times t, it holds that m∈Σ P (X

(t) 0 = m | x (<t) ) = m∈Σ P (X (t) 1 = m | x (<t) ) = 1, then m∈Σ ε(m, x (<t) ) = 0.
We shall be interested in the quantity

d ε (x (<t) ) := m∈Σ |ε(m, x (<t) )|,
which corresponds to the 1 distance between the distributions P (X

(t) 0 = m | x (<t) ) and P (X (t) 1 = m | x (<t)
) given the sequence of previous observations. The bounded family ACDT(ε, δ). We consider a family of instances of ACDT, called ACDT(ε, δ), governed by parameters ε and δ. Specifically, this family contains all instances of ACDT such that for every t, and every history x (<t) , we have:

• d ε (x (<t) ) ≤ ε, and
• for every m ∈ Σ such that ε(m, x (<t) ) = 0, we have δ ≤ P (X

(t) η = m | x (<t) ) for η ∈ {0, 1}.
In the rest of the current section, we show how Theorem 6, that deals with the broadcast-PULL model, follows directly from the next theorem that concerns the adaptive coin distinguishing task, by setting

ε = 2s(1 -δ|Σ|) n .
The actual proof of Theorem 7 appears in Proof of Theorem 7.

Theorem 7. Consider any protocol for any instance of ACDT(ε, δ), The number of samples required to distinguish between a process of type 0 and a process of type 1 with probability of error less than 1 3 is at least ln 2 9

6(δ -ε) 3 δ 3 -δ 2 ε + 3δε 2 -ε 3 δ ε 2 .
In particular, if 10ε < δ, then the number of necessary samples is Ω δ ε 2 .

Proof of Theorem 6 assuming Theorem 7

Consider a rumor spreading protocol P in the broadcast-PULL model. Fix a node u. We first show that running P by all agents, the perspective of node u corresponds to a specific instance of ACDT( 2s(1-δ|Σ|) n , δ) called Π(P, u). We break down the proof of such correspondence into two claims.

The ACDT instance Π(P, u). Recall that we assume that each agent knows the complete neutral initial configuration, the number of sources s, and the shared of random bits sequence r. We avoid writing such parameters as explicit arguments to Π(P, u) in order to simplify notation, however, we stress that what follows assumes that these parameters are fixed. The bounds we show hold for any fixed value of r and hence also when r is randomized.

Each agent is interested in discriminating between two families of charged initial configurations: Those in which the correct opinion is 0 and those in which it is 1 (each of these possibilities occurs with probability 1 2 ). Recall that the correct opinion is determined in the 2nd stage of the charged initial configuration, and is independent form the choice of sources (1st stage).

We next consider the perspective of a generic non-source agent u, and define the instance Π(P, u) as follows. Given the history x (<t) , we set P (X

(t) η = m | x (<t)
), for η ∈ {0, 1}, to be equal to the probability that u observes message m ∈ Σ at time step t of the execution P. For clarity's sake, we remark that the latter probability is conditional on:

• the history of observations being x (<t) ,

• the sequence of random bits r,

• the correct opinion being η ∈ {0, 1},

• the neutral initial configuration,

• the identity of u,

• the algorithm P, and

• the system's parameters (including the distribution f source-state and the number of sources s).

Claim 8. Let P be a correct protocol for the rumor spreading problem in broadcast-PU LL and let u be an agent for which the protocol is guaranteed to produce the correct opinion with probability at least p by some time T (if one exists), for any fixed constant p ∈ (0, 1). Then Π(P, u) can be solved in time T with correctness being guaranteed with probability at least p.

Proof. Conditioning on η ∈ {0, 1} and on the random seed r, the distribution of observations in the Π(P, u) instance follows precisely the distribution of observations as perceived from the perspective of u in broadcast-PU LL. Hence, if the protocol P at u terminates with output j ∈ {0, 1} at round T , after the T -th observation in Π(P, u) we can set Π(P, u)'s output to j as well. Given that the two stochastic processes have the same law, the correctness guarantees are the same.

Claim 9. Π(P, u) ∈ ACDT 2(1-δ|Σ|)s n , δ .
Proof. Since the noise in broadcast-PU LL flips each message m ∈ Σ into any m ∈ Σ with probability at least δ, regardless of the previous history and of η ∈ {0, 1}, at all times t we have

m ∈ Σ =⇒ P (X (t) η = m | x (<t) ) ≥ δ.
Consider a message m ∈ Σ \ Σ (if such a message exists). By definition, such a message could only be received by observing a non-source agent. But given the same history x (<t) , the same sequence of random bits r, and the same initial knowledge, the behavior of a non-source agent is the same, no matter what is the correct opinion η. Hence, for m ∈ Σ \ Σ we have P (X

(t) 0 = m | x (<t) ) = P (X (t) 1 = m | x (<t) ), or in other words, m ∈ Σ \ Σ =⇒ ε(m, x (<t) ) = 0.
It remains to show that d ε (x (<t) ) ≤ 2(1-δ|Σ|)s n . Let us consider two executions of the rumor spreading protocol, with the same neutral initial configuration, same shared sequence of random bits r, same set of sources, except that in the first the correct opinion is 0 while in the other it is 1. Let us condition on the history of observations x (<t) being the same in both processes.

As mentioned, given the same history x (<t) , the behavior of a non-source agent is the same, regardless of the correct opinion η. It follows that the difference in the probability of observing any given message is only due to the event that a source is observed. Recall that the number of sources is s. Therefore, the probability of observing a source is s/n, and we may write as a first approximation ε(m, x (<t) ) ≤ s/n . However, we can be more precise. In fact, ε(m, x (<t) ) is slightly smaller than s/n, because the noise can still affect the message of a source.

We may interpret ε(m, x (<t) ) as the following difference. For a source v ∈ S, let m v η be the message of u assuming the given history x (<t) and that v is of type η ∈ {0, 1} (the message m v η is deterministically determined given the sequence r of random bits, the neutral initial configuration, the parameters of the system, and the identity of v). Let α m ,m be the probability that the noise transforms a message m into a message m. Then

ε(m, x (<t) ) = 1 n v∈S (α m v 1 ,m -α m v 0 ,m ),
and d ε (x (<t) ) = m∈Σ |ε(m, x (<t) )| ≤ 1 n m∈Σ v∈S |α m v 1 ,m -α m v 0 ,m |. (2)
By the definition of ACDT(ε, δ), it follows that either

α m v 1 ,m = α m v 0 ,m (if ε(m, x (<t) ) = 0) or δ ≤ α m v 1 ,m , α m v 0 ,m ≤ 1 -δ (if ε(m, x (<t
) ) = 0). Thus, to bound the right hand side in (2), we can use the following claim (proven in Proof of Claim 10) Claim 10. Let P and Q be two distributions over a universe Σ such that for any element m ∈ Σ,

δ ≤ P (m), Q(m) ≤ 1 -δ. Then m∈Σ |P (m) -Q(m)| ≤ 2(1 -δ|Σ|).
Applying Claim 10 for a fixed v ∈ S to distributions (α m v 0 ,m ) m and (α m v 1 ,m ) m , we obtain

1 n m∈Σ v∈S |α m v 1 ,m -α m v 0 ,m | ≤ 1 n 2 v∈S (1 -δ|Σ|) ≤ 2(1 -δ|Σ|)s n .
Hence, we have Π(P) ∈ ACDT 2(1-δ|Σ|)s n , δ , establishing Claim 9.

Thanks to Claims 8 and 9, Theorem 6 regarding the broadcast-PULL model becomes a direct consequence of Theorem 7 on the adaptive coin distinguishing task, taking

ε = 2(1 -δ|Σ|)s n .
More precisely, the assumption (1-δ|Σ|) δsn ≤ c for some small constant c, ensures that ε δ ≤ c as required by Theorem 7. The lower bound Ω ε 2 δ corresponds to

Ω n 2 δ (1 -δ|Σ|) 2 s 2 .
This concludes the proof of Theorem 6. To establish our results it remains to prove Theorem 7.

Proof of Theorem 7

We start by recalling some facts from Hypothesis Testing.

We use the notation log(•) to denote the base 2 logarithms, i.e., log 2 (•) and for a probability distribution P , use the notation P (x) as a short for P (X = x). First let us recall two standard notions of (pseudo) distances between probability distributions. Given two discrete distributions P 0 , P 1 over a probability space Ω with the same support, the total variation distance is defined as T V (P 0 , P 1 ) := 1 2

x∈Ω |P 0 (x) -P 1 (x)|,
and the Kullback-Leibler divergence KL(P 0 , P 1 ) is defined as KL(P 0 , P 1 ) := x∈Ω P 0 (x) log P 1 (x) P 0 (x) .

The assumption that the support is the same is not necessary but it is sufficient for our purposes, and is thus made for simplicity's sake.

The following lemma shows that, when trying to discriminate between distributions P 0 , P 1 , the total variation relates to the smallest error probability we can hope for.

Lemma 11 (Neyman-Pearson [17, Lemma 5.3 and Proposition 5.4]). Let P 0 , P 1 be two distributions. Let X ∈ Ω be a random variable of distribution either P 0 or P 1 . Consider a (possibly probabilistic) mapping f : Ω → {0, 1} that attempts to "guess" whether the observation X was drawn from P 0 (in which case it outputs 0) or from P 1 (in which case it outputs 1). Then, we have the following lower bound, P 0 (f (X) = 1) + P 1 (f (X) = 0) ≥ 1 -T V (P 0 , P 1 ).

The total variation is related to the KL divergence by the following inequality.

Lemma 12 (Pinsker [START_REF] Rigollet | High dimensional statistics[END_REF]Lemma 5.8]). For any two distributions P 0 , P 1 , T V (P 0 , P 1 ) ≤ KL (P 0 , P 1 ).

We are now ready to prove the theorem.

Recall that we assume ε(m, x (<t) ) P (X

(t) 0 = m | x (<t) ) ≤ ε(m, x (<t) ) δ ≤ ε δ .
We make use of the following claim, which follows from the Taylor expansion of log(1 + u) around 0. More details can be found in the Appendix, Proof of Claim 14.

Claim 14. Let x ∈ [-a, a] for some a ∈ (0, 1). Then |log(1 + x) -x + x 2 /2| ≤ x 3 3(1-a) 3 . Using Claim 14 with a = ε δ , we can bound the inner sum appearing in [START_REF] Doerr | Asymptotically optimal randomized rumor spreading[END_REF] from above and below with

1 ln 2 m∈Σ ε(m, x (<t) ) - 1 2 (ε(m, x (<t) )) 2 P (X (t) 0 = m | x (<t) ) ± δ 3 3(δ -ε) 3 (ε(m, x (<t) )) 3 P (X (t) 0 = m | x (<t) ) 2 . ( 8
)
Since m |ε(m, x (<t) )| ≤ ε, we also have that m (ε(m, x (<t) )) 2 ≤ ε 2 . The latter bound, together with the fact that P (X

(t) 0 = m | x (<t) ) ≥ δ for any m ∈ Σ such that ε( m, x (<t) ) = 0, implies m (ε(m, x (<t) )) 2 P (X (t) 0 = m | x (<t) ) ≤ ε 2 δ . (9)
Finally, we can similarly bound the term m∈Σ (ε(m, x (<t) )) 3 /P (X

(t) 0 = m | x (<t) ) 2 with m∈Σ (ε(m, x (<t) )) 3 /P (X (t) 0 = m | x (<t) ) 2 ≤ ε 3 δ 2 . ( 10
)
Recall that m ε(m, x (<t) ) = 0, thus the first term in (8) disappears. Hence, substituting the bounds ( 9) and ( 10) in ( 8), we have

P (X(t) = m|x(< t)) log 1 + ε(m, x (<t) ) P (X (t) 0 = m | x (<t) ) ≤ 1 ln 2 1 2 ε 2 δ + δε 3 3(δ -ε) 3 • m P (X(t) = m|x(< t)) ≤ 1 ln 2 1 2 + δ 2 ε 3(δ -ε) 3 ε 2 δ • 1. (11)
If we define the right hand side [START_REF] Feinerman | Breathe before speaking: efficient information dissemination despite noisy, limited and anonymous communication[END_REF] to be W (ε, δ) and we substitute the previous bound in (7), we get

KL(P 0 (x (t) | x (<t) ), P 1 (x (t) | x (<t) )) ≤ W (ε, δ),
and combining the previous bound with (5), we can finally conclude that for any integer T , we have

KL(P (≤T ) 0 , P (≤T ) 1 ) ≤ T • W (ε, δ).
Thus, from Theorem 13 and the latter bound, it follows that the error under a uniform prior of the source type, as defined in ( 4), is at least

1 2 P 0 f (X (≤t) ) = 1 + 1 2 P 1 f (X (≤t) ) = 0 ≥ 1 2 - 1 2 KL(P (≤T ) 0 , P (≤T ) 1 ) ≥ 1 2 - 1 2 T • W (ε, δ).
Hence, the number of samples T needs to be greater than

1 9 1 W (ε, δ) = ln 2 9 6(δ -ε) 3 δ 3 -δ 2 ε + 3δε 2 -ε 3 δ ε 2 .
to allow the possibility that the error be less than 1/3.

In particular, if we assume that 10ε < δ, then we can bound

δ 2 ε 3(δ -ε) 3 ≤ δ 3 10 • 1
W (ε, δ) ≥ 0.14 • δ ε 2 = Ω δ ε 2 .
This completes the proof of Theorem 7 and hence of Theorem 6.

Detectable sources

In this section, we aim to prove the following.

Corollary 14.1. Consider the setting in which sources are reliably detectable. Assume that the relaxed δ-uniform noise criterion is satisfied and let T =

n 2 δ s 2 (1-δ|Σ|) 2 1/3 .
• Consider the sequential-PU LL model. Assume that sT ≥ C log n, for a large enough constant C.

Any rumor spreading scheme cannot converge in less than Ω (T ) time steps.

• Let k be an integer. Assume that sT /k ≥ C log n, for a large enough constant C. Any rumor spreading protocol in the parallel-PU LL(k) model cannot converge in less than Ω(T /k) rounds.

Proof. Let us start with the first item of the corollary, namely the lower bound in the sequential-PULL model. For any step t, let S(t) denote the set of sources together with the agents that have directly observed at least one of the sources at some point up to time t. We have S = S(0) ⊆ S(1) ⊆ S(2) ⊆ . . .. The size of the set S(t) is a random variable which is expected to grow at a moderate speed. Specifically, letting s = This last expression is a sum of n•T independent Bernoulli variables with parameter s/n. In other terms, it is a binomial variable with probability s/n and T •n trials. By a standard Chernoff bound the probability that it deviates by a multiplicative factor 11 10 from its mean s • T is less than exp(-Ω(sT )) ≤ n -10 . The last bound holds because we assume sT ≥ C log n for some large enough constant C.

Denote by E the event that |S(t)| ≤ s for every t ≤ T . Using Claim 15, we know that P (E) ≥ 1-n -10 . Our next goal is to prove that the probability ρ that a given agent correctly guesses the correct opinion is low for any given time t ≤ cT , where c is a small constant. For this purpose, we condition on the highly likely event E. Removing this conditioning will amount to adding a negligible term (of order at most n -10 ) to ρ.

In order to bound ρ, we would like to invoke Theorem 6 with the number of sources upper bounded by s . Let us explain why it applies in this context. To begin with, we may adversarially assume (from the perspective of the lower bound) that all agents in S(t) learn the value of the correct bit to spread. Thus, they essentially become "sources" themselves. In this case the number of sources varies with time, but the proof of Theorem 6 can easily be shown to cover this case as long as s (i.e., s here) is an upper bound on the number of sources at all times. We can therefore safely apply Theorem 6 with s . By the choice of T ,

T = Θ n 2 δ (s ) 2 (1 -δ|Σ|) 2 .
Hence, we can set c to be a sufficiently small constant such that for all times t ≤ cT , the probability of guessing correctly, even in this adversarial scenario, is less than 1/3. In other words, we have ρ ≤ 1/3. All together, this yields a lower bound of Ω(T ) on the convergence time.

As for the parallel-PU LL(k) model, the argument is similar. After T = T /k parallel rounds, using a similar claim as Claim 15, we have that with high probability, at most O(ksT ) agents have directly observed one of the s sources by time T . Applying Theorem 6 with s = O(ksT ) = O(sT ) yields a lower bound (in terms of samples in the broadcast model) of

Θ n 2 δ (s ) 2 (1 -δ|Σ|) 2 = Θ n 2 δ s 2 T 2 (1 -δ|Σ|) 2 = Θ(T ).
The last line follows by choice of T . Hence T is a lower bound on the number of samples, which is attained in T rounds of parallel-PU LL(k) model.

  m

11 10 • 1 ( 1 (

 1011 s • T , we obtain: Claim 15. With probability at least 1 -n -10 , we have |S(T )| ≤ s . Proof of Claim 15. The variable S(T ) may be written as a sum of indicator variables S(T ) = n i=1 Agent i observed at least one source before step t) Agent i observes a source on step r).

Proof of Theorem 7. Let us define P η ( • ) = P ( • | "correct distribution is η") for η ∈ {0, 1}. We denote P (≤t) η , η ∈ {0, 1}, the two possible distributions of X (≤t) . We refer to P (≤t) 0 as the distribution of type 0 and to P (≤t) 1 as the distribution of type 1. Furthermore, we define the correct type of a sequence of observations X (≤t) to be 0 if the observations are sampled from P (≤t) 0 , and to be 1 if they are sampled from P (≤t) 1 . After t observations x (≤t) = (x (1) , . . . , x (t) ) we have to decide whether the distribution is of type 0 or 1. Our goal is to maximize the probability of guessing the type of the distribution, observing X (≤t) , which means that we want to minimize

Recall that the correct type is either 0 or 1 with probability 1 2 . Thus, the error probability described in (3) becomes

By combining Lemmas 11 and 12 with X = X (≤t) and P η = P (≤t) η for η = 0, 1, we get the following Theorem. Although for convenience we think of f as a deterministic function, it could in principle be randomized (see A remark about random guess functions, in the Appendix).

Theorem 13. Let f be any guess function. Then

Theorem 13 implies that for the probability of error to be small, it must be the case that the term KL P (≤t) 0 , P (≤t) 1 is large. Our next goal is therefore to show that in order to make this term large, t must be large.

Note that P (≤T ) η for η ∈ {0, 1} cannot be written as the mere product of the marginal distributions of the X (t) s, since the observations at different times may not necessarily be independent. Nevertheless, we can still express the term KL(P

) as a sum, using the Chain Rule for KL divergence. It yields

where

Since we are considering an instance of ACDT (ε, δ), we have

We make use of the previous facts to upper bound the KL divergence terms in the right hand side of (6), as follows (recall that we omit the dependency of p (<t) and ε (<t) on the past observations x (<T ) ∈ Σ t-1 , in the interest of readability).

Missing proofs

Proof of Claim 10

The proof of the Claim could be slightly shortened using the well-known fact that the l 1 distance between two probability distribution equals two times the total variation between them. We derive a proof from first principles for completeness.

In general, define Σ + , Σ -and Σ = respectively as {m :

Note that for any set S, Q(S) ≥ δ|S| using the fact that Q(m) ≥ δ for any m and similarly for P . Thus the last line is bounded by

This completes the proof of Claim 10.

A remark about random guess functions

In this section, we show that we may relax the guessing function f to be probabilistic. In other words, allowing P (f (x (≤t) ) = 1) = p for some x(≤t) and some probability 0 < p < 1, would not allow to reduce further the value of (4). Indeed, if we consider a probabilistic f (x (≤t) ) in ( 4), we could rewrite the latter as a convex combination of deterministic guessing functions. To illustrate this, consider for example the case in which f is random only on a particular input x(≤t) , and define the two deterministic guessing function

and

From the above definition (and the law of total probability), it follows that

The above calculation can be generalized to the case in which f is random on any subset of inputs (possibly all). Thus, our results still hold for probabilistic guess functions. Informally, this means that we can allow agents to take decisions by "flipping a coin". We assume f to be deterministic for the sole purpose of easing the presentation.

Proof of Claim 14

Using a Taylor expansion of log(1 + u) of order 3 around 0 and the Remainder Theorem, we obtain, for any x ∈ [-a, a],

|log(1 + x) -x + x 2 /2| ≤ x 3 6 max y∈[-a,a]

(log(1 + y)) (3) , where f (3) for a function f stands for the third derivative of f . Since (log(1 + y)) (3) = 2 (1 + y) 3 , Claim 14 follows.


Advanced embedded algorithms are growing in complexity and length, related to the growth in autonomy, which allows systems to plan paths of their own. However, this promise cannot happen without proper attention to the considerably stronger operational constraints that safety-critical applications must meet. This paper discusses the formal verification for optimization algorithms with a particular emphasis on receding-horizon controllers. Following a brief historical overview, a prototype autocoder for embedded convex optimization algorithms will be discussed. Options for encoding code properties and proofs, and their applicability and limitations will be detailed as well.

Formal Verification for Embedded Implementation of Convex Optimization Algorithms

Raphael Cohen * , * * Guillaume Davy * * , * * * Eric Feron * Pierre-Loïc Garoche * *

INTRODUCTION

The need for more safety and better performance is currently pushing the introduction of advanced numerical methods into next generations of cyber-physical systems.

While most of the algorithms described in this paper have been known for a long time, their online use within embedded systems is relatively new and opens issues that have to be addressed. Among these methods, we are concerned specifically with numerical optimization algorithms.

Problem: These algorithms solve a constrained optimization problem, defined by an objective function -the cost function -and a set of constraints to be satisfied:

min f 0 (x) s.t. f i (x) ≤ b i for i ∈ [1, m] (1)
This problem searches for x ∈ R n the optimization variable, minimizing f 0 ∈ R n → R, the objective function, while satisfying constraints f i ∈ R n → R, with associated bound b i . An element R n is feasible when it satisfies all the constraints f i . An optimal point is defined by the element having the smallest cost value among all feasible points. An optimization algorithm computes an exact or approximated estimate of the value of the optimal cost, together with one or more feasible points achieving this value. A subclass of these problems can be efficiently solved: convex problems. In these cases, the functions f 0 and f i are required to be convex and the optimal solution is unique. We refer the reader to [START_REF] Boyd | Convex optimization[END_REF]; [START_REF] Nocedal | Numerical Optimization[END_REF] for more details on convex optimization.

Context: Recently this formalism has been used with great success for the guidance of safety-critical application. Such applications include autonomous cars [START_REF] Jerez | Embedded online optimization for model predictive control at megahertz rates[END_REF] and reusable rockets [START_REF] Açikmese | Lossless convexification of nonconvex control bound and pointing constraints of the soft landing optimal control problem[END_REF]; [START_REF] Blackmore | Lossless convexification of control constraints for a class of nonlinear optimal control problems[END_REF]. The latter case has resulted in a spectacular experiments, including landings of SpaceX's Falcon 9 and BlueOrigin's New Shepard. Thus, powerful algorithms solving optimization problems do exist and are already used online. Although those algorithms have been embedded on board and running on an actual system, they still lack the level of qualification required by civil aircraft or manned rocket flight. computing a solution in a given time.

Algorithms: Let us consider the specific case of Linear Programming problems, a subclass of convex optimization. Multiple algorithm methods are available to solve them. First, the simplex method, which has been used successfully since the 1940s, despite its exponential worstcase complexity. Then the Ellipsoid algorithm and the very popular interior-point methods. These latter both exhibit polynomial complexity, though only interior point methods are considered to be practical for large desktop applications. These two methods are also applicable to more general settings such as quadratic programming (QP), affine constraints and quadratic cost, second order conic programming (SOCP), quadratic constraints and linear cost, or semidefinite programs (SDP).

When optimization algorithms are used offline, the soundness of their implementation and the feasibility of the computed optimizers is not as critical. Solutions could be validated a posteriori [START_REF] Roux | Formal proofs of rounding error bounds[END_REF]; [START_REF] Roux | Validating numerical semidefinite programming solvers for polynomial invariants[END_REF]. This paper is concerned with online use of such algorithms, providing evidence of the a priori validity of the computed solution. This paper focuses on the linear programming paradigm and makes the following contributions:

• high level properties defining a sound optimization algorithm are formalized; • these properties are expressed directly on the code artifact as annotations;

• the evidence supporting the algorithms properties is expressed at the code level • the approach is demonstrated on the Ellipsoid algorithm and on an instance of the Interior Point method.

We believe this paper addresses a major certification issue that can considerably influence the use of online optimization methods in safety-critical applications.

Related work. Several authors including [START_REF] Richter | Certification aspects of the fast gradient method for solving the dual of parametric convex programs[END_REF][START_REF] Mcgovern | Computational Analysis of Real-Time Convex Optimization for Control Systems[END_REF]; [START_REF] Mcgovern | Requirements and hard computational bounds for real-time optimization in safety-critical control systems[END_REF] have worked on the certification problem of optimization algorithms for their online uses in control, in particular on worst-case execution time issues. In those cases, the authors have chosen to tackle the problem at a high level of abstraction.

Formally verifying high level properties of numerical algorithm implementations was performed on the other context of linear controller stability [START_REF] Feron | From control systems to control software[END_REF]. This theoretical work was later instantiated concretely, generating code specifications and proving it with respect to the code Herencia-Zapana et al. ( 2012); Wang et al. (2016a). Some other work develop that approach for weaker properties [START_REF] Araiza-Illan | Verification of control systems implemented in simulink with assertion checks and theorem proving: A case study[END_REF]; [START_REF] Pajic | Automatic verification of linear controller software[END_REF].

Our work follows a similar approach with respect to [START_REF] Wang | Credible autocoding of convex optimization algorithms[END_REF] in which interior-point method algorithms are annotated with convergence proof elements in order to show the soundness of the imperative implementation. This work however, remains theoretical, expressed in Matlab level, without any proof performed on the actual code.

Structure. The paper is structured as follows: Section 2 presents backgrounds for linear programming and axiomatic semantics using Hoare triples. Sections 3 and 4 focus on the two considered approaches: Interior Point method and Ellipsoid method. Section 5 presents our early proof results on the code artifact, while Section 7 concludes.

PRELIMINARIES

In order to support the following analysis, this section introduces the notions and notations used throughout the paper. First, we will discuss Linear Programming (LP) problems. Then we introduce axiomatic semantics and Hoare logic.

Linear Programming

Linear Programming is a class of optimization problems.

A linear programming problem is defined by a matrix A of size m × n, and two vectors b, c of respective size m, n. Definition 1. Let us consider A ∈ R m×n , b ∈ R m and c ∈ R n . We will write P (A, b, c) as the linear program:

min Ax≥b c, x . (2)
Definition 2. Let us consider the problem P (A, b, c) and let us assume that an optimal point exists and is unique.

We have the following definitions:

/* @ logic LMat mat_add ( LMat x0, LMat x1 ) ; @ axiom getM_add : \forall LMat A, B ; getM ( mat_add ( A, B ) == getM ( A ) ; @ axiom getN_add : \forall LMat A, B ; getN ( mat_add ( A, B ) == getN ( A ) ; @ axiom get_add : @ \forall LMat A, B ; ( @ \forall integer i, j ; @ mat_get ( mat_add ( A, B ) , i, j ) == mat_get ( A, i, j ) @ + mat_get ( B, i, j ) ) ; */ Fig. 1. Axiomatization of matrix addition

• E f = {x | Ax ≥ b} the feasible set of P , • f (x) = c T x = c, x the cost function and • x * = arg min x∈E f
f the optimal point.

Axiomatic semantics and Hoare Logic

Semantics of programs express their behavior. For the same program, different means can be used to specify it: (i) a denotational semantics, expressing the program as a mathematical function, (ii) an operational semantics, expressing it as a sequence of basic computations, or (iii) an axiomatic semantics. In the latter case, the semantics can be defined in an incomplete way, as a set of projective statements, ie. observations. This idea was formalized by [START_REF] Floyd | Assigning meanings to programs[END_REF] and then [START_REF] Hoare | An axiomatic basis for computer programming[END_REF] as a way to specify the expected behavior of a program through pre-and post-condition, or assume-guarantee contracts.

Hoare Logic. A piece of code C is axiomatically described by a pair of formulas (P, Q) such that if P holds before executing C, then Q should be valid after its execution. This pair acts as a contract for the function and (P, C, Q) is called a Hoare triple. In most uses P and Q are expressed in first order formulas over the variables of the program. Depending on the level of precision of these annotations, the behavior can be fully or partially specified. In our case we are interested in specifying, at code level, algorithm specific properties such as the convergence of the analysis or preservation of feasibility for intermediate iterates.

Software frameworks, such as the Frama-C platform, cf. [START_REF] Cuoq | Frama-c: a software analysis perspective[END_REF], provide means to annotate a source code with these contracts, and tools to reason about these formal specifications. For the C language, ACSL, cf. [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF] (ANSI C Specification language) can be used as source comments to specify function contracts, or local annotations such as loop invariants.

Statement local annotations act as cuts in proofs and are typically required when analyzing loops.

Linear Algebra-based Specification. ACSL also provides means to enrich underlying logic by defining types, functions, and predicates. In its basic version, ACSL contains no predicate related to linear algebra. It is however possible to introduce such notions, supporting the expression of more advanced properties. Figure 1 presents the formalization of matrix addition.

These definitions are straightforward and fully defined. We also wrote a library for operators used in linear /* @ logic real norm ( LMat V, LMat X ) = @ sqrt ( mat_get ( mat_mul ( @ mat_mul ( transpose ( V ) , inv ( hess ( X ) ) ) @ , V ) , 0 , 0) ) ; */ Lemmas can also be defined to support later analyzers to prove annotations.To preserve proof consistency, the Frama-C tool requires to prove these additional lemmas. For example, splitting a complex theorem into several small lemmas usually helps analyzers in proving the whole goal.

INTERIOR-POINT METHOD

We focus here on Interior Point methods, that, starting from a feasible solution, solve iteratively a sequence of linear problems, eg. using Newton methods, to reach the optimal solution. This method presented by [START_REF] Karmarkar | A new polynomial-time algorithm for linear programming[END_REF] and improved by [START_REF] Nesterov | Introductory lectures on convex optimization : a basic course[END_REF]Nemirovski (1988, 1994) is one the most efficient methods to solve linear and convex optimization problems. The presented work directly follows the theory presented in (Nesterov, 2004, S 4 ). To keep the presentation simple and show the feasibility of our approach, we chose a primal method starting from the analytic center, whereas common methods are primal-dual and start from any point of the domain.

Let us introduce the major concepts before presenting the algorithm.

Basic concepts: Barrier functions and Central path

Barrier Function. Finding the optimum of a function is usually achieved by computing its derivative and searching for its zeros. In LP, we know that the optimum value is located on the border of the feasible set, which could be problematic to reach numerically. To solve this issue, we add a barrier function F to the cost function. The rationale of the barrier function F (x) is to diverge when x get closer to the border of the feasible set.

In order to guarantee convergence bounds, the barrier function must be a Lipschitz function. In practice, a self-concordant barrier function [START_REF] Nesterov | Self-Concordant functions and polynomial time methods in convex programming. Materialy po matematicheskomu obespecheniiu EVM[END_REF] satisfies this requirement. The typical choice is a sum of logarithmic functions. Definition 3. (Log barrier function). Let F be the barrier function associated to problem P . We have:

F (x) = m i=1 -ln( a i , x -b i ) (3)
where a i is the i-th row of A. Definition 4. (Analytic center). We define the analytic center as the point x * F verifying:

x * F = arg min x∈E f F (x) (4)
f (t, x) = tf (x) + F (x) (5)
In this function, t is the parameter that weights the relative importance of the initial cost function with respect to the barrier function. When t = 0 the original cost function vanishes and minimizing f leads to the analytic center while when t → +∞ the cost function is predominant and minimizing f leads to the optimal point up to (Figure 3) Definition 6. (Central path). The central path is defined as the curve:

x * (t) = arg min

x∈E f f (t, x) for t ∈ [0, +∞) (6)
Interior-point methods are part of the path-following optimization methods: starting from the analytic center and following the central path to reach the optimal point, as illustrated in Figure 3.

Approximate Centering Condition (ACC). Since the algorithm is doing numerical computation, the computed points cannot be located exactly on the central path. Thus, we need to introduce an approximate centering condition: Definition 7. (ACC). Let us define λ(t, x) as:

λ(t, x) = f (t, x) T F (x) f (t, x) (7)
The approximate centering condition becomes:

λ(t, x) < β 2 (8)
Notice that F (x) induces a norm1 so λ measures how close to 0 the derivative is.

This definition aims to guarantee that the points are close to the central path so that the last point can be close enough to the optimal: Theorem 1. (cf. [START_REF] Nesterov | Introductory lectures on convex optimization : a basic course[END_REF]

, S 4 )). If λ(t, x) < β 2 then c, x -c, x * ≤ ∆ t with ∆ = 1 + (β+1)β 1-β
/* @ requires acc (0 , MatVar ( X,2,1 ) ) ; @ ensures dot ( c, MatVar ( X,2,1 ) ) -dot ( c,sol ( A,b,c ) ) <EPSILON ; @ ensures A * MatVar ( X,2,1 ) > b ; */ void pathfollowing () { t = 0; /* @ loop invariant ensures acc ( t, MatVar ( X, 2 , 1) ) ; loop invariant A * MatVar ( X, 2 , 1) > b ; loop invariant ensures t > LOWER ( l ) ; */ for ( unsigned int l = 0; l<NBR ; l ++) { compute_pre () ; compute_dt () ; compute_dx () ; t = t + dt ; for ( unsigned int i = 0; i<N ; i ++) X 4. Shape of the path following algorithm While t increases, computed points converge to the optimal. The approximate centering condition is the invariant that we maintain along the iterations of the algorithm.

[ i ] = X [ i ] + dx [ i ]; } } Fig.

Algorithm

Since we are targeting embedded systems, we will only use static variables, a common practice in embedded systems. We also define macros representing the different values of the problems, simplifying the definition of annotations.

Figure 4 presents the main function of the algorithm, computing the next iterate while following the central path. It relies on:

• acc, dot and sol, ACSL self-defined operators representing, the approximate centering condition, the dot product and the solution of P (A, b, c) respectively. • several constants representing:

• EPSILON: the desired precision • M x N: the size of the matrix • LOWER, a function associating a lower bound for t at each iteration . • NBR: the number of iterations required so that lower(N BR) > ∆ • subfunctions compute_pre, compute_dt, compute_dx

In this program, we need to prove that at the end of the loop we ensure : dot(c, x) -dot(c, sol(A, b, c)) < EP SILON. This is proved thanks to a lemma we call acc_solution and the main function loop invariants. From the first loop invariant and acc_solution, one deduces that dot(c, x)dot(c, sol(A, b, c)) < ∆ t and from the second loop invariant, that t > ∆ , leading to the expected post-condition.

We will now detail the content and contracts of the three functions used by the main one: compute_dt, compute_dx and compute_pre.

Function compute_pre computes the gradient F (X) and the hessian F (X) of F at point X according to the Equs. 9 and 10.

/* @ requires acc ( t, MatVar ( X,2,1 ) ) ; @ ensures MatVar ( H,2,2 ) == hess ( MatVar ( X,2,1 ) ) ; @ ensures dt > t *0.125; @ ensures dt = (5/36) / norm ( c,MatVar ( X,2,1 ) ) ; @ */ void compute_dt () ;

Fig. 5. Contract on compute_dt (G) j (X) = m-1 i=0 - a i,j a i , X -b i (9) (H) j,k (X) = m-1 i=0 a i,j a i,k a i , X -b i (10)
Function compute_dt computes the update of dt according to Eq. 11. Following [START_REF] Nesterov | Introductory lectures on convex optimization : a basic course[END_REF], we have

dt = 5 36 • 1 c T × H × c (11)
We can deduce that dt > 0.125t. From this, we can lower bound t by a geometric progression: lower(n) = t 0 × (1.125) n . This allow us to determine the fixed number of iterations required to reach a given level of precision.

Function compute_dx updates X using a Newtonian iteration on the derivative of F :

dx = H -1 (tc + G) (12)
Notice that t in this formula is the one already updated (t + dt). We expect from this update that the approximate centering condition remains valid.

THE ELLIPSOID METHOD

The other method that has attracted our attention is the Ellipsoid Method. Despite its relative efficiency with respect to Interior Point methods, the Ellipsoid method benefits from concrete proof elements and could be considered viable option for critical embedded systems where safety is more important than performance. Before recalling the main steps of the algorithm, the needed elements will be presented.

Preliminaries

Ellipsoids in R n . An ellipsoid can be characterized by a positive definite matrix. A symmetric matrix P of R n×n is positive definite (S + n ) when

x T P x > 0 , ∀x = 0 ∈ R n
Furthermore, this definition implies that all matrices P in S + n are invertible. Definition 8. (Ellipsoid Sets). Let x ∈ R n and P ∈ R n×n positive-definite, we denote the Ellipsoid E(P, x) the set :

E(P, x) = {z ∈ R n : (z -x) T P -1 (z -x) ≤ 1} Definition 9. (Euclidean ball). Let n ∈ N we denote V n the unit Euclidean ball in R n . Vol(V n ) denotes its volume.
In addition, we define B(x, R) as the ball of radius R centered on x.

Algorithm

Let us now recall the main steps of the algorithm detailed in [START_REF] Bland | The ellipsoid method: A survey[END_REF]; [START_REF] Grötschel | The ellipsoid method and its consequences in combinatorial optimization[END_REF]; [START_REF] Khachiyan | Polynomial algorithms in linear programming[END_REF]; [START_REF] Boyd | Linear controller design: limits of performance[END_REF]. Remark 1. In the following, we denote E k , x k and P k , the ellipsoid, the corresponding center and positive-definite matrix, respectively computed by the algorithm at the k -th iteration.

e representing the vector normal to the chosen cutting hyperplane from oracle separation.

Fig. 6. Ellipsoid Cut for an iteration

Algorithm. We start the algorithm with an ellipsoid containing the feasible set E f , and therefore the optimal point x * . We iterate by transforming the current ellipsoid E k into a smaller volume ellipsoid E k+1 that also contains x * . Given an ellipsoid E k of center x k , we find a hyperplane containing x k that cuts E k in half, such that one half is known not to contain x * . Finding such a hyperplane is called the oracle separation step, cf. [START_REF] Ben-Tal | Lecture notes, optimizattion i-ii, convex analysis, non-linear programming theory, non-linear programming algorithms[END_REF]. In our LP setting, this step is not computationally expensive. Then, we define the ellipsoid E k+1 by the minimal volume ellipsoid containing the half ellipsoid Êk that is known to contain x * . In addition to that, we can compute an upper bound γ of the ratio of Vol(E k+1 ) to Vol(E k ). The Figures 6 and7 illustrates such ellipsoids cuts.

Ellipsoid transformation. From the oracle separation step, a separating hyperplane, e, that cuts E k in half with the guarantee that x * is localized in Êk has been computed. The following step is the Ellipsoid transformation.

Using this hyperplane e, one can update the Ellipsoid E k to its next iterate E k+1 according to equations 13 and 14.

x k+1 = x k - 1 n + 1 P k ẽ (13)
and

P k+1 = n 2 n 2 -1 (P k - 2 n + 1 P k ẽẽ T P k ) (14)
with: ẽ = e e T P k e . One can now characterize the upper bound of the ratio of Vol(E k+1 ) to Vol(E k ). Property 1. (Reduction ratio). Let k ≥ 0 be an iteration of the algorithm. We have

V ol(E k+1 ) ≤ exp( -1 2 • (n + 1) ) • V ol(E k )
Hypotheses. We recall that we assumed the matrix A and the vector b are such that the feasible set E f is bounded and not empty. Additionally, further information about E f is needed. Indeed, in order to know the number of steps required for the algorithm to return an -optimal solution, two scalars are needed:

• a radius R such that E f is included in B(0 n , R) • another scalar r such that there exist a center c 1 such that B(c 1 , r) is included in E f .
The main result can be stated as: Theorem 2. Let us assume that E f is bounded, not empty and such that R, r are known. Then, for all > 0, c ∈ R n there exists a N such that the algorithm, using N iterations, will return x, such that f (x) ≤ f (x * ) + and x ∈ E f ( optimal)

Algorithm's Structure and Annotations

Figure 8 presents an annotated version of the implementation. Annotations combine ellipsoidal constraints inEllipsoid and volume related properties to express the progress of the algorithm along iterates.

TOWARD PROVING ANNOTATION

FramaC supports the analysis of C code with ACSL annotations through multiple analyzers. The tool WP projects both the ACSL specification -the Hoare triples -and the source code in a predicate encoding on which external tools can be applied.

Low level matrix operation

We used this toolchain to prove that the annotations are valid. The goal is to have as many annotations as possible proven by SMT solvers since a Coq proof requires additional work. One of the difficulties in Coq is to address low level C manipulations such as matrix operations. The memory model axiomatization makes the manipulation of these encodings difficult while the properties are not specifically hard to prove. To avoid this, all matrix operations in C are associated to an ACSL contract, stating the computation performed in linear algebra instead of arrays and pointers. Fig. 10 presents an example of these annotations while Sec. 8 describes the (simplified) generated goal. SMT solvers easily discard this kind of proof objectives.

Lemmas hierarchy

Once these matrix operations have been abstracted by their ACSL axiomatization, all that remains is to prove the regular annotation.

/* @ ensures MatVar ( X, 2 , 1) == \old ( @ mat_add ( MatVar ( X, 2 , 1) , @ MatVar ( dx, 2 , 1) For each postcondition Q, we find what precondition P 1 , ..P k is required to prove Q. We also extract all variables v 1 , ..v l appearing in Q and P 1 , ..P k . With all this information we write an ACSL lemma stating : ∀v 1 , ..v l , P 1 ⇒ P 2 ⇒ ... ⇒ P k ⇒ Q.

) ) ; */ { X [0] = X [0] + dx [0]; X [1] = X [1] + dx [1]; }
Assuming this lemma an SMT will prove easily the postcondition. This permits to remove all C code aspect, to have pure math results to proof.

To ease the proof of this lemmas, one can build intermediary lemmas. This lemmas premits to split the proof into small pieces. Then each pieces can be proved by a SMT solver. We represented as a dependance graph in Fig. 5.2 the current status of the proof. This depency graph is used to tell the SMT which lemmas is useless to prove a given lemmas. This permits to remove a lot of lemmas from the context of the SMT simplifying their task.

We considered some results as Axioms(Blue colored in the figure). Deciding if a results should be an axiom is difficult. For now we decided that all math results we can find in (Nesterov, 2004, S 4 ) should be axioms. This permits to be clear on what we admited or not. Moreover results in Nesterov are usually quite general but remain feasible for SMT.

Maybe at some point, SMT wont be powerfull enough to prove the lemmas. If this happens we can write the proof by hand using a proof assistant as Coq. An other possibility would be to use machine learning to build the proofs, see e.g. [START_REF] Kaliszyk | Learning-assisted theorem proving with millions of lemmas[END_REF].

FLOATING-POINT CONSIDERATIONS

Another issue concerns the necessary rounding on the variables that happens due to the floating point representation of the variable and the impact of this numerical imprecision on the validity of the algorithm implementation.

In the following, we denote by F the set of floating point numbers. Since the sum of two floats is not necessarily a float, each computation will generate noise. One can still bound such numerical errors using the so-called standard model [START_REF] Rump | Verification of positive definiteness[END_REF]; [START_REF] Roux | A generic ellipsoid abstract domain for linear time invariant systems[END_REF]. This analysis can either be performed symbolically, analyzing the compuation performed and summing up the accumulated errors. Another approach, eg. implemented in the tool Fluctuat [START_REF] Delmas | Towards an Industrial Use of FLUCTUAT on Safety-Critical Avionics Software[END_REF], will rely on affine or interval arithmetics to represent such acucmulated errors of the standard models and bound the overall error obtained.

Problem Formulation: Ellipsoid method In the Ellipsoid method algorithm, we iterate the center point and matrix x, P such that the ellipsoid Ell(x, P ) has a decreasing volume and includes the optimal point. The operations being done in F and not R, we have no guarantee that the optimal point lies within the ellipsoid at each iteration (see figure 12). Therefore, this last condition being an important aspect of the program's semantics, the errors due to floating point rounding could have a direct impact on its stability.

Let fl(x) and fl(P ) be, respectively, the floating point values of x and P computed and used in the algorithm. x and P being the values in R obtained if the computations were not rounded. The objective here is to find a λ > 1 such that: Ell(x, P ) ⊂ Ell(fl(x), λ • fl(P ))

By doing that, we widen the ellipsoid Ell(x, P ), as illustrated on Fig. 12, to take into account the rounding.

The convergence issue becomes a problem when reaching the optimal value, ie. the need to prove the decreasing volume of the ellipsoids. Thus, the widening coefficient λ cannot be too big. This last constraint is currenlty under investigation.

CONCLUSION

In this article, annotations for numerical algorithms solving linear programming problems were proposed. We focused on a primal Interior Point algorithm and on the Ellipsoid method. The Interior Point algorithm is fully automated annotated while annotations for the Ellipsoid method are still manual. Annotations can be partitioned in two levels:

• Low-level annotation specifying matrix operations. This separates the handling of memory array accesses from math reasoning, easing the automatic proof of simpler annotations by SMT solvers.

• Math reasoning annotations specify the high level properties of the algorithms: convergence, preservation of feasibility, -optimality of the solution, etc. These annotations are proven either by SMT solvers or with Coq. The finalization of these proofs is still a work in progress.

Our perspectives include the application of this approach to more general settings, eg. generic linear programming problems. We plan to couple code generation with annotations and proof synthesis, providing embedded code together with its functional soundness proof, at the code level. This would include the generation of the corresponding Coq proofs.

Our approach is also compatible with more general convex optimization problems. Finally, we intend to address floating point semantics issues while performing the proofs. This would guarantee a completely sound and bug-free implementation.

Fig. 2 .

 2 Fig. 2. Hessian norm in ACSL. optimization, defining norm, grad, hess, lower.... Figure 2 presents the definition a Hessian-induced norm, parametrized by the Hessian at point X.

Fig. 3 .

 3 Fig. 3. In blue: Points of the central path for different value of t Central path. Using both the barrier function and the initial cost function, one can characterize their linear combination. Definition 5. Let us introduce a new function, f as:

Fig. 7 .

 7 Fig. 7. Ellipsoid Cut for an iteration applied to LP Definition 10. (Volume of ellipsoids). Let E(P, x) be an ellipsoid set in R n . We denote by V ol(E) its volume defined as Vol(E(P, x)) = Vol(E) = det(P ) • Vol(V n )

  Fig. 8. Shape of the Main Ellipsoid Method Algorithm

Fig. 10 .

 10 Fig. 10. Simple example of matrix operation specification

  Fig. 11. Lemma proof hierarchy

By construction, thanks to the use of a self-concordant barrier function.

ACKNOWLEDGEMENTS

This work was partially supported by projects ANR ASTRID VORACE and NSF CPS SORTIES under grant 1446758. The authors would also like to deeply thank Didier Henrion for his participation to this work.


Luca Di Stasio
email: luca.di-stasio@univ-lorraine.fr

Janis Varna
email: bjanis.varna@ltu.se

Zoubir Ayadi
email: czoubir.ayadi@univ-lorraine.fr

EFFECT OF BOUNDARY CONDITIONS ON MICRODAMAGE INITIATION IN THIN PLY COMPOSITE LAMINATES

Keywords: Polymer-matrix Composites (PMCs), Thin-ply, Transverse Failure, Debonding, Finite Element Analysis (FEA)

.

Introduction

Thin ply composites represent today one of the most promising material for advanced applications in the aerospace industry and are attracting the interest of structural designers for use in mission-critical applications such as cryogenic tanks [START_REF] Jackson | Composite Cryotank Technologies and Development 2.4 and 5.5M out of Autoclave Tank Test Results[END_REF] and reusable space launchers' frames [START_REF] Kopp | The Aurora Space Launcher Concept[END_REF].

The experimental evidence on thin-ply laminates [START_REF] Sihn | Experimental studies of thin-ply laminated composites[END_REF][START_REF] Saito | Experimental Evaluation of the Damage Growth Restraining in 90° Layer of Thin-ply CFRP Cross-ply Laminates[END_REF][START_REF] Arteiro | Notched response of non-crimp thin-ply laminates[END_REF][START_REF] Amacher | Thin ply composites: Experimental characterization and modeling of size-effects[END_REF] confirms an earlier result on the influence of thickness and lay-up sequence on the strength and crack suppression behavior of FRP laminates, namely the existence of the so-called thin-ply effect and in-situ strength. At the end of 1970's, the experimental investigation of Bailey and collaborators [START_REF] Garrett | Multiple transverse fracture in 90 cross-ply laminates of a glass fibrereinforced polyester[END_REF][START_REF] Parvizi | On multiple transverse cracking in glass fibre epoxy cross-ply laminates[END_REF] pointed out the relationship between ply thickness and transverse cracking, by showing that reducing the thickness of an inner 90°-ply in a cross-ply laminate results in the delay and even suppression of transverse crack propagation [START_REF] Parvizi | Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates[END_REF]. They identified the fiber-matrix property mismatch, and in particular the difference in Poisson's ratio, as the main mechanical driver of matrix cracks [START_REF] Bailey | On the transverse cracking and longitudinal splitting behaviour of glass and carbon fibre reinforced epoxy cross ply laminates and the effect of Poisson and thermally generated strain[END_REF]. By investigating the micromechanics of transverse cracks, they observed that they originate as debonds or flaws at the fiber-matrix interface, which under increasing loads extend initially along the interface itself, then kink, coalesce and propagate as transverse cracks first through the thickness and finally through the width of the specimen [START_REF] Bailey | On fibre debonding effects and the mechanism of transverse-ply failure in cross-ply laminates of glass fibre/thermoset composites[END_REF].

Experimental evidence thus points to debonding at the fiber/matrix interface as the primary mechanism to investigate in order to achieve a better understanding about the initiation of transverse cracking and the potential of its suppression through an improved laminate design. It is paramount to this end to understand the process of fiber/matrix debonding, kinking and coalescence, as well as the effect of ply thickness, non-uniformity of fiber distribution, material properties' mismatch and thermal strains. Cohesive zone modeling (CZM) has been one of the numerical strategies of choice to model these phenomena, often coupled with an elasto-plastic behavior for the polymeric matrix [START_REF] Kushch | Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading[END_REF][START_REF] Herraez | Transverse cracking of cross-ply laminates: A computational micromechanics perspective[END_REF]. The method provides the possibility of simulating the onset and propagation of transverse cracks starting from a virgin ply, and can be used to record virtual crack path and loading histories which can then be compared with macroscopic experimental tests [START_REF] Canal | Intraply fracture of fiber-reinforced composites: Microscopic mechanisms and modeling[END_REF]. However, several drawbacks exist that cast a serious doubt about its applicability. First and foremost, the bi-and tri-axial stress state observed in inter-fiber regions and the cavitation-like matrix failure taking place at or close to the fiber interface [START_REF] Asp | Effects of a composite-like stress state on the fracture of epoxies[END_REF][START_REF] Asp | A criterion for crack initiation in glassy polymers subjected to a composite-like stress state[END_REF][START_REF] Asp | Prediction of matrix-initiated transverse failure in polymer composites[END_REF] do not agree with the adhesive-like failure mechanism described by the CZM.

A second approach to the problem is the evaluation of the strain and stress state around the debond and consequent estimation of the Energy Release Rate (ERR) at the crack tip in mode I and mode II through the application of the Virtual Crack Closure Technique (VCCT) [START_REF] Krueger | Virtual crack closure technique: History, approach, and applications[END_REF] and/or the J-integral method [START_REF] Rice | A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[END_REF]. Its strength lies in the ability to provide direct comparisons of the magnitude of the effect on crack onset of different geometric, mechanical and thermal factors. However, unless coupled with a failure criterion as in [START_REF] Mantic | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF], it has no predictive power regarding the direction of crack propagation and the size of the propagated crack. This approach was firstly used in the analytical treatment of the fiber/matrix interface problem in FRPC by Toya [START_REF] Toya | A crack along the interface of a circular inclusion embedded in an infinite solid[END_REF], and more recently in numerical analyses using the Boundary Element Method (BEM) [START_REF] París | The fiber-matrix interface crack -a numerical analysis using boundary elements[END_REF] or the Finite Element Method (FEM) [START_REF] Varna | Growth and interaction of debonds in local clusters of fibers in unidirectional composites during transverse loading[END_REF].

The recent advancements in thin-ply technology [START_REF] Kopp | The Aurora Space Launcher Concept[END_REF], where ply thickness can now be reduced to only 2-4 times the diameter of the reinforcement, has given a novel practical relevance to the theoretical analysis of advantages and disadvantages of using thin plies to delay intralaminar cracking. In this work we will compare the ERR related to the fiber/matrix debond growth for a debonded fiber placed in three different locations: a) at the free surface of the UD layer; b) at the interface with a very stiff layer of different orientation; c) as a repeating element in a thick layer. To magnify the expected trends, in cases a) and b) we analyze an extremely thin ply with only one fiber in the thickness direction. The effect of fiber content on fiber/matrix debonding is analyzed through the evaluation of Mode I and Mode II ERRs at the debond's crack tip using the Finite Element Method in conjunction with the VCCT and the J-integral method. Although extreme in its simplification, the model provides valuable insights of the constraining effect of adjacent plies and on the effect of fiber volume fraction, and thus on ply thickness.

The Model

In order to investigate the initiation mechanisms of transverse cracks, we focus our attention on the fiber/matrix interface and we apply the tools of Linear Elastic Fracture Mechanics (LEFM). We consider a 2-dimensional Finite Element (FE) model of a Representative Volume Element (RVE) consisting of a single row of glass fibers in epoxy matrix.

Several fibers have debonds of the same size at the interface with the matrix, indicated in red in Figures 1234. An interaction between debonds is expected at high fiber content (especially between the closest ones), but the interaction will be the same for all the three analyzed cases of the top and bottom surfaces' constraints, which is the main subject for comparison in this paper. The selected configuration of debonds allows the selection of a Repeating Unit Cell (RUC) on which different surface constraints will be applied. The first of such configurations, depicted in Figure 1, represents the case of a single ultra-thin ply with only one layer of fibers in the thickness direction. The bottom and top surfaces are traction-free: free surface in the FEM model (case Free in results). The second layup considered (Figure 2) corresponds to a single inner ultra-thin 90°-layer bounded by stiffer and thicker 0° plies, which is not bending and has a constant strain in the x-direction (the inhomogeneity and damage in the 90°-layer has no effect on strain in the 0°-layer). It is modeled by assigning a constant unknown vertical displacement to account for the Poisson's effect and a linear distribution of x-displacement on the top surface (case Constant v, linear u in results). The final configuration, seen in Figure 3, represents a thick ply with fibers in square-packing arrangement in both the width and thickness direction. This last case is modeled by the application of a constant unknown vertical displacement on the top surface (case Constant v in results).

Pure plane strain conditions in the y-direction are considered and both the glass fiber and the epoxy are assumed to be isotropic materials with properties reported in Table 1. An arc crack is placed at the fiber/matrix interface, referred to as the debond in this article, where the matrix is debonded from the fiber (represented in Figure 4 in the deformed state). The debond is characterized by its angular size, equal to 2∆𝜃 and indicated in blue in Figure 4, and its position 𝜃 , taken as the angular position of debond's mid point. The debond's size is a free parameter, while its position 𝜃 is in this paper assumed equal to 0, i.e. a centered debond symmetric with respect to the xaxis. At large sizes, the debond presents a so-called contact zone of angular size ∆Φ close to the crack tip, where the crack faces are in contact and only Crack Sliding Displacement (CSD) is present. Only half of the RUC is considered and symmetry conditions are applied along the x-axis to model the presence of the lower half. Uniaxial tension is applied to the RUC in the form of an applied constant strain of 1% in the longitudinal direction. integral method [43]. For the latter, the Abaqus built-in function is used; for the VCCT, an in-house Python routine has been developed. The model for a single fiber in an infinite matrix (𝑉 = 0.0079%) is validated against the results obtained in [START_REF]Analysis user's manual[END_REF] using the BEM in the case of a single fiber.

Results & Discussion

As seen in Figure 5, the Mode I Energy Release Rate shows a strong dependence on fiber content. The peak value at 𝑉 = 60% is approximately 4 times the value at 20% and it is shifted towards lower debond's sizes. For ∆𝜃 > 60 -70°, Mode I ERR decays to zero signaling the emergence of the contact zone (see also Figure 7), where Crack Sliding Displacement (CSD) is the only component present. The effect of boundary conditions is less pronounced, although present. In particular, the free surface case presents always lower values of 𝐺 due to the ability of the matrix to comply locally, ability that

G II [J/m 2 ] Δθ [°] Mode II ERR -VF f =60% Free Constant v(x,L) Constant v(x,L), linear u(x,L)
Luca Di Stasio, Janis Varna and Zoubir Ayadi is prevented when the lateral movement of the surface is constrained by the adjacent material. This is an interesting result, because one would expect that closeness of the free specimen's surface would facilitate the debond growth. For the RUC (thick 90°-ply) at high fiber content, the value of 𝐺 is higher as it reflects the presence of another identically debonded fiber on top of the analyzed (Figure 3), which interact with each other at higher 𝑉 (i.e. lower RUC's thicknesses).

The mode II ERR (Figure 6) is on the other hand less sensitive to the change in fiber content, and thus on the presence of adjacent material. Changes in value do not exceed 15 -20% for any combination of fiber content and boundary conditions. As for mode I, increasing the fiber content enhances the effect of the adjacent material on crack initiation. When only the vertical displacement is constrained, mode II ERR is reduced by 15 -20%; when the horizontal displacement is also imposed, a similar reduction in magnitude occurs but a longer plateau of maxima appears. The presence of a free surface favors Mode II, as it is shown by the approximately 20% increase in the value of 𝐺 for a thinner ply. The observed results can be explained in geometric and mechanical terms: under transverse tension, the matrix will contract toward the symmetry line more than the fiber due to the mismatch in Poisson's ratio. When the surface is free, for ∆𝜃 > 45 the applied tension generates a radial Crack Opening Displacement (COD) equal to or smaller than the one created by the Poisson effect (opposite in sign, i.e. closing the crack), thus leaving the crack closed. The tangential components are instead coherent and sum up favoring a mode II crack initiation. When the vertical displacement is constrained, the COD component due to the Poisson's effect is reduced, thus favoring mode I at higher debond's size. The transition to a closed crack with a contact zone is thus delayed, but after the onset of the latter the behavior is unchanged. At larger debond's sizes, in fact, the main driver of crack sliding is the applied tension and not the Poisson effect, which induces displacements approximately normal to the interface and consequently does not have a relevant tangential component anymore. When an horizontal displacement is imposed to the RUC's surface however, it reduces at higher debond's sizes (∆𝜃 > 60°) the constraining effect of the bonded fiber/matrix interface in the near-tip region by forcing the matrix to follow the average applied tensile displacement. This, in turn, translates in a higher local tangential component of the crack faces displacement, which corresponds to higher values of 𝐺 with respect to the other two combinations of boundary conditions. This favors mode II and creates the plateau seen in the picture on the right in Figure 6. The effect of fiber content is on the other hand explained by the de Saint Venant's principle: the lower the fiber content, the bigger the RUC and the less the debond is affected by perturbations in the elastic fields caused by the conditions at the boundary. This mechanical explanation is confirmed by the evolution of the contact zone size (Figure 7). For low fiber content, the debond is not significantly affected by perturbations caused by boundary conditions and its evolution is thus independent from the latter. For high 𝑉 , the constraing effect of adjacent material causes a delay of 30 -40° in contact zone onset and a decrease in size of approximately 40°.

The numerical results presented agree and explain the microscopic observations of Saito & al. [START_REF] Saito | Experimental Evaluation of the Damage Growth Restraining in 90° Layer of Thin-ply CFRP Cross-ply Laminates[END_REF].

Working with [02/90n/02] laminate, with n = 1,2 and 4, they observed that thin plies favored the onset of debonds and their opening, which happened at strains lower than for thick plies. Coalescence of debonds and through-the-thickness propagation were instead delayed and even suppressed in the thinnest laminae.

Conclusions

The initiation mechanisms at the fiber/matrix interface in thin-ply laminates under transverse tensile strain have been investigated evaluating Mode I and Mode II ERRs and the contact zone's size. Three different cases have been considered: a) a debonded fiber close to the specimen free surface; b) a debonded fiber in the middle of a thick 90°-layer; c) a debonded fiber in the proximity of the interface with a stiff ply. These three configurations are modeled by different combinations of boundary conditions: case a) corresponds to a free surface in FEM model; in case b) a constant unknown vertical displacement is applied to upper surface; case c) is modeled by the concurrent application of a constant unknown vertical displacement and a linear horizontal displacement. Several values of fiber volume fraction have been evaluated, and the two extreme cases of 𝑉 = 20% and 𝑉 = 60% have been reported and analyzed to exemplify the trends found.

In accord with de Saint Venant's principle, it has been shown that the main consequence of ply thickness, controlled in the model by the fiber content, is to enhance the effect of the conditions applied to the boundary. This translates in a strong increase in the magnitude of mode I ERR, which is dominated by the applied uniaxial tensile displacement. The presence of a stiff material adjacent to a thin ply reduces mode II ERR, delays the onset of the contact zone and drastically reduces the size of the latter. The increase in mode I ERR with the concurrent decrease in mode II ERR and contact zone size for a constrained thin ply means that opening of the debond is favored for the same level of strain. This result agrees well and provides a mechanical explanation to the microscopic observations reported in [START_REF] Saito | Experimental Evaluation of the Damage Growth Restraining in 90° Layer of Thin-ply CFRP Cross-ply Laminates[END_REF], where it was noted that fiber/matrix debonds in thin plies appeared at lower strain levels than in thick plies.

Figure 1 .

 1 Figure 1. Schematic representation of a single "free" 90° ply with debonded fibers.

Figure 2 .

 2 Figure 2. Schematic representation of the embedded 90° ply bounded by thick UD sections.

Figure 3 .Figure 4 .

 34 Figure 3. Schematic representation of the repeating element in the thick 90° ply.

Figure 5 .

 5 Figure 5. Evolution of 𝐺 with respect to debond's size for the 3 models considered with 20% (left) and 60% (right) fiber content.

Figure 6 .

 6 Figure 6. Evolution of 𝐺 with respect to the debond's size for the 3 models considered with 20%(left) and 60% (right) fiber content.

Figure 7 .

 7 Figure 7. Evolution of contact zone size ∆Φ with respect to debond's size for the 3 models considered with 20% (left) and 60% (right) fiber content.

  Luca Di Stasio, Janis Varna and Zoubir Ayadi

Table 1 .

 1 Material properties for glass fiber and epoxy.

	Material	E [GPa]	ν[-]
	Glass fiber	70.0	0.2
	Epoxy	3.5	0.4

Luca Di Stasio, Janis Varna and Zoubir Ayadi

Luca Di Stasio, Janis Varna and Zoubir Ayadi

Acknowledgments

Luca Di Stasio gratefully acknowledges the support of the European School of Materials (EUSMAT) through the DocMASE Doctoral Program and the European Commission through the Erasmus Mundus Program.