Nazim Fatès
email: nazim.fates@inria.fr

Remarks on the cellular automaton global synchronisation problem

Keywords: inverse problems, SAT solving, stochastic vs. deterministic solutions

The global synchronisation problem consists in making a cellular automaton converge to a homogeneous blinking state from any initial condition. We here study this inverse problem for one-dimensional binary systems with periodic boundary conditions (i.e., rings). For small neighbourhoods, we present results obtained with the formulation of the problem as a SAT problem and the use of SAT solvers. Our observations suggest that it is not possible to solve this problem perfectly with deterministic systems. In contrast, the problem can easily be solved with stochastic rules.

Introduction

The study of inverse problems is becoming a fertile field in the research on cellular automata (CA). Among the recent achievements, we mention the construction of an exact solution to the parity problem [?], the construction of exact or approached solutions on the density classification problem [?,?,?,?] or the reinterpretation of the solution to the Firing squad problem with fields [?].

A common feature of these problems is the need to reach a global consensus: there exists a moment where all cells, or a large fraction, must agree on a given state that is interpreted as the output of the algorithm. The difficulty is related to the propagation of this information from a local to a global scale: in a decentralised framework such as cellular automata, how do the cells "agree" on a common state while they have only a local view of the system?

We here study the global synchronisation problem: in its original form, the question is to find a CA rule such that, from any initial condition, the system reaches a "blinking state" in which the two homogeneous configurations alternate. This problem can be generalised to more states but we will here restrict our study to the binary case.

Since its formulation by Das et al. in 1994 [?], the problem has received only a limited attention. This lack of interest is probably due to the fact that it is much easier to solve than other inverse problems such as the density classification problem. In fact, solutions with "100% success rates" were presented in the very paper where the problem was formulated [?]. The authors used genetically engineered solutions to show that it was possible to obtain a performance of "100%" for ring sizes going up to 999. Their interest was to find rules which attain a consensus by removing the "defects" that separate the non-synchronised regions of the system.

Our purpose in this note is to go one step forward and to ask if perfect solutions do exist. We will thus request that all initial conditions lead to the blinking cycle and not only a sample of configurations, drawn at random. After presenting the formal definitions of the problem (Sec. ??), we will present a simple "manual proof" that no ECA solves the problem (Sec. ??). This construction will guide us to formulate the problem as a SAT problem (Sec. ??) and to obtain a first set of results for larger neighbourhoods (Sec. ??). We then show that, in contrast, perfect solutions can easily been constructed (Sec. ??). We conclude by formulating a few questions.

Fundamentals

CA definitions

We here consider finite binary cellular automata with periodic boundary conditions. The basic components of our systems are the cells; each cell can hold one of the two states: 0 or 1. The variable n is used to denote the number of cells that compose the system, the cells are arranged in a ring and the set of cells is denoted by L = Z/nZ.

A configuration x = (x i ) i∈L represents the state of the system at a given time. The set of configurations is denoted by E n = {0, 1} L . The interactions between the cells are local, that is, each cell can only "see" a finite subset of the cells of the system, the neighbourhood. Without loss of generality, we can consider that the neighbourhood N of our cellular automata are formed of discrete intervals: N = {-l, . . . , r}, with l ∈ N and r ∈ N * = {1, 2, . . . }. The width of this interval is called the size of the neighbourhood and is denoted by k, with k = l + r + 1.

The evolution of a cell follows a function f : {0, 1} k → {0, 1}, called the local rule. Following Wolfram's notation, a local rule of size k, that is, defined on a neighbourhood of size k, is assigned a code W which is an integer between 0 and 2 k -1. This code is given by the formula:

W = 2 k -1 i=0 f (b k (i), . . . , b 1 (i)) • 2 i
where b j (i) is the value of the j-th bit of the binary representation of i.

For a given ring size n, the global transition function F : E n → E n associated to the ring size n is the function that maps a configuration x t to a configuration x t+1 = F (x t ) such that x 0 = x and:

∀i ∈ L, x t+1 i = f (x i-l , . . . , x i+r ).
(1)

Note that to be perfectly rigorous, we should denote F with indices showing that it depends on f and n. We however drop these elements for the sake of clarity since f and n will be made clear from the context.

Formulation of the problem

We denote by 0 = 0 L and 1 = 1 L the two uniform configurations. In the following, we will require that 0 and 1 are two non-quiescent states, that is, f (0, . . . , 0) = 1 and f (1, . . . , 1) = 0. We call this condition the blinking condition.

We define the height h of a configuration x ∈ E n as the time needed to reach one of the two uniform configurations:

h(x) = min{t ∈ N, x t = F t (x) ∈ {0, 1}}.
(2)

with the convention that h(x) = ∞ if x t does not reach 0 or 1. Similarly, the height of a given configuration space E n is the maximum height of the configurations of E n . We say that F synchronises a configuration x ∈ E n if h(x) is finite. We also say that x is synchronised on 0 (resp. on 1) if the first homogeneous configuration that is met is 0 (resp. 1). Similarly, we say that F synchronises the size n if it synchronises all the configurations of E n . We can now formulate the global synchronisation problem: Does there exist a local rule f such that for all n ∈ N * , the associated global function F synchronises the size n?

Elementary properties

Let f be a local rule of size k. The reflexion R(f ) and the conjugate C(f ) are the local rules respectively obtained by the exchange of the left and right directions and by the exchange of the 0 and 1 states. Formally, ∀(q 1 , . . . , q k ) ∈ {0, 1} k , R(f )(q 1 , . . . , q k ) = f (q k , . . . , q 1 ), and C(f )(q 1 , . . . , q k ) = f (q 1 , . . . , q k ), where q = 1 -q denotes the inversion of states. Similarly, the reflexion-conjugate rule RC(f ) is the local rule obtained by composing the two previous symmetries:

RC(f ) = R • C(f ) = C • R(f ).
Proposition 1 (rule symmetries). f is a solution to the global synchronisation problem if and only if R(f ), C(f ), RC(f ) are also solutions.

Proof. Clearly, the property of synchronising a given size is preserved by the reflection and conjugation symmetries: for a given size n, if F, F r , F c , F rc are the global functions respectively associated to f, R(f ), C(f ) and RC(f ), then F synchronises the size n if and only if F r , F c , F rc synchronise the size n.

Let σ denote the (left) shift operator, that is, a function σ : E n → E n such that ∀i ∈ E, σ(x) i = x i+1 . We call the rotations of x the set of configurations that are obtained by applying a positive number of shifts on x; this set is denoted by

[x] = {σ k (x), k ∈ N}.

Proposition 2 (configuration symmetries). A global rule F synchronises a configuration x if and only if it synchronises all the configurations of [x].

This simply results from the fact that: (a) F commutes with the shift and (b) [0] = {0} and [1] = {1}. Note that the rotation [•] defines an equivalence class: we say that a configuration y is equivalent to x if y is a rotation of x. It can be easily verified that this is an equivalence relation.

The next proposition states that the iterates of a configuration can not be contained in the rotations of this configuration.

Proposition 3 (images of a configuration). If F synchronises a configura- tion x, then ∪ k≥1 F k ([x]) [x] = ∅, that is, ∀x ∈ [x], ∀k ∈ N * , F k (x ) / ∈ [x].
Proof. By contradiction, let us assume that there exists k ∈ N * and i ∈ N such that F k (x) = σ i (x). By recurrence, using the commutation of the shift with F , we have: F kn (x) = σ in (x). Since the space is a ring of size n, we have σ in (x) = x, which implies x = F kn (x). The configuration x thus evolves on a cycle of length k. The two homogeneous states 0 and 1 are excluded from this cycle -otherwise x would be reachable from these two states -and x can not be synchronised, which contradicts the hypothesis.

An immediate consequence of this proposition is that if F synchronises x, then x can not be a fixed point

(F (x) = x) or a blinking point (F (x) = x and F (x) = x) or a translating point (F (x) = σ i x with i ∈ N * ). Proposition 4 (color-discernation). If a local rule f is a solution to the problem, then it is not color-blind, that is, C(f ) = f .
Proof. By contradiction, let us assume that we have C(f ) = f . Let us take an even size n = 2m with m ∈ N and consider the configuration x = (01) m . Without loss of generality we can assume that x is synchronised on 0. Formally, if we denote by T = h(x) the height of x, this reads:

F T (x) = 0.
Then, if we consider x = (x i ) i∈L , we have:

F T c (x) = F T (x) = 0 = 1.
On the other hand: F T (σx) = σF T (x) = σ0 = 0. Since x = σx, we remark that these two equations are contradicting.

Direct inspection of the ECA space

As a starting point, let us consider Elementary Cellular Automata (ECA), that is, binary CA with N = {-1, 0, 1}.

Proposition 5. There exists no ECA which solves the synchronisation problem.

Proof. We have 256 rules to consider, these rule are defined with:

a = f (0, 0, 0) b = f (0, 0, 1) c = f (1, 0, 0) d = f (1, 0, 1) f = f (0, 1, 1) e = f (0, 1, 0) g = f (1, 1, 0) h = f (1, 1, 1
) . We now have to determine whether there exists an assignment of these eight boolean variables a, . . . , h which satisfies the problem.

Case n = 1: Given the specification of the problem, the states 0 and 1 are not quiescent: a = 1 and h = 0. We are thus left with 64 rules to search.

Case n = 2: Proposition ?? implies d = e. Indeed, if the two variables are not equal, 01 would be either a fixed point or a blinking point. We are now left with 32 rules.

Case n = 3: By noting that the configuration 001 can not be a fixed point and can not be translated (by Prop. ??), we obtain:

s 1 = b + c + e = 1 (CondA).
By symmetry, we have

s 2 = d + f + g = 2 (CondB).
The case where a configuration of [001] is transformed into a configuration of [011] and vice-versa is also impossible, otherwise the uniform configurations would never be reached (Prop. ??). We thus have (s 1 , s 2 ) = (2, 1) (CondC). It is then easy to verify that these conditions are sufficient to achieve the synchronisation of size 3.

We find the following set of remaining 8 rules, divided into three sets of rules: (1, 127), (9, 65, 111, 125), (19, 55). The rules are grouped by their equivalence with the conjugation and reflexion symmetries (see Prop. ??).

Case n = 4: It is easy to check that none of the remaining rules allows a synchronisation of size 4. For instance, we remark that 0001 is a translating point for rule 1 and 9 and that 0110 is a translating point for rule 19.

The synchronisation problem as a SAT problem

The previous method required a methodical inspection of the space but this approach can not easily been generalised to larger neighbourhoods. We now propose to do an "automated filtering" of the rules by transforming the problem into a SAT problem. The idea is to model the transition table as sequence of boolean variables and to express a boolean formula to state that a configuration is synchronised. We want to examine more and more initial conditions until we reach a point where we find that the problem is not satisfiable. If, on the contrary, we find that there exists a rule which synchronises all the initial conditions considered, then we would have a good "candidate" for solving the problem.

General SAT formulation

A SAT problem consists in finding a assignment to boolean variables that satisfies a given boolean formula. A clause is a conjunction of literals, that is, a boolean formula defined only with the or operator (e.g., a ∨ b ∨ c). The convention is that a SAT problem is formulated in a conjunctive normal form (CNF): it is a disjunction of clauses. In the sequel, we call a CNF formula any of such disjunction of clauses.

Let k be the size of the neighbourhood and M = 2 k the number of transitions1 of this neighbourhood. We introduce M boolean variables t 0 , . . . , t M -1 to encode the transitions of a rule f ; the convention is that t i is true if and only if f b k (i), . . . , b 1 (i) = 1. For an initial condition x ∈ E n , we also introduce (τ + 1) • n additional variables, denoted by (ξ(t, i)) t∈{0,...,τ },i∈{0,...,n-1} , which correspond to the values (x t i ) ∈ {0, 1} taken by the cells in the evolution of x.

Blinking condition

The blinking condition f (q, . . . , q) = q is simply expressed by a CNF formula with two atomic clauses: F bl = t 0 ∧ ¬t M -1 .

Initial state

Let us now see how to encode the states of an initial condition x ∈ E n in a formula. The operation simply consists in "translating" the initial condition x into the CNF formula:

F ic (x) = i∈{0,...,n-1} 1{ξ(0, i), x 0 i }. ( 3
)
where 1{V, q} is a function which associates to the boolean variable V and to a cell state q ∈ {0, 1} the variable V if q = 1 and the variable ¬V otherwise.

For instance if we have x = 0011 as an initial condition, the associated formula will be:

F ic (x) = ¬ξ(0, 0) ∧ ¬ξ(0, 1) ∧ ξ(0, 2) ∧ ξ(0, 3).
(4)

Synchronisation condition

Let τ be the maximum number of time steps to achieve the synchronisation. To express the condition that x is synchronised in at most τ steps, we can write x τ ∈ {0, 1}. Unfortunately, this can not be translated in a straightforward way into a CNF formula. Indeed, if we write:

F sc (x, τ ) =   i∈{0,...,n-1} ξ(τ, i)     i∈{0,...,n-1} ¬ξ(τ, i)   , (5)
we need to distribute the and operator over the orto obtain a CNF. This is why we prefer to formulate this condition as:

F sc (x, τ ) = i∈{0,...,n-2} ξ(τ, i) = ξ(τ, i + 1), (6)
which simply represents the fact that all the states of x τ are equal. By noting that a = b is equivalent to (a ∨ ¬b) ∧ (¬a ∨ b), F sc becomes:

F sc (x, τ ) = i∈{0,...,n-2} (ξ(τ, i) ∨ ¬ξ(τ, i + 1)) ∧ (¬ξ(τ, i) ∨ ξ(τ, i + 1)). (7)

Consistency conditions

We now need to write a CNF formula for the condition: x t+1 = F (x t ) for t ∈ {0, . . . , τ -1}. We call this formulation the "consistency condition", as it expresses the fact a given boolean formula in consistent with the evolution of the cellular automaton. We now give a precise description of this CNF formula.

In order to ease the notations, let us detail this operation for the specific case of N = {-1, 0, 1}; it is easy to generalise it to other neighbourhoods. Locally, our condition is expressed by

∀i ∈ L, x t+1 i = f (x t i-1 , x t i , x t i+1 ) (8)
which is translated as:

∀i ∈ {0, . . . , n -1}, ϕ(ξ(t + 1, i), ξ(t, i -), ξ(t, i), ξ(t, i + )) (9)
where i -= (i -1) mod n and i + = (i + 1) mod n, and where ϕ is a function that remains to be found. Our goal is to find ϕ such that ϕ(y , x, y, z) is a CNF formula that expresses that y is the result of transition f (x, y, z). In a usual programming environment, one would need simply to calculate i = x + 2y + 4z and then to read the value of t i and assign it to y . Unfortunately, there is no direct way of "coding" these operations in a SAT formula. We thus need to enumerate all the possible values for the variables y , x, y, z and then write a consistency condition that expresses that y equals t i where i is the index which corresponds to the transition (x, y, z).

For example, if we take transition f (0, 1, 1) = 1, we have i = 3 and we write the formula with five clauses:

ϕ 3 = ¬x ∧ y ∧ z ∧ t 3 ∧ (y ∨ ¬t 3 ) ∧ (¬y ∨ t 3 ), ( 10
)
where the two last clauses stand for y = t 3 . Formally, we write: ϕ(y , x, y, z) = λ∈{0,...,7} ϕ λ with:

ϕ λ = 1{x, λ 1 } ∧ 1{y, λ 2 } ∧ 1{z, λ 3 }∧ 1{t λ , f (x, y, z)} ∧ (y ∨ ¬t λ ) ∧ (¬y ∨ t λ ), (11)
where λ i = b i (λ) is the value of the i-th bit of the binary representation of λ.

By distributing the and operator over the or, ϕ λ becomes:

ϕ λ = (1{x, λ 1 } ∨ 1{y, λ 2 } ∨ 1{z, λ 3 } ∨ y ∨ ¬t λ ) (1{x, λ 1 } ∨ 1{y, λ 2 } ∨ 1{z, λ 3 } ∨ ¬y ∨ t λ ). ( 12
)
Each elementary transition of a given cell at a given time step is thus encoded with a formula ϕ which contains 2M = 2 k+1 = 16 clauses. As there are nτ such elementary conditions, the consistency condition is given by F e with 2 k+1 • nτ clauses:

F e (x, τ ) = t∈{0,...,τ -1}, i∈{0,...,n} ϕ ξ(t + 1, i), ξ(t, i -), ξ(t, i), ξ(t, i + ) . (13)

Combining initial conditions

The last step that remains is to combine various initial conditions in order to: (a) either find out that the problem is not solvable for a given setting or (b) exhibit a good candidate to solve the problem. We proceed iteratively by increasing the size of the initial conditions to synchronise. From Prop. ??, we know that we do not need to consider all the initial conditions: for each size n, it is sufficient to select only one initial condition in each possible set of rotations. Formally, we say that a set of configurations is representative if the rotations of its members form a partition of the configuration space. Formally, let us denote by χ(n) the set of representative conditions; we write:

χ(n) = {X ⊂ E n , ∀x, y ∈ X, [x] ∩ [y] = ∅, x∈X [x] = E n }. ( 14
)
The following table shows the growth of these sets To sum up, if we fix a set of ring sizes S = {n 1 , . . . , n s }, and a maximum synchronisation time τ , we construct a sequence of sets of configurations: X 1 , . . . , X s such that X i ∈ χ(n i ) and build the general formula:

F synch (S, τ ) = F bl i∈{1,...,s} x∈Xi F ic (x, τ ) ∧ F sc (x, τ ) ∧ F e (x, τ ). ( 15
)
This is the final formula; it expresses the fact that all the configurations of size n ∈ S are synchronised in at most τ time steps.

First experimental results

The use of SAT solvers is a well-explored field of research in Computer science.

As we are not a specialist of these questions, we did not endeavour to optimise the search for a solution by any means. We simply used the minisat solver3 and generated the formulae with our cellular automata simulation program FiatLux4 .

ECA space

We take k = 3. The results with the SAT solver confirm the results of Sec. ??: there exists no rule which is a solution for S = {2, 3, 4}. For S = {2, 3}, we find that: -(19, 55) and (9, 65, 111, 125) have height of 2, and no rule has height of 3.

-(1, 127) have height of 1,
We can also explore the synchronisation for other sets of ring sizes S. For S = {4}, we find that:

-(37, 91) have a height of 3, -(25, 67, 103, 61) and (45, 101, 75, 89) have height of 4 and, no rule has a height of 5.

For ECA 61, the synchronisation process is presented on Fig. ??, p. ??. Surprisingly, for S = {5}, we find that (9,65,111,125) synchronises with a height of 5.

The k = 4 space

We now examine a neighbourhood with one more cell: we take k = 4 and N = {-1, 0, 1, 2}. The space contains 2 2 4 = 2 16 = 65536 rules.

By testing increasing values of n and setting the maximum synchronisation time τ equal to χ(n) (see Prop. ??), we found that the maximum synchronisation length of this neighbourhood is 6, that is, for S = {2, . . . , 7}, no solution is found.

For S = {2, . . . , 6}, we find that:

The k = 5 space

We now examine a neighbourhood with one more cell: we take k = 5 and N = {-2, -1, 0, 1, 2}. The space contains 2 2 5 = 2 32 ∼ 4.10 10 rules. At this point, we reach the limits of our approach: The CNF formula of the problem has 74768 variables and 4563060 clauses. By progressively increasing the values of n and τ , our best result was to find a rule which synchronises the size interval S = {2, . . . , 11}: rule 1078270911. This rule has a height of 18 (see Fig. ??, p. ??). There are probably other rules which solve the problem for τ higher than 18, but we leave this exploration for future work. For S = {2, . . . , 12}, no solution was found for τ = 30. Surprisingly, the non-satisfiability of the formula is given rapidly, which suggests that the maximum synchronisation length for k = 5 is equal to 12.

We also tested the SAT solver for k = 6, but τ = 12 is sufficient to generate SAT problems that are not solved after more than two hours of computation. We thus leave the exploration of these greater spaces for future work.

"Perfect" stochastic solutions

We now propose to examine what is the situation of the stochastic rules. In fact, if we allow randomness in the transitions of the rule, it becomes difficult not to solve the problem! For the sake of simplicity we restrict our study to the probabilistic ECA case. We thus take N = {-1, 0, 1} and define a local transition function φ : Q 3 → [0, 1], which associates to each neighbourhood state its probability to be updated to 1.

Formally, starting from a configuration x, the system can be described by a stochastic process (x t ) t∈N . The sequence (x t ) now denotes a sequence of random variables, which is constructed recursively with: x 0 = x (with probability 1) and

∀t ∈ N, ∀i ∈ L, x t+1 i = B(φ(x t i-1 , x t i , x t i+1 )), (16)
where B(p) is the Bernoulli random variables, which equals 1 with probability p and 0 with probability 1 -p. Note that strictly speaking, the definition above is more a characterisation than a definition and that a "proper" definition would require the use of tools from measure theory (see e.g. [?]). We now need to redefine what it means to solve the problem perfectly. The blinking condition is easily translated to φ(0, 0, 0) = 1 and φ(1, 1, 1) = 0. For a rule which verifies the blinking condition and a given configuration x ∈ E n , we define the synchronisation time T (x) as the random variable which corresponds to the number of steps needed to attain one of the two homogeneous configurations: T (x) = min{t, x t ∈ {0, 1}}. The average synchronisation time of x is the expectancy of T (x), denoted by E{T (x)}. For a given size n, we define the worst expected synchronisation time (WEST) of x and the expected average synchronisation time (EAST) of x as:

WEST(n) = max x∈En E{T (x)}, ( 17
) EAST(n) = 1 2 n x∈En E{T (x)}. ( 18
)
We say that f synchronises the size n if WEST(n) is finite. Clearly, this is equivalent as having a finite EAST(n). By extension, f is a solution to the global synchronisation problem if f synchronises all sizes n ∈ N.

Let us now examine how to build a solution. For a function φ, we introduce the variables: p 0 = φ(0, 0, 0), p 1 = φ(0, 0, 1), . . . , p 7 = φ(1, 1, 1). Proposition 6. Let φ be a probabilistic ECA such that p 0 = 1, p 7 = 0 and ∀i ∈ {1, . . . , 6}, p i ∈]0, 1[, then φ is a solution to the global synchronisation problem.

Proof (Sketch). To show that T (x) is finite for every x ∈ E, it is sufficient to show that there is a non-zero probability to reach 0 from x in a finite number of steps.

Let Z(x) be the function such that Z(x) i = 1 if (x i-1 , x i , x i+1 ) = (0, 0, 0), 0 otherwise . The effect of Z is to change a cell to a 0 whenever there is a non-zero probability that this cells updates to 0. Let us denote by m the greatest time needed to shrink the regions of ones; clearly, m = n/2 . It can then easily be checked that:

(a) ∀x = 1, Z m (x) ∈ {0, 1}, and (b) there is a non-zero probability to have ∀t ∈ {1, . . . , m}, x t = Z t (x). In other words, there is a non-zero probability that x is absorbed in the 0-1 cycle in at most m time steps. Viewing the system as a Markov chain, this corresponds to the fact that the only two recurrent states are 0 and 1 and all the other states are transient5 .

Note that the proposition above only gives sufficient conditions to solve the problem; moreover, it does not give any idea on the time needed to converge to a fixed point. This means in particular that the function WEST(n) may scale exponentially with n, which is not what is expected for an efficient solution to the problem.

Proposition 7. Consider the probabilistic ECA φ defined with the relationship: φ(x, y, z) = (y + z)/2; this rule is such that: W EST (n) = Θ(n 2 ).

Proof (Sketch). By construction, the rule is the composition of two (commutative) operations: (a) an α-asynchronous left shift with α = 1/2 and (b) a global inversion (x → x). In the α-asynchronous updating each cell independently applies the local rule with probability α and keeps its state with probability 1 -α. This rule was analysed in a previous work [?] and it was shown that its WECT scales quadratically with n.

It can be verified that the global inversion preserves the dynamics of the asynchronous shift. Indeed, the effect of φ is simply to shift the interfaces between regions of 0s and regions of 1s. To formalise this, one could use a coupling between the original process and the asynchronous shift.

The same construction can be applied to the rules proposed by Fukś and Schüle to solve the density classification problem with stochastic rules. These rules were also analysed and were shown to have a quadratic scaling of their convergence time [?].

It is however interesting to note that the Traffic-Majority rule [?], whose convergence is conjectured to be linear with the ring size, does not obey the invariance by global inversion. Intuitively, the reason of this non-symmetry is that the Traffic rule (i.e., ECA 184) treats the 0s and 1s differently and that a global inversion also reverses the direction in which each state is "translated". It is an open question to find a rule with a synchronisation time that has a linear scaling.

Questions

We studied the global synchronisation problem and listed a few simple properties the potential solutions. By formulating the problem as a SAT problem, we could perform a first systematic exploration of the existence of perfect solutions and give a different point of view than the techniques that have been used so far (see e.g. Ref. [?]).

We also noted the existence of a huge gap between deterministic and stochastic systems: the use of randomness allows one to easily obtain a "perfect" solution in the sense that any configuration will be almost surely synchronised in finite time. The precise estimation of the average time of convergence is a delicate operation in all generality but insights could be given for some precise rules.

We end this note with a list of questions and indications:

Question 1. Does there a exist a deterministic rule which synchronises all sizes?

As many researchers do, we believe that the answer to this question is negative.

Question 2. If the answer is no, given a neighbourhood of size k, what is the maximum ring size that can be synchronised? Is this function computable?

We have absolutely no hint on how to answer this question.

Question 3. Given a neighbourhood of size k and maximum synchronisation time τ , what is a good algorithm to find all the rules with a height less or equal to τ ? What is the complexity of this problem?

The work presented here with the use of SAT solvers can largely be improved. This is only a first attempt to use such techniques. Research could be continued by looking for other symmetries of the problem or other ways to add the constraints expressed in Prop. ??. Naturally, other paths have also to be searched.

Question 4. Is there a stochastic solution the global synchronisation problem whose WEST scales linearly with the ring size?

At the moment, we do not see how such a rule could be constructed.

Question 5. To which general context is it worth to generalise the global synchronisation problem?

One may think of higher dimensions, more states, non-homogeneous rules6 , Boolean networks, etc.

Fig. 1 .

 1 Fig. 1. Four space-time diagrams for rule 1078270911 with k = 5 (see page ??). Time goes from bottom to top, white and blue squares represent cells in state 0 and 1, respectively. (a), (b): synchronised initial conditions with n = 11, (c) and (d): initial conditions with n = 12, the last one is not synchronised.

Fig. 2 .

 2 Fig. 2. Transition graphs of ECA 61 for n = 3 and n = 4. An oriented link between a configuration x and y represents the relationship y = F (x).

Fig. 3 .

 3 Fig. 3. Neighbourhood of size k = 4: transition graphs of rule 5419 for n = 4 and n = 5. For the sake of readability configurations have been represented by a number which corresponds to the decimal conversion of their bits.

A transition is a tuple of size k which is given as an input to the local rule.

It corresponds to sequence A000031 in the Online Encyclopedia of integer sequences. One reads: "In music, |χ(n)| is the number of distinct classes of scales and chords in an n-note equal-tempered tuning system", see: https://oeis.org/A000031

see: http://minisat.se/

see http://fiatlux.loria.fr

The definitions of recurrent and transient can be found in the introductions to Markov chain theory. Informally, a recurrent state corresponds to a state who is returned to an infinite number of times with probability 1 and a transient state is a state which is not recurrent: it will then be "leaved" definitively with probability 1.

Readers may test ECA 3 with two boundary cells that have a state 0.

Acknowledgements

The author is grateful to Irène Marcovici, Jean Mairesse and Sukanta Das for stimulating discussions on the topic and to the anonymous reviewers for their precious remarks and suggestions.


Paolino Tona

Hoai-Nam Nguyen

Guillaume Sabiron

Yann Creff

An Efficiency-Aware Model Predictive Control Strategy for a Heaving Buoy Wave Energy Converter

Keywords: Wave energy converters, model predictive control, short-term wave prediction, conversion efficiency

In this paper we consider wave energy converters (WECs) of the heaving buoy type. We assume that the available power take-off (PTO) mechanism is able to transform mechanical to electrical power (generator mode) and electrical to mechanical power (motor mode). In this context, a control system can be designed so as to take off as much energy as possible from the waves to the network, for a broad range of sea states. Among the implementable strategies surveyed in the literature, model predictive control (MPC) has the potential to lead to almost optimal performance, because of its ability to deal with many design objectives and performance criteria and to handle constraints on states and control inputs. However, a critical factor is often overlooked in the design and validation of MPC (as well as of most other WEC control strategies): the imperfect power conversion in the PTO. The MPC strategy presented here does take it into account, in the form of a conversion efficiency function in the criterion to be minimised. First, an optimal control problem is formulated and solved offline, which maximises the average net power output of the wave energy converter. Its performance is tested on five scenarios corresponding to different sea states and compared to that of a conventional control, based on proportional-integral (PI) velocity feedback. Then a nonlinear MPC strategy is introduced and its performance is tested on the same scenarios. It is shown that it yields roughly the same results as the optimal control, and can help harvesting up to 50% more energy than the reference PI control.

I. INTRODUCTION

Among the various Wave Energy Converter (WEC) concepts, many are based on buoyant bodies placed in the sea (see e.g. [START_REF] Drew | A review of wave energy converter technology[END_REF] for a survey). The oscillations caused by the waves passing the body can be converted into electricity by a mechanism called Power Take-Off (PTO). Some PTOs allow four-quadrant operation, i.e. they are able not only to transform mechanical to electrical power (generator mode) but also electrical to mechanical power (motor mode). This feature provides for a great deal of flexibility to control the body motion and bring it into resonance with the incident wave, in order to maximise the net harvested energy.

Among the implementable control strategies that have been proposed in this framework (see [START_REF] Hals | A comparison of selected strategies for adaptive control of wave energy converters[END_REF]), such as proportional velocity feedback control, approximated complex-conjugated control or model predictive control, the latter has the potential to lead to almost optimal performance, because of its ability to deal with many design objectives and performance criteria and to handle constraints on states and control inputs.

Examples of model predictive control (MPC) strategies, based on linear or nonlinear models, convex or non-convex criteria, can be found in [START_REF] Cretel | An application of model predictive control to a wave energy point absorber[END_REF]- [START_REF] Li | Model predictive control of sea wave energy converters part i: A convex approach for the case of a single device[END_REF]. Most strategies are formulated so as to maximise the net produced energy over a given horizon, which corresponds (for large enough horizons) to the ultimate goal any WEC control system should strive for. However, to the best of our knowledge, no MPC strategy has been proposed so far which takes directly into account the imperfect power conversion in the PTO by introducing a (nonlinear) conversion efficiency, function of the instantaneous power, in the optimisation criterion. Notice that, even though this practice is seldom followed even in the most recent literature, conversion efficiency should be always taken into account in the simulation model to obtain realistic results and a fair comparison among controllers, as shown in [START_REF] Sanchez | Control performance assessment and design of optimal control to harvest ocean energy[END_REF] for PI velocity feedback control and in [START_REF] Saupe | Latching control strategies for a heaving buoy wave energy generator in a random sea[END_REF] for latching control.

Section II describes the system studied, and introduces the mathematical model and control parameterisation used in the following. In Section III, an optimal control problem which takes into account PTO conversion efficiency is formulated and solved offline, assuming complete knowledge of the sea state over the whole optimisation interval. Optimal control results on five different sea states are compared with those obtained with a standard proportional-integral (PI) velocity feedback control strategy, whose parameters are optimised off-line for each sea state. In order to eliminate the need for a far-ahead prediction of the wave force, and overcome the computational burden of the optimal control, a movinghorizon model predictive control (MPC) strategy is introduced in Section IV. Using once again the PI control as a benchmark, it is shown that the MPC strategy can extract up to 50% more power from the waves. The last section gathers the main open issues and future topics worth investigating.

II. SYSTEM DESCRIPTION A. Wave Energy Converter Modeling

A wave energy converter of the heaving buoy type, outlined in Fig. 1, is considered. Its mechanics is modelled as a (second order) rotational mass spring damper system. The spring represents the hydrostatic restoring force. The input is the sum of all acting forces (actuator force, radiation force and wave excitation force) and the output is the velocity. Wave interaction is considered by modelling a radiation resistance force f r which is calculated from the velocity using order five dynamics. Actuator dynamics is modelled using a second order low-pass filter giving the relationship between the commanded force u c of the PTO system and the actual force u. The maximum and minimum amplitudes of the force u are limited by physical limitations at ±u max .

A block diagram of the complete model is depicted in Fig. 2. The following system representations are introduced: the actuator dynamics

M c a ẋa = A c a x a + B c a u c (1) u = C c a x a . (2)
and the system dynamics

M c s ẋs = A c s x s + B c s (w -u) (3) v = C c s x s . ( 4
)
which comprises the model of the mechanics and the radiation force. w denotes the wave force which is assumed to be a measurable exogenous signal. In the following, five different realistic sea states with spectra shown (after smoothing) in Fig. 3, are considered for the generation of w.

B. Control Parameterisation and System Discretisation

In the optimal control framework, one often works with parameterised controls. In this work, the standard approach of a piecewise constant control signal is chosen. Accordingly,

x a,i+1 = A a x a,i + B a u c,i (5)
u i = C a x a,i + D a u c,i . (6)
For the actuator dynamics, a zero-order hold discretisation is chosen in accordance to the piecewise constant parametrisation of the commanded control. Discretised with a zero order hold, the actuator dynamics have no direct feed-through (i. e. D a = 0). But other choices are possible, so a generic D a is used in the following.

The system dynamics are discretised using Tustin's rule. Both the wave force and the PTO force (the filtered commanded force) are continuous signals s. t. Tustin's rule yields a much better approximation for a given sample time than the zero-order hold. A discretisation using Tustin's rule results in a system representation which is not strictly proper (i. e. D s = 0, where D s is the feed-through matrix for M sys ).

A sample time of 0.1 sec is selected for the discretisation of both models.

For a certain horizon of n samples, starting to count at sample 0 (i. e. now or t = 0), the sampled vectors of commanded/actual control, velocity and wave force are given as

u T c = u c,0 , u c,1 , . . . , u c,n-1 , (7)
u T = u 0 , u 1 , . . . , u n-1 , (8)
v T = v 0 , v 1 , . . . , v n-1 , (9)
w T = w 0 , w 1 , . . . , w n-1 .

Using ( 5)-( 6), the expressions ( 8)- [START_REF] Falnes | A review of wave-energy extraction[END_REF], and the initial values of the actuator and system dynamics, relationships between u and u c and between v and (wu) can be established. The actuator dynamics can be written as

u = Bc→u        D a 0 0 • • • 0 C a B a D a 0 • • • 0 C a A a B a C a B a D a • • • 0 . . . . . . . . . . . . 0 C a A n-1 a B a C a A n-2 a B a C a A n-3 a B a • • • D a        u c + +        C a C a A a C a A 2 a . . . C a A n a       
x a,0 U0 [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] In analogy to B c→u and U 0 , the matrices B u→v and V 0,x are defined for the system dynamics. Using these expressions, the system dynamics can be written as

v = B u→v (w -u) + V 0,x (12) = B u→v (w -(B c→u u c + U 0 )) + V 0,x . (13)
With ( 11) and ( 13), the PTO force and the velocity are expressed in terms of the commanded force, the incident wave force and the initial states.

III. OFF-LINE OPTIMAL CONTROL

A. Control Objective

The control objective is to take off the maximum amount of energy/mean power from the considered machine by optimising the PTO control strategy. Since this mean power is taken from the system, it is given a negative sign and thus, the control objective corresponds to the minimisation of the mean power. The mean power P c m is easily calculated as the normed integral over the PTO force multiplied with the device velocity and the actuator efficiency η:

P c m = - 1 T T t=0 ηuvdt. (14)
In [START_REF] Falnes | A review of wave-energy extraction[END_REF], the mean power output is characterised by assuming an almost perfect PTO with losses that are only proportional to the device velocity and η = 1. In this case, the average mean power that can be taken off the machine can be calculated based on the mean losses, the integral [START_REF] Shook | Identification for long-range predictive control[END_REF] does not need to be evaluated. This characterisation allows for an elegant expression of the mean power output as a convex quadratic function of the control force. It is straightforward to use such an objective function in a standard MPC formulation which can be solved using readily available software tools as done in e. g. [START_REF] Hals | Constrained optimal control of a heaving buoy wave-energy converter[END_REF].

Unfortunately, this formulation of the objective function is unrealistic because the PTO is never one-hundred-percent efficient. This makes power that is taken from the grid more expensive and reduces the value of the generated power. This phenomenon has a great impact on the optimal control and thus cannot be neglected. The phenomenon can be modelled by choosing an efficiency which is a function of the ideal instantaneous power uv:

η(uv) = η 0 if uv ≥ 0 1 η0 if uv < 0 , 0 ≤ η 0 ≤ 1 (15)
The coefficient η 0 corresponds to the product of the average efficiencies of each PTO stage (actuator, hydraulic transmission if any, generator/motor, inverters). It depends on the PTO and may even differ between the two cases uv ≥ 0 and uv < 0. Considering ( 15), the losses depend on the instantaneous power and the elegant formulation from [START_REF] Falnes | A review of wave-energy extraction[END_REF] is no longer applicable s. t. it is mandatory to optimise ( 14) directly.

In order to avoid issues in gradient based optimisation, a smoothed approximation of the function [START_REF] Fusco | Short-term wave forecasting for real-time control of wave energy converters[END_REF], continuous in uv = 0, can be used.

B. Formulation of the Optimal Control Problem

In accordance with the previous section, the following optimal control problem has to be solved to harvest the maximum amount of energy from the considered machine:

min uc P c m s. t. ( 16
) (1) -(4) u min ≤ u ≤ u max (17) x min ≤ x ≤ x max (18)
The mean power P c m is a functional of the unknown control signal u c which can be considered to be an arbitrary function of time. To facilitate the numerical solution of the optimisation problem, the control signal u c is restricted to being piecewise constant (on an equidistant grid of time) as already mentioned in Section II-B. Instead of the "function" u c , the sampled vector of the commanded force u c is the new vector of decision variables.

C. Reformulation of the Objective Function

Introducing η(u, v) = diag (η(u 0 v 0 ), η(u 1 v 1 ), . . . , η(u n-1 v n-1 )) , (19)
the integral in the formulation of the instantaneous power P c m (15) is reformulated into an Euler integration

P m = -1/T u T η(u, v)v. (20)
The optimal control problem can only be solved for a given sea state w (with its sampled counterpart w). It is pointed out explicitly, that the sea state must be known for the complete interval over which the optimal control is to be calculated. The formulation of the velocity vector ( 13) is rewritten accordingly to highlight the interdependency of velocity and commanded force.

v = B u→v (w -(B c→u u c + U 0 )) + V 0,x (21) = -B u→v B c→u Bc→v u c + B u→v w -B u→v U 0 + V 0,x V 0 (22) = B c→v u c + V 0 . ( 23
)
Inserting the discrete-time formulation of the dynamics, i. e. ( 11) and ( 23), into (20), yields

P m = -1/T (B c→u u c + U 0 ) T η(u c ) (B c→v u c + V 0 ) (24) = -1/T u T c B T c→u η(u c )B c→v u c + -1/T U T 0 η(u c )B c→v + V T 0 η(u c )B c→u f (uc) T u c + U T 0 η(u c ) V 0 c0(uc) (25)
In ( 25), the function η(u, v) has been replaced with the function η(u c ). The two functions are related as follows

η(u c ) = η(B c→u u c + U 0 , B c→v u c + V 0 ). ( 26
)
If the efficiency of the PTO η(u c ) is neglected, i. e. if η(u c ) = 1, the formulation (25) corresponds to a quadratic cost function. In order to bring it to the standard form

1/2u T c H u c + f T u c + c 0 , H symmetric, ( 27
)
the first term in (25) has to be cast into an expression with a symmetric inner matrix H. Keeping in mind that P m is a scalar and using the following property of any scalar s

s = (s T + s)/2, (28)
the quadratic term is rewritten as

1/2u T c H(uc) 1/T -B T c→v η(u c ) T B c→u -B T c→u η(u c )B c→v u c (29)
In the simple case η(u c ) = 1, the convexity of the objective function can be checked by looking at the definiteness of the inner matrix H(u c ) in (29). This matrix only depends on the actuator and system dynamics. For the problem data of the considered machine, the objective function is non-convex, because H(u c ) is not positive semi-definite.

For an efficient solution of the optimal control problem, the gradient of the objective function

P m = 1/2u T c H(u c )u c + f (u c ) T u c + c 0 (u c ) (30)
is easily calculated as

∂P m ∂ u c = u T c H(u c ) + 1/2u T c ∂H(u c ) ∂ u c u c + f (u c ) T + + ∂f (u c ) T ∂ u c u c + ∂c 0 (u c ) T ∂ u c . ( 31
)
If a non-smooth function for the PTO efficiency is used, the partial derivative in (31) is zero except for the points for which it is not defined. If a smooth approximation of the efficiency η a (uv) is used instead, the partial derivative can be calculated over the whole range of the function.

D. Reformulation of the Constraints

The constraints on the control

u min ≤ u ≤ u max (32)
which must hold for all samples of u, i. e. for all entries of u, can be equivalently formulated introducing constraints on the decision variable u c , through [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF]. With 1 n×1 denoting an n × 1 matrix of ones, the resulting set of constraints is

1 n×1 u min -U 0 ≤ B c→u u c ≤ 1 n×1 u max -U 0 . ( 33
)
The state constraints for the mechanical subsystem can be treated analogously.

E. Smoothing of the Control Signal

In order to avoid unnecessary oscillations in the commanded PTO force u c , a penalty term for variations of u c can be added to the control objective. Such a penalty term P var is easily formulated using backward differences. The vector of command variations is

∆u c = u c,0 u c,1 -u c,0 u c,2 -u c,1 . . . u c,n-1 -u c,n-1 T (34) =         1 0 0 . . . 0 0 -1 1 0 . . . 0 0 0 -1 1 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 0 . . . -1 1         M∆ u c . ( 35
)
Introducing the weight w var which allows to tune the tradeoff between desired performance and smoothness of the command, the penalty term P var is chosen as

P var = w var u max n ∆u c T w var u max n ∆u c (36) = w var u max n u T c M T ∆ M ∆ u c (37) = 1/2u T c H var u c (38)
The normalisation with u max n helps to choose the weight w var . In essence, if the command jumps from its maximum to its minimum value from sample to sample, a weight w var = 1 would lead to an overall penalty term of 2 2 . For the results given in this paper, w var is set to 50. The effect of the smoothing is depicted in Fig. 4. Both versions effectively yield the same power output. However, the smoothed control is much less aggressive and thus to be preferred over the non-smoothed. Furthermore, the solution to the optimisation problem with smoothing converges faster.

F. Scaling of the Decision Variables and the Objective Function

In order to avoid numerical issues, some scalings are introduced for the decision variables u c and the objective function. On the one hand, the decision variables are scaled to the bounds of the control signal (u max ). On the other hand, the objective function average mean power is scaled to the unit kW which corresponds to an overall scaling factor of 1/1000.

G. Formulation of the Optimisation Problem

Using the expressions from the previous sections, the overall optimisation problem, which is referred to as the optimal control problem, is

min uc u T c (H(u c ) + H var ) u c + f (u c ) T u c + c 0 (u c ) s. t. (39) 1 n×1 u min -U 0 < B c→u u c < 1 n×1 u max -U 0 . (40)
The dynamic constraints are implicitly contained in the objective function. State constraints for the mechanical subsystem are not considered to avoid slowing down the computations. If desired, they could be incorporated in analogy to (33).

H. Choice of the Optimisation Algorithm

Since the control objective of (39) is nonlinear and nonconvex with a high number of unknowns, a large-scale nonlinear optimisation algorithm is needed. Because the objective function is continuous and an (approximated) gradient can be determined at low computational cost, a gradient-based optimisation algorithm is favoured. Among the freely available solvers available for Matlab falling into the desired category, the interior point algorithm IPOPT [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF], interfaced via the OPTI toolbox [START_REF] Currie | OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User[END_REF], was chosen, for its performance.

I. Consistency Check of the Optimal Solution

Two different scenarios are considered in order to verify the solution found with IPOPT, i. e. whether a global or a local solution is found, and in order to quantify the possible performance loss. Distribution of average power harvested with optimal controls calculated using random initial values.

1) Random Initial Guesses for the Decision Variables:

In the first scenario, the same optimal control problem is solved for different random initial choices of u c . Results from optimisation runs with 50 random initial values for u c are depicted in Fig. 5. It is clear that the optimisation problem converges toward the same solution for all considered random initial values.

At the same time, the corresponding objective values that are obtained are extremely close as shown in Fig. 6. All results are within a band of about 0.5 % of the same objective value.

2) Test of the Optimum Principle: The second scenario corresponds to a verification of Bellman's optimum principle. It states that an optimal trajectory has the property that at an intermediate point, no matter how it was reached, the rest of the trajectory must coincide with an optimal trajectory as computed from this intermediate point as the initial, see [START_REF] Bellman | On the Theory of Dynamic Programming[END_REF].

In this scenario, two problems are solved. The second problem corresponds to a subproblem of the first, and the optimal control is calculated on an interval starting at an intermediate point from the first problem to the last point of the first problem. Fig. 7 shows that this test is also passed.

J. Application Results and Comparison with Reference Control

The optimal control scenario is solved offline. The execution time of the code is thus unimportant s. t. the problem can be solved over long time horizons. In this case study, a horizon of n = 2000 grid points is chosen which corresponds to 200 sec. This section focuses on the optimal control scenario for the sea state 3, i. e. for waves of medium height. The maximum CPU time is set to 15 min. During this time, on a recent but not particularly fast PC (with an Intel i7 CPU at 2.8GHz), the solver can perform about 50 iterations and the solution converges.

As a reference for the performance of the optimal control a standard PI velocity feedback control (also referred to as BK-control in [START_REF] Sanchez | Control performance assessment and design of optimal control to harvest ocean energy[END_REF]) is used:

u c = -k p t t0 v dτ -k v v. (41)
The proportional gains k p and k v are chosen using a gridding approach: for each sea state and time interval considered, the closed loop is simulated for a grid of gains and the combination leading to the best average power output is selected. If well tuned, the PI control law brings the device into resonance with the incident wave.

The optimal control problem is solved for different values of the efficiency η 0 ∈ {1, 0.9, 0.8, 0.7} in order to study its impact on the results. For the nominal case η 0 = 0.7, the obtained optimal trajectory is compared with the results from the PI velocity control, see Fig. 8. In the considered time interval, the optimal control results in higher amplitudes of displacement and velocity.

All results are also gathered in Table I. The average mean power that can be harvested clearly drops significantly if the PTO efficiency is reduced. At the same time, the optimal control avoids to take power from the grid as seen clearly in Fig. 9 (uv is seldom negative). The ratio between the maximum and minimum instantaneous uv is much lower for the high efficiencies.

IV. MOVING HORIZON MODEL PREDICTIVE CONTROL

In order to overcome the computational burden of the optimal control and to eliminate the need for a far-ahead prediction of the wave force, a moving-horizon MPC can be used. The basic idea of moving-horizon MPC is the following: Instantaneous power with optimal controls for different PTO efficiencies. 1) At current sample i measure (or estimate) the system state and use it as an initial value to calculate an optimal control over a limited horizon with n samples starting at sample i. This yields a time series of optimal control values u c,i of length n. 2) Apply the first value of u c,i (corresponding to the control for the current sample) as a command for the system and hold it for the duration of the sample period. 3) At the next sample i + 1, measure the system state and use it again as an initial value to calculate an optimal control over a horizon with n samples starting at i + 1. Initialise the vector of decision variables using the results from the previous step. This yields u c,i+1 This principle is depicted in Fig. 10.

The main advantages of the moving-horizon MPC scheme is that it is actually implementable in real word applications as an online control law if a short-term wave prediction is available.

A. Penalty for the First Control Value

The transition between the optimal control problem over a single horizon and the moving-horizon problem mainly comes down to the problem of reinitialising the optimisation problem at each time step. To do this correctly, the objective has to be extended in order to penalise the difference between the control for the current step and the previous step u c,pre in analogy to (38)

w var u max n (u c,0 -u c,pre ) (u c,0 -u c,pre ) (42) = w var u max n u T c 1 0 . . . 0 T 1 0 . . . 0 u c -2u c,pre 1 0 . . . 0 u c + u 2 c,pre . (43)
The latter term must be added to the objective function of the optimisation problem (39) which is solved in each step.

B. Compliance with Real-Time Computation Constraints

The prototype control, which is presented in this paper, has not been coded with respect to real-time computation constraints. The control, i. e. the optimisation problem, is solved in a for-loop over the desired time span. In each time step, the optimisation algorithm is terminated after 10 iterations1 and the first entry of the corresponding u c,i is applied in a continuous time simulation of the model in Simulink. The final values of the simulation are recovered (state feedback) and used to re-initialise the optimisation -all according to the general MPC procedure described in Section IV. Despite this simple implementation, the computational time, including the simulation, is only 2.5 times slower than real-time (on the same not particularly fast PC with a 2.8 GHz processor). Using the profiler option of Matlab, it can be found that only about 45 % of this time is spent solving the MPC problem. Since these results have been obtained without any code optimisation, there is some room to develop a real-time capable version of the control law.

C. Short-Term Wave Prediction

In order to predict the wave force with the methods presented in this section, it must be available as a measurable external signal of the system as assumed in Section II-A (which is quite a strong assumption).

The approach used to predict the waves corresponds to the online identification of an autoregressive (AR) model of the order n AR :

y i = n AR k=1 a k y i-k . (44)
The parameter vector

a = a 1 a 2 . . . a n AR T (45)
characterises the AR model. An equivalent representation of the AR model (44) is found with

        y i y i-1
. . .

y i-n AR +1         =            a 1 a 2 . . . a n AR -1 a n AR 1 0 . . . 0 0 0 1 . . . 0 0 0 0 . . . 0 0 0 0 . . . 1 0                    y i-1 y i-2 . . . y i-n AR         , (46)
x i = A ar (a)x i-1 , (47)
y i = 1 0 . . . 0 Car A ar (a)x i-1 . (48)
Once a suitable model of the type ( 46) is identified for the wave signal, it may be extrapolated as

y i+l-1|i-1 (a) = 1 0 . . . 0 Car A ar (a) l x i-1 . (49)
The notation y i+l-1|i-1 (a) specifies that the l-step-ahead prediction of y is calculated based on measurements up to the step i -1 and that this prediction depends on the choice of a.

To cover different sea states, the parameter vector a can be optimised online in an iterative way or in a batch job based on recently recorded data. The AR model identification is performed with the objective of minimising a multi-stepprediction over multiple (past) horizons as proposed in [START_REF] Saupe | Latching control strategies for a heaving buoy wave energy generator in a random sea[END_REF], [START_REF] Shook | Identification for long-range predictive control[END_REF], [START_REF] Fusco | Short-term wave forecasting for real-time control of wave energy converters[END_REF]. With a finite set of samples T up to the current one, such a criterion can be formulated as

J lrp (a) = k∈T l j=1 y k -y k|k-j (a) 2 . ( 50
)
Because of (49), this criterion is a nonlinear least squares problem. No constraints on the decision variables a are enforced s. t. the optimisation problem simply reads

min a J lrp (a). (51)
It is solved using the adaptive nonlinear least-squares algorithm (NL2SOL) interfaced by the OPTI toolbox. Several options are important to obtain good results for the prediction:

• sample time of the AR model ∆t AR ,

• pre-filtering of the wave signal,

• order of the AR model n AR ,

• wave prediction horizon (in samples) n w,pred . 1) Sample Time of the AR Model: The sample time of the AR model has a huge impact on the prediction quality due to its filtering properties and its impact on the effective prediction horizon T w,pred = ∆t AR n w,pred . It must be chosen according to the dominant wave frequency 2 f w,dom .

2) Pre-filtering of the Wave Signal: A pre-filtering of the wave signal can greatly increase the prediction performance, see e. g. [START_REF] Fusco | Short-term wave forecasting for real-time control of wave energy converters[END_REF]. The use of non-causal filtering techniques, to be applied off-line, has been suggested in order to avoid a phase lag of the estimated signal, see [START_REF] Fischer | Online-algortihm using adaptive filters for short-term wave prediction and ist implementation[END_REF]. However, since the issues concerning the trade off between non-causality (which can possibly be compensated by a longer prediction horizon) and phase lag have not been solved yet, no pre-filtering of the wave signal is applied in this work.

3) Order of the AR model: The order of the AR model corresponds to twice the number of harmonics that are modelled. There is a trade-off between model complexity (which increases the difficulty to solve the identification problem) and prediction performance. In this work, n AR is set to 32. 2 The frequency with the highest power. 4) Prediction Horizon: Like the order of the AR model, this parameter represents a trade-off between the complexity of the identification problem and the prediction performance. A relatively large prediction horizon of n w,pred = 50 is chosen. Of that long horizon, only a few first steps approximate the future wave well. The prediction horizon considered in the identification is thus much higher than the prediction horizon needed for the MPC.

Once an AR model is identified, it remains valid for several minutes. In the current implementation, the AR model is updated every 300 sec. Real-time constraints are thus not an issue. The obtained prediction performance for wave 3 is depicted in Fig. 12. It is obvious that it is far from perfect for more than 1 sec into the future. In the final section of this chapter, the performance loss of the MPC control due to an imperfect wave prediction is quantified.

D. State Observation

For the relatively simple linear model considered in this case study, the impact of an imperfect state feedback seems to be negligible. If only the measurable position and velocity are fed back instead of the full state, the performance loss vanishes. The performance loss is established in the subsequent section.

E. Application Results and Comparison with PI Control

The moving-horizon MPC control is applied to the considered machine over the complete time intervals of the available wave data. The nominal case with η 0 = 0.7 is considered. The results are then compared to the results achieved with PI velocity feedback controllers especially optimised for the respective waves. It should be noted that, in a realistic setting, a gain-scheduling policy has to be implemented in order to adjust the PI gains to the current sea state, which has to be estimated at regular intervals. In contrast to that, the MPC strategy remains the same for all five sea states. Fig. 13, Fig. 14 and Fig. 15 depict the results, i. e. float position and velocity, PTO force and power output of the MPC and the PI control laws, for wave states 1, 3 and 5 respectively. The time intervals correspond to the intervals during which the maximum deflections occur.

A summary of the results is found in Table II. The gain in harvested energy lies between 16 % for high waves and 50 % for low waves. For the representative time interval [301 sec, 501 sec], Fig. 16 depicts the performance gain (with respect to the PI control) of the optimal control, an MPC with state feedback and perfect wave prediction (MPC 1 ), an MPC with output feedback and perfect wave prediction (MPC 2 ) and an MPC with output feedback and AR identification based wave prediction (MPC 3 ).

The perfect MPC and the optimal control practically yield the same results. Most of the performance is lost due to the imperfect wave prediction based on the AR model identification. Compared to the optimal control or perfect MPC, the performance loss with respect to the PI control lies between 3 % and 7 %. V. CONCLUSIONS A nonlinear MPC strategy for a wave energy converter of the heaving buoy type has been presented, which takes explicitly into account the conversion efficiency of the PTO system to maximise the mean electrical power extracted from the WEC. It has been shown that the MPC strategy yields results very close to the solution of the corresponding optimal control problem solved offline and has the potential to harvest, depending on the sea state, from 15 to 50% more energy than a conventional controller based on PI velocity feedback.

Nevertheless, a few issues need to be solved and/or further investigated.

A. Short-Term Wave Prediction

The performance of the MPC strategy clearly deteriorates with imperfect wave force prediction. Moreover, the prediction approach adopted here requires periodic analysis of past measurements in order to adapt not only the parameters of the AR model to the current sea state, but also its sample time. Alternative approaches exist, such as in [START_REF] Fischer | Online-algortihm using adaptive filters for short-term wave prediction and ist implementation[END_REF], and should be tested in the future to verify if they yield superior performance in combination with the MPC scheme.

B. Availability and Quality of Wave Force Signal

In order to predict future wave force from current and past values, those values must be available. However, wave force is not directly measurable online in a realistic set-up. It can be approximated from wave height measurements, but it must be recalled that a direct measurement of wave height at the center of the float is impossible during normal operation. Some approaches (usually model-based) have already been proposed to estimate it from commonly available sensors, but their performance has yet to be proven for a broad range of operating conditions.

C. Computation Time

Though it seems very likely that, optimising the code, a realtime capable version of the control law should be realisable, the proposed strategy remains intrinsically demanding in terms of computations. To speed up its convergence, it may be advantageous to compute and use an analytic Hessian. Another line of research, which would require some non-obvious theoretical development, consists in trying to convexify the criterion to be minimised.

D. Robustness

The robustness of the design with respect to unmodelled dynamics need to be investigated. It must be established if there are dominant nonlinearities to be considered in the system and how they affect prediction quality and computation time.

VI. ACKOWLEDMENTS

Fig. 1 .

 1 Fig. 1. Simplified model of the considered machine.

Fig. 2 .

 2 Fig. 2. Structure of the considered mathematical model of the considered machine.

Fig. 3 .

 3 Fig. 3. The (smoothed) spectra of the five considered sea states.

Fig. 4 .

 4 Fig. 4. Impact of the smoothing on the control.

control / 10 6 NFig. 5 .

 65 Fig. 5. Optimal control trajectories calculated using random initial values.

Fig. 7 .

 7 Fig. 7. Verification of the optimum principle: Solution on a sub-interval yields the same trajectory as obtained on the complete interval.

Fig. 8 .

 8 Fig. 8. Comparison of optimal control and PI control for the wave 3.

Fig. 9 .

 9 Fig. 9.Instantaneous power with optimal controls for different PTO efficiencies.

Fig. 10 . 1 Fig. 11 . 4 )

 101114 Fig. 10. The principle of sliding horizon MPC.

Fig. 12 .

 12 Fig. 12. Wave prediction quality.

Fig. 13 .

 13 Fig. 13. Comparison of MPC and PI control for the wave 1.

Fig. 14 .

 14 Fig. 14. Comparison of MPC and PI control for the wave 3.

Fig. 15 .

 15 Fig. 15. Comparison of MPC and PI control for the wave 5.

Fig. 16 .

 16 Fig. 16. Performance gain (with respect to the PI control) of the optimal control, an MPC with state feedback and perfect wave prediction (MPC 1 ), an MPC with output feedback and perfect wave prediction (MPC 2 ) and an MPC with output feedback and AR identification based wave prediction (MPC 3 ). Only the representative time interval [301 sec, 501 sec] is considered.

TABLE I AVERAGE

 I MEAN POWER OUTPUT OF THE OPTIMAL CONTROL WITH DIFFERENT EFFICIENCIES.

			Optimal control		PI
	η 0 / -	1	0.9	0.8	0.7	0.7
	Pm / kW	-27.86	-22.26	-19.06	-16.32	-12.09

TABLE II AVERAGE

 II MEAN POWER OUTPUT OF THE MPC CONTROL VS. THE PI CONTROL OVER THE WHOLE TIME INTERVALS OF THE WAVE DATA.

		Pm PI	Pm MPC gain Pm
		/ kW	/ kW
	wave 1	-2.83	-4.26	50.43 %
	wave 2	-8.64	-11.89	37.54 %
	wave 3	-16.42	-20.76	26.44 %
	wave 4	-24.44	-29.26	19.75 %
	wave 5	-32.14	-37.14	15.55 %

A termination of the algorithm after only 10 steps is possible because an interior points algorithm is applied which yields a series of feasible points.

The authors wish to thank Florian Saupe, now consultant at McKinsey & Co., as this paper draws extensively on the research work he did while at IFP Energies nouvelles.


F Saupe

J.-C Gilloteaux

P Bozonnet

Y Creff

P Tona

Latching Control Strategies for a Heaving Buoy Wave Energy Generator in a Random Sea

This paper addresses the latching control of wave energy converters. The principle of this control approach is to bring the wave energy generator into resonance with the incident wave using a clamping mechanism. Maximum energy extraction is the control objective. The main challenge in any latching control scheme is to calculate the precise time when to release the device after it has been locked at zero velocity. At the example of a generic heaving buoy device and using real wave data, three latching strategies are compared to a PI velocity control. The simplest considered latching strategy releases the device as soon as the wave force reaches a certain threshold. The other strategies use a short-term wave prediction in order to calculate the latching timing. Imperfect wave predictions based on AR models and imperfect mechanical to electrical and electrical to mechanical power conversions are taken into account. While the imperfect wave predictions impact the achievable performance with the predictive latching strategies, the imperfect power conversions have a high impact on the PI velocity control due to its reactive power flow.

INTRODUCTION

The control of wave energy converters via a discontinuous control law, commonly referred to as latching, is addressed in this paper. The principle of latching control is to bring the wave energy generator into resonance with the incident wave by using a clamping mechanism. It is ensured that the latching control needs to draw no significant amounts of power from the net by only activating the clamping mechanism if the velocity of the device is zero. The main aspect of any latching control scheme is the calculation of the precise time to unlatch the device. The optimal timing depends on the incident wave. In a regular sea with monochromatic waves, the problem of determining the optimal timing is easily solved, see e. g. Babarit and Clément (2006). However, in a more realistic irregular sea, the same problem is much harder to solve and has been an active field of research for decades, see Budal et al. (1982); Babarit and Clément (2006); Babarit et al. (2004);Falcão (2008).

Several approaches to latching control in an irregular sea are considered. The first is a simple but effective threshold unlatching which releases the device if the incident wave force reaches a certain threshold, Falcão (2008). The second latching scheme uses a prediction of the incident wave force to calculate a latching timing which is optimal on the considered prediction horizon. The calculations are performed using the full WEC model, Babarit et al. (2003Babarit et al. ( , 2004)). In a novel approach, proposed in this paper, the time to unlatch the device is also determined based on a short-term prediction of the incident wave. It is simpler than the approach from Babarit et al. (2003,2004) as the available horizon of the wave force prediction is only searched for the time interval with the highest mean up/down force. An eigenvalue analysis of the wave energy generator is sufficient to determine the length of the interval.

Proportional integral (PI) velocity feedback is considered as a reference control. This approach is simple to apply effectively in irregular waves and thus popular in the industry, see e. g. Hansen and Kramer (2011). The downside of a PI velocity control is its reactive power flow, i. e. the fact that power is also taken from the net. The different control strategies are presented in Section 3.

In the literature on short-term wave prediction, algorithms based on auto regressive (AR) models have been established as the state of the art, see e. g. Fusco and Ringwood (2010b,a); Schoen et al. (2011). Accordingly, the robustness of the latching schemes which use a prediction of the incident wave force is analysed by replacing perfect predictions with more realistic predictions based on AR models. Section 4 gives a description of the applied wave prediction algorithm.

The performances of the considered control schemes are benchmarked based on a scenario using wave data recorded at Galway Bay by the Irish Marine Institute. Section 2 presents a generic model of a heaving buoy wave energy generator which is used for the simulations. Imperfect mechanical to electrical and electrical to mechanical power conversion efficiencies are taken into account. Since this effect makes the power which is taken from the net more expensive compared to the produced power, these efficiencies have a great impact on the achievable performance with the PI velocity control and accordingly, the advantages of latching show off nicely. The results of the simulation study are presented in Section 5.

MODELLING APPROACH

This section introduces the mathematical model of a generic heaving buoy point absorber wave energy converter and gives a short characterisation of the considered sea state.

Model of a Generic Wave Energy Converter

A wave energy converter of the heaving buoy type is considered. The incident wave excites oscillations of the buoy. The installed power take off device (PTO) is used to convert the kinetic energy of the device into useful electric energy. The buoy and the wave interact via the hydrodynamic restoring force (modelled as a spring), a simple radiation force (modelled as a damper) and a constant added mass. The system is depicted schematically in Figure 1. This simple modeling approach is chosen over an approach which takes into account dynamic radiation forces because it allows to change the important system characteristics in an uncomplicated way. Being able to take into account systems with different resonance frequencies allows to establish interesting insights into the suitability of the considered control approaches. It is pointed out that all considered control approaches are also applicable (without modification) to the more general case of a system with a dynamic radiation force and a dynamic added mass term.

With the position q of the oscillating body, the eigenfrequency ω 0 , the damping ζ, the incident wave force w and the PTO force u, the equations of motion of the wave energy converter model are q + 2ζω 0 q + ω 2 0 q = w + u.

(1) The relation between the physical parameters mass (including the added mass) m, damping coefficient d, stiffness coefficient k and the parameters parameterising the equations of motion ω 0 , ζ are

ω 0 = k m , ζ = d 2ω 0 m . ( 2
)
Two example systems are considered in this study. Both have a damping of ζ = 0.1 but differ in their eigenfrequencies which lie at ω 0 = 0.22 Hz and ω 0 = 0.33 Hz.

The instantaneous power output1 P i of the device is modelled as a nonlinear function involving the product of the PTO force u and the velocity q and an efficiency whose value (η or 1/η) depends on whether energy is taken from or fed into the net.

P i = -ηu q if u q ≥ 0 -1 η u q if u q < 0 , 0 ≤ η ≤ 1 (3)
In essence, this formulation of the instataneous power considers the important power conversion losses in both directions, from mechanical to electrical power and vice versa. This makes power taken from the net more expensive from a control point of view.

Wave Force

A wave recorded by the Irish Marine Institute (wave rider buoy data from Galway Bay) is used to represent an irregular sea state. The wave is sampled at a frequency of 2.56 Hz. An excerpt of the signal and its spectrum are depicted in Figure 2. The wave force w is considered to be directly proportional to the wave elevation.

CONTROL OF THE WEC

This section states the control objective and presents different control strategies. Among the considered strategies are a PI velocity control, force threshold unlatching and two latching control schemes based on a short term wave prediction.

Control Objective

The control objective is to maximise the energy capture.

It is pointed out that mean power which is taken from the system (fed to the net), has a negative sign and thus, the control objective corresponds to the minimisation of the mean power. The mean power P m is easily calculated as the normed integral over the instantaneous power P i (3)

P m = 1 T T t=0 P i dt. ( 4
)

PI Velocity Control

The considered reference control is a PI feedback of the PTO velocity 2 . u = -k v qk p q.

(5)

Due to the simplicity of the control law, a brute force grid-based search method can be applied to find the gains yielding the optimal mean power output achievable with (5). The PI velocity feedback leads to reactive power flows which are more or less significant depending on the chosen gains. The ratio of peak to average power reaches 31 for the optimal gains in the case of perfect PTO efficiency (η = 1). For PTO efficiencies less than 1, high reactive power flows severly decrease the mean power output s. t. the optimal PI gains depend drastically on η. The main advantage of the PI velocity control is that it does not need measurements of the wave force.

Latching Control

In a latching control strategy, the moving parts of the wave energy converter are successively latched and unlatched. This way, the system is brought in phase (and thus in resonance) with the sea. The approach allows to achieve high velocities which in turn allow to take off high amounts of energy.

If a latching control is applied to the system (1) , the closed loop is a switching system with two phases.

(1) The unlatched system with the equations of motion (1) and the PTO force control law u = -k v,l q. (6) (2) The latched system with q = 0 and q = 0 (no movement).

The system can change from phase 1 to phase 2 only if q = 0. This ensures that no energy is consumed for braking the device. Once at rest, the device can be unlatched at any time.

The biggest benefit of latching control is that it is passive, i. e. that it does not need to draw energy from the net (k v,l is negative s. t. -u q is negative which corresponds to an energy flow to the net). Active controls, i. e. controls that allow for reversible energy flows, need PTOs of higher complexity which are more expensive.

The challenge in latching control consists of finding the precise times at which to unlatch the body -especially in a realistic irregular sea state. But before the strategies to determine the latching timing are discussed, a choice for the PTO force control law (6) is made.

Choice of the PTO Force Control Law

The gain k v,l in the PTO force control law ( 6) is often considered as given (see e. g. Babarit et al. (2004); Babarit and Clément (2006)) or, for simpler schemes, found by gridding over all possible values (see e. g. Falcão ( 2008)). In order to systematically determine a good practise for the choise of k v,l , the optimal latching control problem in monochromatic waves 2 This corresponds to a full state feedback for the simple model considered in this paper

w = A w cos(ωt + ϕ) (7)
is considered in this paper. In an irregular sea with the dominant frequency ω d , the gain which is optimal for a monochromatic wave with the frequency ω d is implemented.

The optimal latching control in monochromatic waves is determined using the knowlege that the system approaches a periodic trajectory on which the system switches between the two phases 'latched' and 'unlatched'. In the highest and lowest point of the periodic trajectory, when the velocity is zero, the system is latched. The phase ϕ of the wave when the system is to be unlatched can be determined based on the following considerations, see e. g. Babarit and Clément (2006); Babarit et al. (2004). The initial state of the system in the moment of unlatching is the position q 0 = ζ 0 and velocity q0 = 0. From the periodic nature of the final trajectory, it can be deduced that after half a period of the controlled motion, the system is latched again in the position q 1 = -ζ 0 when its velocity is q1 = 0. Instead of directly maximising the mean power output, the problem is usually simplified to finding the ϕ which leads to the maximum amplitude ζ 0 . In the following treatise, x = [q q] T .

With the undamped eigenfrequency

µ = ω 2 0 -ζ 2 ω 2 0 (8)
and a, b (which depend on the initial values), the solution of the wave energy generator equations of motion (1) to a monochromatic wave ( 7) is given by

q(t) = (a cos (µt) + b sin (µt)) e -ζω0t + + |H(ω)| A w cos (ϕ + ωt + H (ω)) . (9)
Using this explicit solution of the wave energy converter's equations of motion (9) (without loss of generality, the initial time t 0 is set to 0), the problem which is to be solved in order to find the optimal φ can be formulated as max

ϕ ζ 0 s. t. (10) x(t 0 ) = [ζ 0 0] T , x(t 1 ) = [-ζ 0 0] T and (9) (11)
There are 5 unknowns ϕ, a, b, t 1 and ζ 0 which can be solved for using the 4 boundary conditions (11) and the optimum condition max ζ 0 . Because the equations are highly nonlinear, a semi-numerical solution is called for e. g.

(1) eliminate a and b using x(t 0 ) = [ζ 0 0] T , (2) using a and b from step 1, solve x(0) = -x(t 1 ) for ζ 0 (now a function of t 1 and ϕ), (3) find the maximum ζ 0 over a grid of t 1 and ϕ.

It is not surprising that, for excitation frequencies slower than the undamped eigenfrequency µ of the device, the optimal t 1 equals half a natural period of the undamped eigenfrequency T = 2π/µ. (12) The control approach brings the system in resonance with the wave force. The resulting velocity of the wave energy converter is in phase with the excitation force. Similar calculations are also possible for systems of higher complexity, e. g. for systems with a dynamic radiation force/added mass, see Babarit and Clément (2006). Using the above results, the gain k v,l which leads to optimal energy extraction can be found for each wave frequency by performing a line search consisting of two steps.

(1) The damping of the system is modified to account for the proportional gain

k v,l s. t. d = d + k v,l .
Starting from k v,l = 0.05ζω 0 , k v,l is increased and the optimisation problem ( 11) is solved for each value of k v,l . Based on the obtained latching timings, the mean power output with the individual k v,l over t 1 is calculated. This is repeated until the obtained mean power decreases for the first time in the k th step. (2) A golden section search is started between the (k -2) th and the k th step in order to find the k v,l yielding the highest power output.

The results of this approach for the considered wave energy converter are depicted in Figure 3 for the case ω 0 = 0.22 Hz.

A gain k v,l obtained by the procedure described above leads to good results in irregular waves as shown in the following section. From the characteristic of the optimal k v,l over the frequency, given in Figure 3, it is already clear that k v,l does not change significantly over a wide band of relevant frequencies around the eigenfrequency of the device.

Threshold Unlatching This simple but effective scheme to determine the latching timing in an irregular sea was proposed e. g. in Falcão (2008). The device is unlatched as soon as the incident wave force exceeds a certain threshold. Accordingly, this strategy merely needs a measurement of the incident wave force and no prediction. The choice of the force threshold depends on the sea state. Again, the control law can easily be optimized by using a grid based search (the threshold and the gain k v,l are the only parameters).

At the example of the threshold unlatching control it is shown that the choice of the gain k v,l , which is proposed in Section 3.3.1, leads to good results in irregular waves. Simulations are performed with the optimal k v,l , obtained from the grid based search, and the k v,l which results from the procedure given in Section 3.3.1. The two resulting energy outputs are almost equal: the difference is below 0.5 %.

Model Based Predictive Latching

The model based predictive latching is proposed in Babarit et al. (2004). Each time the system is latched after reaching zero velocity, the prediction horizon is gridded leading to a sequence of possible unlatching times t i . For each t i , the energy generated between t i and the next zero-crossing of the velocity is calculated by simulating the system response.

In the end, the system is unlatched at the t i leading to the highest energy output.

It is pointed out explicitly, that even though the latching timing is selected optimally with respect to the limited prediction horizon, it is not globally optimal. Because of the discrete nature of the control law, a decision at time t a may have undesired effects for a later time t b (which is not part of the prediction horizon). Obvious disadvantages of the scheme are its dependence on the full model of the WEC and the fact that a full state feedback is needed to initialize the simulations used to determine the energy output predictions.

Mean Wave Force Predictive Latching

An alternative approach, newly proposed in this paper, is to unlatch the system during a period when the predicted mean wave force has the correct sign and a high absolute value. This approach to calculate the latching timing with a prediction of the wave force is sketched in Figure 4. From Section 3.3.1, it is known that the unlatched motion should roughly last for about half of the period T (12) of the wave energy converter. For models of higher complexity, this corresponds to half a period of the undamped frequency of the system's dominating mode. The predicted wave force is sampled (sample time T s ) s.t. an integer multiple of the sample time corresponds to T /2 = N T s . In the example of Figure 4, N = 4. Starting from each sample of the prediction horizon, windows of length T /2 = N T s are integrated in order to determine their corresponding mean wave force. The integration is performed e. g. using a simple Euler rule.

I i = N +i-1 k=i T s w pred (kT s ) (13)
Based on the information about the predicted mean wave force in all the windows I i , the timing of the unlatching is determined. The window with the maximum mean force is the window with the maximum integral, denoted I max (with index i max ). Should the device be in its lower latched position, it is a good practise to unlatch it at t = (i max -1)T s (the initial time of the window with the maximum mean up-force). This is the scenario depicted in Figure 4. The marked samples represent the period with the maximum mean up-force, i. e. the time during which the device should be unlatched. Respectively, if the device is in its upper latched position, the window with the minimum integral is chosen to time the unlatching. The necessary length of the prediction horizon mainly depends on the characteristics of the incoming wave. A prediction horizon of 1.3 times the dominating period of the wave already yields good results.

This scheme is similar to other schemes which try to synchronize the maximum wave force with the maximum device velocity as e. g. Eidsmoen (1998). Because of the averageing, the proposed approach is more robust with respect to noisy, imperfect wave force predictions.

SHORT-TERM WAVE PREDICTION

The approach most commonly used in the literature (Fusco and Ringwood (2010b,a); Schoen et al. ( 2011)) is the online identification of an AR model of the order n AR :

y i = n AR k=1 a k y i-k . ( 14
)
The parameter vector

a = [a 1 a 2 . . . a n AR ] T (15) characterises the AR model.
An equivalent representation of the AR model ( 14) is found with

    y i y i-1
. . .

y i-n AR +1     =       a 1 a 2 . . . a n AR -1 a n AR 1 0 . . . 0 0 0 1 . . . 0 0 0 0 . . . 0 0 0 0 . . . 1 0           y i-1 y i-2 . . . y i-n AR     , (16)
x i = A ar (a)x i-1 , ( 17
)
y i = [1 0 . . . 0] Car A ar (a)x i-1 . ( 18
)
Once a suitable model of the type ( 16) is identified for the wave signal, it may be extrapolated as

y i+l-1|i-1 (a) = [1 0 . . . 0] Car A ar (a) l x i-1 . ( 19
)
The notation y i+l-1|i-1 (a) specifies that the l-step-ahead prediction of y is calculated based on measurements up to the step i -1 and that this prediction depends on the choice of a.

The AR model identification, i. e. the identification of the parameter vector a, is performed with the objective to minimise a multi-step-prediction over multiple (past) horizons as proposed in Fusco and Ringwood (2010a). With a finite set of samples T up to the current one, such a criterion can be formulated as

J lrp (a) = k∈T l j=1 y k -y k|k-j (a) 2 . ( 20
)
Because of ( 19), this criterion is a nonlinear least squares problem. No constraints on decision variables a are enforced s. t. the optimisation problem simply reads min a J lrp (a). ( 21)

Once an AR model is identified, it remains valid for several minutes. This allows to perform the identification in a batch job using past data, even in real-time applications.

In this paper, a Levenberg-Marquardt nonlinear least squares algorithm is applied. It is noted that in order to predict the wave force with the methods presented in this section, the wave force must be a measurable external signal of the system.

Several options are important to obtain good results for the prediction:

• sample time of the AR model ∆t AR ,

• pre-filtering of the wave signal,

• order of the AR model n AR ,

• wave prediction horizon (in samples) n w,pred ,

• number of horizons considered in the identification n hor .

These options are discussed in the following paragraphs.

Sample Time of the AR Model

The sample time of the AR model has a huge impact on the prediction quality due to its filtering properties and its impact on the effective prediction horizon T w,pred = ∆t AR n w,pred . In this paper, it is chosen depending on the dominant wave frequency 3 ω d s. t.

∆t AR = 2π 13ω d . (22)
Pre-filtering of the Wave Signal A pre-filtering of the wave signal can greatly increase the prediction performance, see e. g. Fusco and Ringwood (2010a). Non-causal filtering techniques are often used in order to avoid a phase lag of the estimated signal, see Fischer et al. (2012).

The issues concerning the trade off between non-causality (which can possibly be compensated by a longer prediction horizon) and phase lag remain an open problem. In this paper, the filter was chosen following the recommendations of Fusco and Ringwood (2010a).

Order of the AR Model

The order of the AR model corresponds to twice the number of harmonics that are modelled. There is a trade-off between model complexity (which increases the difficulty to solve the identification problem) and prediction performance. In this work, n AR is set to 32.

Prediction Horizon Like the order of the AR model, this parameter represents a trade-off between the complexity of the identification problem and the prediction performance. On the one hand, the higher n w,pred , the 'more nonlinear' the identification problem due to (19). On the other hand, a low n w,pred often results in an unstable AR model. This is why a relatively high prediction horizon of n w,pred = 50 is chosen. Of that long horizon, only a few first steps approximate the future wave well. The prediction horizon considered in the identification is thus much higher than the prediction horizon needed for the latching control strategies.

Number of Horizons Considered

The number of horizons which is considered smooths the results and thus increases their long term validity. The more horizons are considered, the higher the computational effort for the evaluation of the objective function of the identification (20). In practise, every 25 th value over a horizon of the last 300 s is considered in the set T .

In the considered application, the first 500 s of wave data are used for the identification of the AR model which is then used to predict the wave. The resulting prediction quality for a 6 s-ahead prediction after 800 s of the wave signal (training window and prediction window are well separated) is depicted in Figure 5.

The poles of the optimised AR model are predominantly on the unit circle (0-damping line) representing the harmonic nature of the identified signal.

SIMULATION RESULTS

A simulation study is performed using the two configurations of the generic WEC which differ in their resonance frequencies (0.22 Hz and 0.33 Hz). The energy produced by all control approaches considered in this report is calculated on a predefined interval of the wave recorded at Galway Bay (between 750 s and 1000 s). This allows to perform the identification of the AR model for the wave prediction based on past data -an important real-time constraint.

Some simulation results for the case ω 0 = 0.22 Hz are given in the form of system trajectories. Figure 6 captures the product of control and velocity (instantaneous power) of the PI velocity controls for different values for the PTO efficiency. The lower the efficiency of the PTO, the less power is taken from the net. Obviously, the PI velocity control scales badly with a lower PTO efficiency since both negative effects show their impact: power taken from the net becomes more expensive and the produced power is reduced. This is a big advantage of the latching controls which do not need to draw power from the net such that they are only effected by the lower mechanical to electrical power conversion effiency. This also means that the same latching control law is applicable to different efficiency coefficients η -in contrast to the PI velocity control which strongly depends on η. The Figures 7 and8 show the trajectories of the device velocity together with the incident wave force for the PI velocity control and the threshold unlatching control. The trajectory of wave force and device velocity for the two predictive latching schemes are given in Figure 9. Even though the mean force predictive control is much simpler, it leads to the same η = 0.8 η = 0.9 η = 1 time / s normed u q / -960 980 1000 -1 0 Fig. 6. The product of control and velocity (instantaneous power) of the PI velocity controls for different values for the PTO efficiency.

device motion and thus to the same average power output (in fact it is even slightly higher).

The resulting mean power outputs are captured by Table 1. The power output is always normalized row-wise (for each efficiency coefficient η) with respect to the power output of the PI velocity control. Normalized mean power outputs greater than 1 indicate that the corresponding control does better than the PI velocity control. It is quite interesting to note the different results for the two different WEC.

The only scenario in which the PI velocity control performs better than the latching controls is the case in which the eigenfrequency of the WEC is much higher than the dominant wave frequency (ω 0 = 0.33) with perfect PTO efficiency. By using a lot of reactive power, the PI velocity control can bring the device in resonance with the wave and extract a high amount of energy. As soon as the PTO efficiency drops lower than 1 however, the passive latching control schemes outperform the PI velocity control.

Both predictive latching schemes yield almost the same power output. In the case of the slower WEC (ω 0 = 0.22), the mean force latching control performs slightly better than the model based latching. On first thought, this seems astonishing. However, the reader is reminded that in the end both schemes are globally suboptimal due to the limited predicion horizon. Clearly, the imperfect prediction leads to higher performance losses in the case of the faster WEC (ω 0 = 0.33). The reason for this is that the latching timing is determined based on shorter windows (proportional to the period of the systems dominant eigenfrequency) of the predicion horizons s. t. local errors have a higher impact. The threshold unlatching even outperforms the predictive latching schemes using the AR model predictions in the case of the faster WEC.

The performance gap between threshold unlatching and predictive unlatching is smaller for the faster WEC because the WEC eigenfrequency is further away from the dominant frequencies in the wave spectrum which allows to use higher force thresholds.

SUMMARY

Different predictive and non-predictive latching control schemes were benchmarked against a PI velocity control in an irregular sea. Considering imperfect mechanical to electrical and electrical to mechanical power conversion shows the benefit of the latching controls with respect to the PI velocity feedback. Depending on the ratio between system resonance frequency and dominant wave frequency, the predictive latching schemes are more or less sensitive to imperfect wave force predictions. For fast WEC, whose controls are more sensitive to imperfect wave predictions, the considered non-predictive latching scheme can show off its advantages.

ACKNOWLEGEMENTS

The authors would like to thank the Irish Marine Institute for providing the wave elevation measurements used in this study.

  Fig. 1. Simplified model of a wave energy generator.

Fig. 2 .

 2 Fig. 2. Wave recorded at Galway Bay in the time and frequency domain.

Fig. 3 .

 3 Fig. 3. Results of the line search for k v,l (case ω 0 = 0.22 Hz).

Fig. 4 .

 4 Fig. 4. The principle of mean wave force predictive latching.

  Fig. 5. 6 s-ahead wave prediction based on the identified AR model.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. The incident wave force and the resulting velocity of the PTO: PI velocity control.

Table 1 .

 1 A Babarit and AH Clément. Optimal latching control of a wave energy device in regular and irregular waves. Applied Ocean Research, 28(2):77-91, 2006. Auréelien Babarit, Gaielle Duclos, and AH Clemént. Benefit of latching control for a heaving wave energy device in random sea. In 13th Int Offshore and Polar Engineering Conference, volume 1, pages 341-348, 2003. Aurélien Babarit, Gaelle Duclos, and Alain H Clément. Comparison of latching control strategies for a heaving wave energy device in random sea. Applied Ocean Research, 26(5):227-238, 2004. K. Budal, J. Falnes, L.C. Iversen, P.M. Lillebekken, G. Oltedal, T. Hals, and T. Onshus. The norwegian wave-power buoy project. In The Second International Symposium on Wave Energy Utilization, 1982. H. Eidsmoen. Tight-moored amplitude-limited heavingbuoy wave-energy converter with phase control. Applied Ocean Research, 20(3):157 -161, 1998. ISSN 0141-1187. A. F. O. Falcão. Phase control through load control of oscillating-body wave energy converters with hydraulic pto system. Ocean Engineering, 35(3):358-366, 2008. B. Fischer, P. Kracht, and S. Perez-Becker. Onlinealgortihm using adaptive filters for short-term wave prediction and ist implementation. In Proceedings of the ICOE 2012, 2012. Francesco Fusco and John V Ringwood. Short-term wave forecasting for real-time control of wave energy converters. Sustainable Energy, IEEE Transactions on, 1(2):99-106, 2010a. Francesco Fusco and John V Ringwood. Short-term wave forecasting with AR models in real-time optimal control of wave energy converters. In Industrial Electronics (ISIE), 2010 IEEE International Symposium on, pages 2475-2480. IEEE, 2010b. Rico Hjerm Hansen and Morten M. Kramer. Modelling and Control of the Wavestar Prototype. University of Southampton, 2011. M.P. Schoen, J. Hals, and T. Moan. Wave prediction and robust control of heaving wave energy devices for irregular waves. Energy Conversion, IEEE Transactions on, 26(2):627-638, 2011. ISSN 0885-8969. Performance Comparison.

	WEC model ω 0	PTO eff. η	PI-velocity	thresh. latch.	model based	latch.	model based	latch. AR pred.	mean force	latch.	mean force	latch. AR pred.
		1	1	1.06	1.20	1.11	1.24	1.17
	0.22	0.9	1	1.15	1.31	1.21	1.35	1.28
		0.8	1	1.23	1.40	1.29	1.44	1.36
		1	1	0.86	0.89	0.76	0.89	0.76
	0.33	0.9	1	1.12	1.16	1.00	1.17	1.00
		0.8	1	1.34	1.39	1.19	1.39	1.19

The instantaneous power output can be positive and negative, i. e. power can be generated and consumed.