Vu Trung Hieu

On the solution existence and stability of polynomial optimization problems

Keywords: Polynomial optimization, Regularity condition, Asymptotic cone, Frank-Wolfe type theorem, Eaves type theorem, Upper semicontinuity, Local upper-Hölder stability, Genericity Mathematics Subject Classification (2000) 90C30, 14P10

In this paper, we introduce and investigate a new regularity condition in the asymptotic sense for optimization problems whose objective functions are polynomial. The normalization argument in asymptotic analysis enables us to study the existence as well as the stability of solutions of these problems. We prove a Frank-Wolfe type theorem for regular optimization problems and an Eaves type theorem for non-regular pseudoconvex optimization problems. Moreover, under the regularity condition, we show results on the stability such as upper semicontinuity and local upper-Hölder stability of the solution map of polynomial optimization problems. At the end of the paper, we discuss the genericity of the regularity condition.

Introduction

We consider the following optimization problem minimize f (x) subject to x ∈ K, where K is a nonempty, closed subset of R n and f : R n → R is a polynomial in n variables of degree d ≥ 2. The problem and its solution set are denoted by OP(K, f ) and Sol(K, f ) respectively. Let f d be the homogeneous component of degree d of f , and let K ∞ be the asymptotic cone of K that will be introduced in Section 2. We say that OP(K, f ) is regular if the solution set of the asymptotic problem OP(K ∞ , f d ) is bounded, and the problem is nonregular otherwise. The regularity condition has appeared in studies about the solution existence and stability in quadratic programming (see, e.g., [START_REF] Lee | Quadratic Programming and Affine Variational Inequalities: A Qualitative Study[END_REF][START_REF] Tam | On the solution existence and stability of quadratically constrained nonconvex quadratic programs[END_REF] and the references therein).

Asymptotic cones and functions play an important role in optimization and variational inequalities [START_REF] Auslender | Asymptotic Cones and Functions in Optimization and Variational Inequalities[END_REF]. The normalization argument in asymptotic analysis enables us to study the existence and stability of solutions not only for quadratic programming, linear complementarity problems, and affine variational inequalities (see, e.g., [START_REF] Lee | Quadratic Programming and Affine Variational Inequalities: A Qualitative Study[END_REF][START_REF] Cottle | The Linear Complementarity Problem[END_REF]), but also for polynomial complementarity problems and polynomial variational inequalities that have unbounded constraint sets (see, e.g., [START_REF] Gowda | Polynomial complementarity problems[END_REF][START_REF] Hieu | Solution maps of polynomial variational inequalities[END_REF]). In this paper, the normalization argument is used as the main technique to investigate the existence as well as the stability of solutions to polynomial optimization problems.

In 1956, Frank and Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF] proved that if K is polyhedral and f is quadratic and bounded from below on K, then Sol(K, f ) is nonempty. Several versions of the Frank-Wolfe theorem for quadratic, cubic, and polynomial optimization problems have been shown in [START_REF] Lee | Quadratic Programming and Affine Variational Inequalities: A Qualitative Study[END_REF][START_REF] Tam | On the solution existence and stability of quadratically constrained nonconvex quadratic programs[END_REF][START_REF] Luo | On extensions of the Frank-Wolfe theorems[END_REF][START_REF] Belousov | A Frank-Wolfe type theorem for convex polynomial programs[END_REF][START_REF] Obuchowska | On generalizations of the Frank-Wolfe theorem to convex and quasiconvex programmes[END_REF][START_REF] Dinh | A Frank-Wolfe type theorem for nondegenerate polynomial programs[END_REF][START_REF] Klatte | On a Frank-Wolfe type theorem in cubic optimization[END_REF]. Belousov and Klatte [START_REF] Belousov | A Frank-Wolfe type theorem for convex polynomial programs[END_REF], and Obuchowska [START_REF] Obuchowska | On generalizations of the Frank-Wolfe theorem to convex and quasiconvex programmes[END_REF] have proved Frank-Wolfe type theorems for convex and quasiconvex polynomial optimization problems. Recently, by using a technique from semi-algebraic geometry, Dinh, Ha and Pham [START_REF] Dinh | A Frank-Wolfe type theorem for nondegenerate polynomial programs[END_REF] have shown a Frank-Wolfe type theorem for nondegenerate problems. The present paper gives another Frank-Wolfe type theorem, which says that if OP(K, f ) is regular and f is bounded from below on K, then the problem has a solution. Besides, the Eaves theorem [START_REF] Eaves | On quadratic programming[END_REF] provides us with another criterion for the existence of solutions to quadratic optimization problems. Extensions of this theorem for quadratically constrained quadratic problems have been investigated in [START_REF] Lee | Quadratic Programming and Affine Variational Inequalities: A Qualitative Study[END_REF][START_REF] Tam | On the solution existence and stability of quadratically constrained nonconvex quadratic programs[END_REF][START_REF] Kim | Solution existence and stability of quadratically constrained convex quadratic programs[END_REF][START_REF] Nguyen | Extension of Eaves Theorem for determining the boundedness of convex quadratic programming problems[END_REF]. This paper introduces an Eaves type theorem for non-regular pseudoconvex optimization problems, where the constraint sets are convex.

Under the assumption that the constraint set K is compact and semialgebraic, some stability and genericity results for polynomial optimization problems have been shown by Lee and Pham [START_REF] Lee | Stability and genericity for semi-algebraic compact programs[END_REF]. If K is compact, then its asymptotic cone is trivial, i.e., K ∞ = {0}; Hence that OP(K, f ) satisfies the regularity condition obviously. In the present paper, K may be unbounded. Under the regularity condition, we prove several local properties of the solution map of polynomial optimization problems such as local boundedness and upper semicontinuity. Furthermore, based on an error bound for a polynomial system in [START_REF] Li | New fractional error bounds for polynomial systems with applications to Holderian stability in optimization and spectral theory of tensors[END_REF], we prove the local upper-Hölder stability of the solution map.

We denote by R d [x] the space of all polynomials of degree at most d and by R d the set of all polynomials g of degree

d such that OP(K, g) is regular. The set R d is an open cone in R d [x]. At the end of this work, K is defined by convex polynomials, we prove that R d is generic in R d [x].
The organization of the paper is as follows. Section 2 gives a brief introduction to asymptotic cones, polynomials, and the regularity condition. Section 3 proves two criteria of the solution existence. Section 4 investigates properties of the solution map. The last section discusses the genericity of the regularity condition.

Preliminaries

Recall that the asymptotic cone [START_REF] Auslender | Asymptotic Cones and Functions in Optimization and Variational Inequalities[END_REF] of a nonempty closed subset S in R n is defined and denoted by

S ∞ = v ∈ R n : ∃t k → +∞, ∃x k ∈ S with lim k→∞ x k t k = v .
Clearly, the cone S ∞ is closed and contains 0. The set S is bounded if and only if S ∞ is trivial. Furthermore, if S is convex then S ∞ is a closed convex cone and S ∞ = 0 + S, where 0 + S is the recession cone of S, that consists of all vectors v ∈ R n such that x + tv ∈ S for any x ∈ S and t ≥ 0. Thus, one has S = S + S ∞ when S is convex. Let d ≥ 2 be given. The dimension of the space R d [x] is finite; its dimension is denoted by ρ. Let X(x) be the vector consisting of ρ monomials of degree at most d which is listed by lexicographic ordering

X(x) := (1, x 1 , x 2 , . . . , x n , x 2 1 , x 1 x 2 , . . . , x 1 x n , . . . , x d 1 , x d-1 1 x 2 , . . . , x d n ) T . For every g ∈ R d [x]
, there exists a unique vector a = (a 1 , . . . , a ρ ) ∈ R ρ such that g(x) = a T X(x). We denote by g the 2 -norm of the polynomial g, namely

g := a = a 2 1 + • • • + a 2 ρ . The Cauchy-Schwarz inequality yields |g(x)| ≤ X(x) g . Furthermore, if {g k } is a convergent sequence in R d [x] with g k → g, then g k d → g d .
Throughout the paper, we assume that the constraint set K ⊂ R n is nonempty and closed, and the objective function f : R n → R is a polynomial of degree d ≥ 2.

We say that OP(K, f ) is a polynomial optimization problem if K is given by polynomials. With the given set K and the given integer d, the solution map of polynomial optimization problems OP(K, g), where g ∈ R d [x], is defined by

Sol K (•) : R d [x] ⇒ R n , g → Sol(K, g). Assume that g ∈ R d [x] with deg g = d and g = g d + • • • + g 1 + g 0 ,
where g l is a homogeneous polynomial of degree l, i.e., g l (tx) = t l g l (x) for all t ≥ 0 and x ∈ R n , l ∈ [d] := {1, . . . , d}, and g 0 ∈ R. Then, g d is the leading term (or the recession polynomial) of the polynomial g (of degree d). Clearly, one has

g d (x) = lim λ→+∞ g(λx) λ d , ∀x ∈ R n .
For the pair (K, f ), the asymptotic pair (K ∞ , f d ) is unique. The asymptotic optimization problem OP(K ∞ , f d ) plays a vital role in the investigation of behavior of OP(K, f ) at infinity. The following remarks point out (without proof) the basic properties of the asymptotic problem.

Remark 2.1 Since f d is a homogeneous polynomial and K ∞ is a closed cone, the asymptotic optimization problem OP(K ∞ , f d ) has a solution if and only if

f d is non-negative on K ∞ . Remark 2.2 Assume that Sol(K ∞ , f d ) is nonempty. Then, this set is a closed cone with 0 ∈ Sol(K ∞ , f d ). In addition, Sol(K ∞ , f d ) coincides with the zero set of f d in K ∞ , i.e., Sol(K ∞ , f d ) = {x ∈ K ∞ : f d (x) = 0}.
Now, we introduce the regularity notion concerning the boundedness of the solution set of OP(K ∞ , f d ).

Definition 2.1 The problem OP(K, f ) is said to be regular if Sol(K ∞ , f d ) is bounded and non-regular otherwise.

Denote by E d (resp., O d , U d ) the set of all polynomials g of degree d such that Sol(K ∞ , f d ) is the empty set (resp., the trivial cone, an unbounded cone). Clearly, R d = E d ∪ O d , and one has the following disjoint union:

R d [x] = R d-1 [x] ∪ E d ∪ O d ∪ U d .
(2.1)

Remark 2.3
The boundedness of Sol(K ∞ , f d ) implies that of Sol(K, f ). Indeed, assume to the contrary that Sol(K, f ) is unbounded. There exists an unbounded sequence {x k } ⊂ Sol(K, f ). Without loss of generality, we can assume x k is nonzero for all k, x k → +∞, and

x k -1 x k → x for some x ∈ R n with x = 1. Note that f (x k ) = f * , where f * ∈ R is the minimum of f over K, for all k.
By dividing the last equation by x k d and letting k → +∞, we obtain

f d (x) = 0. It follows that x ∈ Sol(K ∞ , f d ). Since x = 0, the cone Sol(K ∞ , f d )
is unbounded, which contradicts our assumption. Thus, the claim is proved.

Remark 2.4 We observe that the set R d is nonempty. If K is bounded then R d coincides with the set of all g ∈ R d [x] such that deg g = d.
Hence, we can suppose that K is unbounded. Clearly, the cone K ∞ also is unbounded. Let x ∈ K ∞ be nonzero. There exists l ∈ [n] such that xl = 0. Let us define a homogeneous polynomial of degree d as f (x) := -(x l x l ) d . For any t > 0, one has tx ∈ K ∞ , and

f (tx) = -(x 2 l ) d t d → -∞ as t → +∞. Then, f is not bounded from below on K ∞ . This yields Sol(K ∞ , f ) = ∅ and f ∈ R d . Example 2.1 Consider the case that n = 1, K = R and d = 2. One has R 2 [x] = {a 2 x 2 +a 1 x+a 0 : (a 2 , a 1 , a 0 ) ∈ R 3 }. Since K ∞ = R, an easy computation shows that E 2 = {a 2 x 2 + a 1 x + a 0 : a 2 < 0, a 1 ∈ R, a 0 ∈ R}, O 2 = {a 2 x 2 + a 1 x + a 0 : a 2 > 0, a 1 ∈ R, a 0 ∈ R}, and R 2 = {a 2 x 2 + a 1 x + a 0 : a 2 = 0, a 1 ∈ R, a 0 ∈ R}.

Two criteria for the solution existence

We now introduce two criteria for the solution existence of OP(K, f ). In the proofs, the normalization argument in asymptotic analysis plays a vital role; meanwhile, the semi-algebraicity of K is not required.

A Frank-Wolfe type theorem for regular problems

The following theorem provides us a criterion for the solution existence of regular optimization problems. Theorem 3.1 (Frank-Wolfe type theorem) If OP(K, f ) is regular and f is bounded from below on K, then its solution set is nonempty and compact.

Proof Suppose that f ∈ R d , i.e. f ∈ E d ∪ O d , and there exists γ ∈ R such that γ ≤ f (x), for all x ∈ K. For any given v ∈ K ∞ , there are two sequences

{t k } ⊂ R + and {x k } ⊂ K such that t k → +∞ and t -1 k x k → v as k → +∞. For any k, one has γ ≤ f (x k ). Dividing both sides of the last inequality by t d k and letting k → +∞, we obtain 0 ≤ f d (v). Thus, f d is non-negative over K ∞ . It follows from Remark 2.1 that f does not belong to E d ; hence, we conclude that f must be in O d .
Let x ∈ K be given, and

M := {x ∈ K : f (x) ≤ f (x)}. It is easy to check that Sol(M, f ) = Sol(K, f ).
Hence, we need only to prove that Sol(M, f ) is nonempty and compact.

Clearly, M is closed. We claim that M is bounded. On the contrary, we suppose that there exists an unbounded sequence {x k } ⊂ M such that x k is nonzero for all k, x k → +∞, and

x k -1 x k → v for some v ∈ R n with v = 1. One has γ ≤ f (x k ) ≤ f (x), (3.1)
for all k. Dividing the values in (3.1) by x k d and letting k → +∞, we get

f d (v) = 0. This yields v ∈ Sol(K ∞ , f d ). Because of v = 0, Sol(K ∞ , f d ) is
unbounded. This contradicts our assumption and, thus, the claim is proved.

The compactness of M and Bolzano-Weierstrass' Theorem allow us to conclude that Sol(M, f ) is nonempty and compact. Remark 3.1 From the proof of Theorem 3.1, we see that if f is bounded from below on

K then Sol(K ∞ , f d ) is nonempty, i.e. f ∈ O d ∪ U d . Hence, if f ∈ E d then OP(K, f ) has no solution. Corollary 3.1 Assume that f = α 1 x d 1 +• • •+α n x d n +p where d is even, α > 0 for all ∈ [n], p is a polynomial with deg p < d. Then, Sol(K, f ) is nonempty and compact. Proof Clearly, f d = α 1 x d 1 + • • • + α n x d n is non-negative over R n . It follows that f d is also non-negative over K ∞ . From Remarks 2.1 and 2.2, it is clear that Sol(K ∞ , f d ) is nonempty and Sol(K ∞ , f d ) = {x ∈ K ∞ : α 1 x d 1 + • • • + α n x d n = 0} = {0}. This means that f ∈ O d .
Clearly, f is bounded from below on K, and the condition of Theorem 3.1 holds. Therefore, Sol(K, f ) is nonempty and compact.

The following example illustrates Theorem 3.1, in which the constraint set is neither convex nor semi-algebraic.

Example 3.1 Consider the optimization problem OP(K, f ), where the polynomial f is given by f (x 1 , x 2 ) = x 3 2 -x 1 x 2 and the constraint set K is given by K

= {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 -x 1 ≥ 0, e x1 -x 2 ≥ 0}. Since f 3 (x 1 , x 2 ) = x 3 2 and K ∞ = {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 -x 1 ≥ 0}, one has Sol(K ∞ , f 3 ) = {(0, 0)}.
According to Theorem 3.1, Sol(K, f ) is nonempty and compact.

An Eaves type theorem for non-regular problems

In this subsection, we investigate the solution existence of non-regular optimization problems, where the objective functions are pseudoconvex on the constraint sets.

Assume that U is an open subset of R n . One says the polynomial f is pseudoconvex on U if, for any x, y ∈ U such that ∇f (x), y -x ≥ 0, here ∇f is the gradient of f , we have f (y) ≥ f (x). Recall that f is pseudoconvex on U if and only if ∇f is pseudomonotone on U [17, Theorem 3.1], i.e. if, for any x, y ∈ U such that ∇f (x), y -x ≥ 0, we have ∇f (y), y -x ≥ 0.

Lemma 3.1 Assume that K is convex and f is pseudoconvex on an open set U containing K. If x 0 ∈ Sol(K, f ), then ∇f (x), x -x 0 ≥ 0 for all x ∈ K.
Proof Since f is pseudoconvex on the set U , the gradient ∇f is pseudomonotone on U . Suppose that x 0 ∈ Sol(K, f ), one has ∇f (x 0 ), x -x 0 ≥ 0 for all x ∈ K (see, e.g., [START_REF] Ansari | Generalized Convexity, Nonsmooth Variational Inequalities, and Nonsmooth Optimization[END_REF]Proposition 5.2]). The pseudomonotonicity of the gradient implies that ∇f (x), x -x 0 ≥ 0 for all x ∈ K. The lemma is proved. Proof Suppose that OP(K, f ) is non-regular. We prove (a) ⇒ (b). Assume that (a) holds. For each k ∈ N, we denote

K k = {x ∈ R n : x ∈ K, x ≤ k}.
Clearly, K k is compact and convex. Without loss of generality, we can assume that K k is nonempty. According to Bolzano-Weierstrass' Theorem, OP(K k , f ) has a solution, denoted by x k .

We assert that the sequence {x k } is bounded. Indeed, suppose on the contrary that {x k } is unbounded, x k = 0 for all k, x k → +∞, and

x k -1 x k → v, where v ∈ K ∞ and v = 1. For each k, one has f (x k ) ≤ f (x), ∀x ∈ K k . (3.2)
Let y ∈ K be given. For k large enough, y ∈ K k and f (x k ) ≤ f (y). By dividing two sides of the last inequality by x k d and letting k → +∞, we obtain

f d (v) ≤ 0. This leads to v ∈ Sol(K ∞ , f d ) \ {0}. Furthermore, since f is pseudoconvex on K k , from Lemma 3.1 we have ∇f (y), y -x k ≥ 0. (3.3)
Dividing both sides of the inequality in (3.3) by x k and letting k → +∞, we obtain ∇f (y), v ≤ 0. The conclusion holds for any x ∈ K, i.e., ∇f (x), v ≤ 0 for all x ∈ K. This contradicts to our assumption. Hence, {x k } is bounded. We can assume that x k → x. From (3.2), by the continuity of f , it easy to check that x solves OP(K, f ), so Sol(K, f ) is nonempty.

To prove the compactness of the solution set, we can repeat the previous argument by supposing that there is an unbounded solution sequence {x k }, and can show that there exists v ∈ Sol(K ∞ , f d ) \ {0} such that ∇f (x), v ≤ 0 for all x ∈ K. This contradicts to (b).

(b) ⇒ (a) Since K is convex, one has K ∞ = 0 + K and K = K + K ∞ . Suppose that Sol(K, f ) is nonempty and compact, but (b) is wrong, i.e. there exists v ∈ Sol(K ∞ , f d ) \ {0} such that ∇f (x), v ≤ 0 for all x ∈ K. Let x be a solution of OP(K, f ). For any t ≥ 0, one has x+tv ∈ K and ∇f (x+tv), v ≤ 0. Thus, we have ∇f (x + tv), x -(x + tv) ≥ 0.

The pseudoconvexity of f yields f (x) ≥ f (x 0 + tv). Hence, x 0 + tv belongs to Sol(K, f ), for any t ≥ 0. This shows that Sol(K, f ) is unbounded which contradicts to our assumption. Thus (a) holds, and the proof is complete. Proof Since the convexity implies the pseudoconvexity, by applying Theorem 3.2 for the convex polynomial f , we have the assertion.

The following example illustrates Corollary 3.2.

Example 3.2 Consider the polynomial optimization problem OP(K, f ) with the objective function f (x 1 , x 2 ) = 1 6 x 3 2 + 1 2 x 2 1 -x 1 x 2 and the constraint set

K = {(x 1 , x 2 ) ∈ R 2 : x 1 x 2 ≥ 1, x 2 ≥ 2}.
The gradient and the Hessian matrix of f , respectively, are given by

∇f = x 1 -x 2 -x 1 + 1 2 x 2 2 , H = 1 -1 -1 x 2 .
It is easy to check that K is convex and

H is positive semidefinite on the open set U = {(x 1 , x 2 ) ∈ R 2 : x 1 x 2 > 0, x 2 > 1} ⊃ K; hence f is convex on K. One has K ∞ = R 2 + and f 3 (x 1 , x 2 ) = 1 6 x 3 2 . This yields Sol(K ∞ , f 3 ) = {(x 1 , x 2 ) ∈ R 2 : x 1 ≥ 0, x 2 = 0}.
For every v = (α, 0) in Sol(K ∞ , f 3 )\{0}, one has α > 0. By choosing the point x = (3, 2) in the constraint set, we get ∇f (x), v = α > 0. Finally, according to Corollary 3.2, the solution set of OP(K, f ) is nonempty and compact.

Stability of the solution map

We investigate the local boundedness, the upper semicontinuity, and the local upper-Hölder stability of the solution map under the regularity condition.

Upper semicontinuity of the solution map

To prove the local boundedness and the upper semicontinuity of the solution map, we need the following lemma.

Lemma 4.1 The set R d is open in R d [x].
Proof To prove the openness of R d , we only need to show that the complement

R d [x]\R d is closed. Clearly, R d [x]\R d = R d-1 [x]∪U d . Let {g k } be a sequence in R d [x] \ R d such that g k → g. From the definition of R d , if deg g < d, i.e., g ∈ R d-1 [x], then g ∈ R d [x] \ R d .
} such that x k ∈ Sol(K ∞ , g k d ), x k → +∞, x k -1 x k → x with x = 1. Let v ∈ K ∞ be given. One has x k v ∈ K ∞ and g k d ( x k v) ≥ g k d (x k
), for any k. Dividing the last inequality by x k d and letting k → +∞, one has

g k d (v) ≥ g k d (x). This conclusion holds for every v ∈ K ∞ . This yields x ∈ Sol(K ∞ , g d ). As x = 1, we have x = 0. It follows that g belongs to U d . The closedness of R d [x] \ R d is proved.
Recall that a set-valued map Ψ : R m ⇒ R n is locally bounded at ū if there exists an open neighborhood U of ū such that ∪ u∈U Ψ (u) is bounded [START_REF] Rockafellar | Variational Analysis[END_REF]Definition 5.14]. The map Ψ is upper semicontinuous at ū ∈ T if for any open set V ⊂ R n such that Ψ (ū) ⊂ V there exists a neighborhood U of ū such that Ψ (u) ⊂ V for all u ∈ U . Recall that if Ψ is closed, namely, the graph Assume to the contrary that O ε is unbounded. Then, there exists an unbounded sequence {x k } and a sequence {g k } ⊂ B(ε, d) such that x k solves OP(K, f + g k ) with x k = 0 for every k, x k → +∞, and x k -1 x k → x with x = 1. By the compactness of B(ε, d), without loss of generality, we can assume that g k → g with g ∈ B(ε, d).

gph(Ψ ) := (u, v) ∈ R m × R n : v ∈ Ψ (u) is closed in R m × R n ,
From assumptions, for every k, one has

(f + g k )(y) ≥ (f + g k )(x k ), (4.3)
for any y ∈ K. Let y ∈ K be fixed and assume that v ∈ K ∞ . By the convexity of K, one has y + x k v ∈ K for any k. From (4.3), we conclude that

(f + g k )(y + x k v) ≥ (f + g k )(x k ).
Dividing this inequality by x k d and taking k → +∞, we obtain

(f + g) d (v) ≥ (f + g) d (x).
The conclusion hold for any v

∈ K ∞ . It follows that x ∈ Sol (K ∞ , (f + g) d ). Because of (4.2), Sol (K ∞ , (f + g) d ) is contained in {0}, which contradicts to x = 1. Hence, O ε must be bounded. (b) It is not difficult to prove that the graph gph(Sol) := (g, x) ∈ R d [x] × R n : x ∈ Sol(K, g) is closed in R d [x] × R n . Since Sol K (•) is locally bounded on R d , according to [19, Theorem 5.19], Sol K (•) is upper semicontinuous at f .

Local upper-Hölder stability of the solution map

When the constraint set K is convex and given by polynomials, we can investigate the local upper-Hölder stability of the solution map under the regular condition. To prove the stability, we need the following lemma.

Lemma 4.2 ([20]

) Let U be a semi-algebraic subset in R n , represented by

U = {x ∈ R n : u i (x) = 0, i ∈ [l], v j (x) ≤ 0, j ∈ [m]} ,
where u i (x), i ∈ [l], and v j (x), j ∈ [m], are polynomials. For any compact set V ⊂ R n , there are constants c > 0 and H > 0 such that

d(x, U ) ≤ c l i=1 |u i (x)| + m j=1 [v j (x)] + H ,
for all x ∈ V , here [r] + := max{0, r} and d(x, U ) the usual distance from x to the set U .

Theorem 4.1 Assume that OP(K, f ) is regular and K is a convex set given by

K = {x ∈ R n : p i (x) = 0, i ∈ [l], q j (x) ≤ 0, j ∈ [m]} ,
where all p i , q j are polynomials. If Sol(K, f ) is nonempty, then the map Sol K (•) is locally upper-Hölder stable at f , i.e., there exist > 0, H > 0 and ε > 0

such that Sol(K, g) ⊂ Sol(K, f ) + g -f H B, (4.4
)

for all g ∈ R d [x] satisfying g -f < ε, where B is the closed unit ball in R n .
Proof Suppose Sol(K, f ) is nonempty and its optimal value is f * . Since OP(K, f ) is regular and K is convex, according to Proposition 4.1, there exists ε > 0 such that Sol(K, f ) ⊂ O ε , defined by (4.1), is bounded. Let V be the closure of O ε . It follows that V is a nonempty compact set. By the assumptions, we see that

Sol(K, f ) = {x ∈ R n : f (x) -f * = 0, p i (x) = 0, i ∈ [l], q j (x) ≤ 0, j ∈ [m]}.
From this equality, by applying Lemma 4.2 for U = Sol(K, f ) and the compact set V , there are constants c 0 > 0 and H > 0 such that

d(x, Sol(K, f )) ≤ c 0 A(x) H ∀x ∈ V, (4.5)
where

A(x) := |f (x) -f * | + l i=1 |p i (x)| + m j=1 [q j (x)] + . Let g ∈ R d [x]
be arbitrary given such that g -f < ε. From the definition of V , Sol(K, f ) and Sol(K, g) are subsets of V . Here, Sol(K, g) may be empty. By the compactness of V , we define the constant L := max{ X(x) : x ∈ V }. Hence, one has

|g(x) -f (x)| ≤ X(x) g -f ≤ L g -f ∀x ∈ V. (4.6)
If Sol(K, g) is empty, then (4.4) is obvious. Thus, we consider the case that Sol(K, g) = ∅. Since both Sol(K, f ) and Sol(K, g) are nonempty and compact, for any x g ∈ Sol(K, g), there is x f ∈ Sol(K, f ) such that

x g -x f = d(x g , Sol(K, f )). ( 4.7)
Because of p i (x g ) = 0 for i ∈ [l] and q j (x g ) ≤ 0 for j ∈ [m], from the definition of A(x), one has A(x g ) = |f (x g ) -f * |. By (4.7) and (4.5), we see that

x g -x f ≤ c 0 A(x g ) H = c 0 |f (x g ) -f * | H . Since x f ∈ Sol(K, f ), we have f (x f ) = f * ≤ f (x g )
. Therefore, we obtain

x g -x f ≤ c 0 |f (x g ) -f * | H = c 0 (f (x g ) -f (x f )) H . (4.8)
It follows from x g ∈ Sol(K, g) that g(x g ) -g(x f ) ≤ 0. Since x g , x f ∈ V , we conclude from (4.6) that

f (x g ) -f (x f ) = (f (x g ) -g(x g )) + (g(x g ) -g(x f )) + (g(x f ) -f (x f )) ≤ (f (x g ) -g(x g )) + (g(x f ) -f (x f )) ≤ 2L g -f .
The inequality (4.8) and the last result lead to

x g -x f ≤ c 0 (2L) H g -f H , consequently, d(x g , Sol(K, f )) = x g -x f ≤ g -f H ,
where = c 0 (2L) H . The conclusion holds for any x g in Sol(K, g). Hence, the inclusion (4.4) of the theorem is proved.

Genericity of the regularity condition

In this section, we discuss the genericity of the regularity condition of polynomial optimization problems.

A subset

A is called generic in R m if A contains a countable intersection of dense and open sets in R m . If A is generic in R m and A ⊂ B then B also is generic in R m . Let T be a topological space. It is known that if h : R m → T is a homeomorphism and A is generic in R m , then h(A) is generic in T .
Let U ⊂ R m be a semi-algebraic set. Then, there exists a decomposition of U into a disjoint union [START_REF] Bochnak | Real Algebraic Geometry[END_REF]Theorem 2.3.6], U = ∪ s i=1 U i , where each U i is semi-algebraically homeomorphic to (0, 1) di . Here, let (0, 1) 0 be a point, (0, 1) k ⊂ R k be the set of points x = (x 1 , . . . , x k ) such that x j ∈ (0, 1) for all j ∈ [k]. The dimension of U is defined by dim(U ) := max{d 1 , . . . , d s }. The dimension is well-defined and does not depends on the decomposition of U . Recall that if U is nonempty and dim(U ) is zero, then U has finitely many points. Furthermore, if dim(R m \U ) < m, then U is generic in R m (see, e.g. [START_REF] Dang | Well-posedness in unconstrained polynomial optimization problems[END_REF]Lemma 2.3]).

The space generated by all monomials of degree d listed by lexicographic ordering {x d 1 , x d-1

1 x 2 , x d-1 1 x 3 , . . . , x d n } is denoted by H d . One has the direct sum R d [x] = H d ⊕ R d-1 [x]. The dimension of H d is denoted by η. For every homogeneous polynomial h ∈ H d , one has a unique b ∈ R η , such that h(x) = b T X d (x), where X T d (x) = (x d 1 , x d-1 1 x 2 , x d-1 1 x 3 , . . . , x d n ). Here, ∇(b T X d (x)) is the gradient vector of b T X d (x) and D b [∇(b T X d (x))] is the Jacobian matrix of b T X d (x) with respect to b. Lemma 5.1 One has rank(D b [∇(b T X d (x))]) = n for all x ∈ R n \{0}.
Proof In the proof, we are only interested in the monomials x d-1 i

x j , where i, j ∈ [n]. Hence, for convenience, we rewrite X T d (x) and b T respectively as follows:

(x d 1 , x d-1 1 x 2 , . . . , x d-1 1 x n ; . . . ; x d-1 n x 1 , x d-
T X d (x) = j∈[n] b 1j x d-1 1 x j + • • • + j∈[n] b nj x d-1 n x j + Q, (5.1)
where Q is a homogeneous polynomial of degree d. From (5.1), an easy computation shows that

∂(b T X d (x)) ∂x i = db ii x d-1 i + (d -1) j =i b ij x d-2 i x j + j =i b ji x d-1 j + ∂Q ∂x i , and the n × η-matrix D b [∇(b T X d (x))
] can be described as follows

D b [∇(b T X d (x))] = M 1 , M 2 , • • • , M n , • • • ,
where the submatrix M i , for i ∈ [n], is defined by

M i =    x d-1 i I i-1 O i×1 O i×(n-i) L 1×(i-1) dx d-1 i R 1×(n-i) O i×(n-i) O (n-i)×1 x d-1 i I n-i    with I k being the unit k × k-matrix, O k×s being the zero k × s-matrix, L 1×(i-1) = (d -1)x d-2 i x 1 , . . . , (d -1)x d-2 i x i-1 , and R 1×(n-i) = (d -1)x d-2 i x i+1 , . . . , (d -1)x d-2 i x n .
We observe that det(

M i ) = dx d(d-1) i
, for all i ∈ [n]. Since x = 0, there exists l ∈ [n] such that x l = 0. This implies that rank(M l ) = n. Hence, the rank of

D b [∇(b T X d (x))] is n, for any x = 0.
Suppose that C is a polyhedral cone given by

C = {x ∈ R n : Ax ≥ 0} , (5.2)
where A = (a ij ) ∈ R p×n . Let KKT(C, g), where g ∈ R d [x], be the set of the Karush-Kuhn-Tucker points of OP(C, g), i.e., x ∈ KKT(C, g) if and only if there exists λ ∈ R p such that

∇g(x) -A T λ = 0, λ T (Ax) = 0, λ ≥ 0, Ax ≥ 0. ( 5.3)
From the Karush-Kuhn-Tucker conditions, we see that Sol(C, g) ⊂ KKT(C, g) for all g ∈ R d [x].

For each index set α ⊂ [p], we associate the pseudo-face C α of C, which is denoted and defined by

C α := x ∈ R n : n j=1 a ij x j = 0 ∀i ∈ α, n j=1 a ij x j > 0 ∀i ∈ [p] \ α ,
where a ij is the element in the i-th row and the j-th column of A. The number of pseudo-faces of C is finite. These pseudo-faces establish a disjoint decomposition of C. So, we obtain KKT(C, g) =

Perspectives

The regularity condition enables us to investigate the stability of the optimal value function of polynomial optimization problems. Furthermore, the regularity condition is useful to study the connectedness of the solution sets of convex polynomial vector optimization problems.

Theorem 3 . 2 (

 32 Eaves type theorem) Assume that K is convex and f is pseudoconvex on an open set containing K. If OP(K, f ) is non-regular, then the following statements are equivalent: (a) If v ∈ Sol(K ∞ , f d ) \ {0}, then there exists x ∈ K such that ∇f (x), v > 0; (b) Sol(K, f ) is nonempty and compact.

Corollary 3 . 2

 32 Assume that K is convex and f is convex on an open convex set containing K. If OP(K, f ) is non-regular, then the following statements are equivalent: (a) If v ∈ Sol(K ∞ , f d ) \ {0}, then there exists x ∈ K such that ∇f (x), v > 0; (b) Sol(K, f ) is nonempty and compact.

  and locally bounded at u, then Ψ is upper semicontinuous at u [19, Theorem 5.19]. Proposition 4.1 Assume that K is convex. If OP(K, f ) is regular, then the following statements hold: (a) The solution map Sol K (•) is locally bounded at f , i.e., there exists ε > 0 such that the set O ε := g∈B(ε,d) Sol(K, f + g), (4.1) where B(ε, d) is the open ball in R d [x] with center 0 and radius ε, is bounded. (b) The solution map Sol K (•) is upper semicontinuous at f . Proof (a) According to Lemma 4.1, R d is open in R d [x]. There is a closed ball B(ε, d) such that f + B(ε, d) ⊂ R d . (4.2)

From Lemma 5 . 1 ,Example 5 . 1

 5151 α⊂[p](KKT(C, g) ∩ C α ) ,(5.4)The following proposition shows that the Karush-Kuhn-Tucker set-valued map of homogeneous polynomial optimization problemsKKT C : R η ⇒ R n , b → KKT C (b) = KKT(C, b T X d (x)), is finite-valued, i.e., the cardinal # KKT C (b) is finite, on a generic semialgebraic set of R η .Proposition 5.1 Assume that C is a polyhedral cone given by (5.2) and the matrix A is full rank. Then, there exists a generic semi-algebraic setS ⊂ R η such that # KKT C (b) < ∞ for any b ∈ S.Proof Let C α be a nonempty pseudo-face of C and 0 / ∈ C α . This implies that X d (x) is nonzero on this pseudo-face. We consider the functionΦ α : R η ×C α × R |α| + → R n+|α| ,which is defined byΦ α (b, x, λ α ) = ∇(b T X d (x)) + i∈α λ i A i , A α x ,where A α x = (A i1 x, . . . , A i |α| x ), i j ∈ α. Clearly, C α is smooth and Φ α is a semi-algebraic function of class C ∞ . The Jacobian matrix of Φ α is determined as followsDΦ α = D b [∇(b T X d (x))] for all x ∈ C α , the rank of D b [∇(b T X d (x))] is n.Since the rank of A α is |α|, we conclude that the rank of the matrix DΦ α is n + |α| for all x ∈ C α . Therefore, 0 ∈ R n+|α|+|J| is a regular value of Φ α . According to the Sard Theorem with parameter [22, Theorem 2.4], there exists a generic semi-algebraic setS α ⊂ R η such that if b ∈ S α then 0 is a regular value of the map Φ α,b : C α × R |α| → R n+|α| , Φ α,b (x, λ α ) = Φ α (b, x, λ α ).Proof From Remark 5.1, the recession cone K ∞ is a nonempty polyhedral cone, where K ∞ = {x ∈ R n : Ax ≥ 0}. According to Corollary 5.1, there exists a generic setG d in H d such that # Sol(K ∞ , g) < ∞ for any g ∈ G d . Because of the direct sum R d [x] = H d ⊕ R d-1 [x], the set G d ⊕ R d-1 [x] is generic in R d [x]. It is easy to check that G d ⊕ R d-1 [x] ⊂ R d .Hence, R d is generic in R d [x]. Consider the problem OP(K, f ) given in Example 2.1, we see that R 2 = {a 2 x 2 + a 1 x + a 0 : a 2 = 0, a 1 ∈ R, a 0 ∈ R} is open and dense in R 2 [x].

  Thus, we can suppose that deg g = d. One has g k d → g d , here g k d is the component of degree d of g k . We now prove that g belongs to U d . For each k, Sol(K ∞ , g k d ) is unbounded. There exists an unbounded sequence {x k

  1 n x 2 , . . . , x d n ; . . . ) and (b 11 , b 12 , . . . , b 1n ; . . . ; b n1 , b n2 , . . . , b nn ; . . . ).

	Then, we have
	b

Acknowledgements The author would like to thank the anonymous referees for their corrections and comments. This work has been supported by European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Actions, grant agreement 813211 (POEMA).

We see that Ω(α, b) := Φ -1 α,b (0) is a semi-algebraic set. From the Regular Level Set Theorem [START_REF] Tu | An Introduction to Manifolds[END_REF]Theorem 9.9], we can claim that if Ω(α, b) is nonempty then the set is a 0-dimensional semi-algebraic set. It follows that Ω(α, b) is a finite set. Moreover, from (5.3)

We consider the case that 0 ∈ C α and define U := C α \ {0}. As is clear, U is semi-algebraic. From (5.3), we see that 0 ∈ KKT C (b). Hence,

where q 1 , . . . , q m are convex polynomials, then the recession cone of K is a nonempty polyhedral cone. We denote

For each j ∈ [m], K j is closed convex set, and K j ∞ is polyhedral (see [9, p.39]).

If follows that K ∞ is a nonempty polyhedral cone. Hence, there exists a matrix


Peshal Nayak

DS-CDMA/CI Based Communication Technique for Underwater Unmanned Vehicles

An increasing demand to explore the vast ocean floor for various purposes has made it necessary to deploy networks of underwater unmanned vehicles (UUVs) capable of carrying out multiple tasks. For information transfer among these vehicles, it is essential to design efficient communication systems and techniques which reduce the transmission error thereby improving system performance. The network consists of underwater components which perform activities such as communicating with each other using sound waves, gathering information and transmitting this information to the surface buoys which forward the data to a receiver onshore or on board some ship. This report studies the use of Carrier Interferometry (CI) spreading codes in underwater networks which use DS-CDMA technique for transmitting data. A network model has been presented and performance of DS-CDMA/CI has been studied for this model through simulations. The performance of this system has also been compared to those networks using Walsh-Hadamard (WH) codes to prove that DS-CDMA/CI is a better transmission scheme as compared to DS-CDMA/WH. Furthermore, multirate transmission in underwater environment has also been explored with the use of DS-CDMA/CI simulations. The results have shown that DS-CDMA/CI is an excellent candidate for future single rate as well as multirate transmission schemes.

I. INTRODUCTION

There has been a continuous growth in the demand for communication networks designed for activities involving exploration as well as exploitation of the ocean resources such as information collection, environment monitoring, predicting natural disasters by observation of seismic activity in the ocean, providing assistance to navigating ships, etc. Recent search efforts for the missing Malaysian Airlines flight MH370 included an underwater search for its debris with the help of an Unmanned vehicle named Bluefin-21. Due to availability of reliable acoustic modems, use of unmanned underwater vehicles (UUV) has been considered as a potential candidate for such tasks [START_REF] Stojanovic | Multiuser undersea acoustic communications in the presence of multipath propagation[END_REF]. An underwater network consists of various devices which may be submerged or floating. Each of them have been briefly explained below:

An Unmanned Underwater Vehicle (UUV) is a device which gets driven by a propulsion system and is controlled and instructed by an onboard processing device. It has the ability to maneuver itself in all the three directions. As a result of this level of control, the vehicle is able to precisely follow numerous programmed trajectories at all possible times. The vehicle is equipped with sensors which sample the ocean as the vehicle moves through it making measurements in both space and time.

There are various modes in which a UUV operates. Commonly used ones are passive floatation at a specified depth or density layer in the sea or actively loiter in the vicinity of the desired location, vary depth as it moves through the water or swim at a constant pressure level. Thus, they are a highly productive means of performing seafloor surveys using acoustic or optical imaging systems.

As the need for sensor as well as control power increases, UUV's become the preferred choice for such operations. A UUV can stay underwater for a sufficient amount of time in the range of 8-50 hours. Most vehicles can vary their velocity between 0.5 and 2.5 m/s. The optimum speed and the resulting greatest range of the vehicle occurs when its hotel load (all required power except propulsion) is twice the propulsion load. For most vehicles, this occurs at a velocity of around 1.5 m/s. Most UUVs are also equipped for surface communication in order to deal with cases wherein they get lost. During such transmissions, the antenna used in these systems communicates with the nearest surface component at 2.4 GHz so it can be traced and the information can be retrieved.

Typically, there is a network operating underwater to gather information which is passed on to a floating device mounted on a surface buoy which then transmits the signal over the air. The surface buoy enables flotation of all the devices necessary for reception of acoustic signal from the underwater components, transducers for conversion of sound signal into electrical signal which then gets radiated as an electromagnetic wave and other miscellaneous components. Advanced terminals set up on such buoys also contain hardware to provide the facility of setting up a Wireless Local Area Network (WLAN) 1 with other buoys or nearby receivers [START_REF] Hodgkiss | A high frequency phase conjugation array[END_REF]. The Wireless LAN operates at a frequency of 2.4 GHz. However, it is necessary to ensure that the weight and size of components loaded on the buoy does not exceed specific limits. Size reduction may also facilitate installation of additional hardware devices thereby increasing the facilities on the surface component. In this sense, the antenna in [START_REF] Nayak | A novel compact tri-band antenna design for WiMax, WLAN and bluetooth applications[END_REF] becomes important from the point of view of size reduction. The antenna is compact, operates at 2.4 GHz and has good radiation properties as well as radiation efficiency which meets all the requirements for the antennas used in surface buoys. Another benefit of using the antenna is size reduction which enables addition of extra hardware in the space saved and thus additional facilities on the buoy 2 .

Due to the broadcast nature of the channel, underwater acoustic communication among multiple devices requires a multiple access communication techniques. In a multiple access channel, multiple transmitters send signals to a common receiver. Examples of multiple access communication systems include mobile to base station communication of cellular systems or the telephone lines which are time multiplexed between multiple voice users. The transmitter power for each of the transmitter may be different. The main objective of multiple user communications is to facilitate efficient sharing of limited resources while maintaining lower levels of interference to each other. Most common techniques used are Frequency-Division Multiple Access (FDMA), Time-Division Multiple Access (TDMA) and Code-Division Multiple Access (CDMA). These techniques allocate resources to individual communication devices such that an efficient information transfer can be achieved.

In case of FDMA systems, the frequency band available is divided into several small sub-bands and each device communicates within the sub-band allocated to it. The channels are typically narrowband so they experience flat-fading as a result of which there is no requirement for compensation techniques for Inter Symbol Interference (ISI). Considering the low coherence bandwidth and total communication bandwidth of a UWA system [START_REF] Sozer | Underwater acoustic networks[END_REF], it is likely that as the number of device increase the transmission bandwidth for individual devices will fall below the coherence bandwidth of the channel thus making the device vulnerable to fading which is prominent in an underwater acoustic channel. Also, due to fixed bandwidth allocation, such systems prove inefficient in scenarios involving bursty traffic [START_REF] Rappaport | Wireless Communications[END_REF].

TDMA, on the other hand, allocates time slots for individual devices by dividing the total time frame. Thus, every user is assigned a cyclically repeating time slot for transmission and reception. As a result, transmission is not continuous and digital transmission techniques are required to allow for buffering. Although the resultant higher transmission bandwidth makes the system more resistant to frequency selective fading, inclusion of guard intervals to prevent packet collision in adjacent time slots results in greater overhead in comparison to FDMA. Thus, in case of underwater communications, CDMA gains an advantage over these techniques as it has the ability to exploit multipath-diversity with the use of RAKE receiver. Due to this feature, it has become the chief focus of engineers working on underwater acoustic communication [START_REF] Tsimenidis | Spreadspectrum based adaptive array receiver algorithms for the shallow-water acoustic channel[END_REF].

In this report, the use of Direct Sequence (DS) CDMA for achieving an efficient communication among several underwater unmanned vehicles with the use of Carrier Interferometry spreading codes has been studied.

II. UNDERWATER ACOUSTIC CHANNEL

Underwater channels are generally known to be the most difficult communication medium in use today. Use of acoustic waves is not the only way of transmitting information wirelessly in an underwater medium. Attempts have also been made in the past to use radio and optical waves. However, in order for a radio wave to travel a sufficient amount of distance through the saline water of the sea, its frequency has to be extra low. The typical range of such frequencies is 30 Hz to 300 Hz. Transmission at such low frequencies requires large antenna and higher transmission powers. Optical waves, which is the second option, do not suffer a lot of attenuation. However, they are seriously affected by scattering. Consequently, the narrow laser beam being used for optical signals needs high pointing precision. Although, laser technology is still being perfected for practical use in underwater environment, acoustic waves remain the best solution for signal transmission underwater.

A. Channel Characteristics

Being an extremely complex communication media, shallow water acoustic environment is characterized by frequency dependent attenuation, multipath propagation and a speed of sound which is five orders of magnitude lower (approximately 1500m/s) as compared to that of the electromagnetic waves used in radio transmissions [START_REF] Jiang | Underwater Acoustic Networks-Issues and Solutions[END_REF]. Furthermore, every tap gain in the sparse channel response exhibits a time-varying nature accompanied by motion induced Doppler effect. Thus, the underwater channel combines the worst properties of radio channels such as poor quality of physical link as in the case of a mobile terrestrial communication channel and large latency which characterises a satellite channel.

The frequency dependent attenuation in an underwater acoustic channel results from the loss of acoustic energy in the form of heat which gets absorbed by the channel. The loss is given by

γ(l, f ) = K 0 × l m × a(f ) l
where K 0 is a scaling constant, l is the transmission distance, a(f ) the absorption constant and m is an exponent which models the spreading loss with a typical value between 1 and 2 [START_REF] Stojanovic | Underwater acoustic communication channels: Propagation models and statistical characterization[END_REF].

The noise in an underwater channel consists of two types: ambient noise and site specific noise. Ambient noise is present in the background of a quiet deep sea all the time while sitespecific noise is present only in specific places. Examples of such locations are the polar regions which experience ice cracking resulting in acoustic noise or the snapping shrimp in the case of warm water. Although noise is often approximated as Gaussian, in case of underwater channel it is not white. Multipath propagation in UWA channel occurs primarily due to reflection from the surface and the bottom as well as sound refraction occurring due to variation in sound speed with temperature, salinity and pressure which change depending on the depth and the location. Changes in the underwater medium due to natural phenomena such as those caused by surface waves and other sources such as transmitter and receiver Fig. 1. Impulse response of the underwater channel [START_REF] Tu | OFDMA for underwater acoustic communications[END_REF] motion result in a highly time varying environment. These effects also cause changes in the position of the reflection point from the surface of the sea resulting in Doppler spreading.

B. Channel Model

Unlike the radio channel, there is no standardized statistical characterization of underwater acoustic channels. Through experimental results it has been proved that some channels may be characterized as deterministic while others have shown to exhibit Rice, Rayleigh or even K-distributed fading [START_REF] Chitre | A High-Frequency Warm Shallow Water Acoustic Communications Channel Model and Measurements[END_REF], [START_REF] Yang | High-Frequency Channel Characterization for M-ary Frequency-Shift-Keying Underwater Acoustic Communications[END_REF].

In general, the frequency response of the p th propagation path between transmitter and the receiver is represented by

H p (f ) = Γ p γ(l p , f )
Where Γ p is the cumulative reflection coefficient along the p th path [START_REF] Stojanovic | Underwater acoustic communication channels: Propagation models and statistical characterization[END_REF]. The overall impulse response of the channel is given as

h(t) = p h p (t)δ(t -τ p )
where τ p are the path delays and h p is the inverse fourier transform of H p .

As UUVs travel with a typical speed of a few meters per second, the random channel for simulation purposes in this paper has been modelled as a three tap channel [START_REF] Tu | OFDMA for underwater acoustic communications[END_REF] with channel taps generated by using a Rician distribution with expected values of tap gains as shown in fig. 1 andK factor in the range of 2-6 dB [START_REF] Radosevic | Statistical characterization and capacity of shallow water acoustic channels[END_REF].

III. NETWORK MODEL

Consider a network of unmanned vehicles deployed in a shallow underwater environment. The network primarily 2. The first type is a swimming survey vehicle that performs a wide scale mapping of a particular area to carry out a specific task such as detection of an object. In addition to the swimming vehicle, the network consists of a set of crawling vehicles which accompany the swimming vehicle in performing detailed inspection within its allocated zone. The swimming vehicle can be interfaced with a surface device which can be communicate wirelessly to a base station on the ground. These type of networks are used for scientific missions such as the ones meant for oceanographic sampling [START_REF] Freitag | Acoustic communication in very shallow water: Results from the 1999 AUV Fest[END_REF] in an autonomous fashion.

The network considered in this paper consists of a swimming vehicle with K crawling vehicles assisting it. The purpose of using multiple vehicles is to ensure faster task execution as compared to the case of a single vehicle. The area to be surveyed is divided into several small adjacent zones by the swimming vehicle and a crawling vehicle is assigned per zone. Thus, the zone allocation ensures that no two zones for crawling vehicles overlap each other. The data gathered by each crawling vehicle from within its zone is sent to its corresponding swimming vehicle. This ensures continuous progress monitoring as well as efficient terrain evaluation and target detection.

The crawling vehicles communicate with their corresponding swimming vehicles simultaneously and use DS-CDMA with a M-chip spreading code (explained later). For the swimming vehicle, the n th data transmitted by its k th crawling vehicle is given by

x n,k (t) = Re s n,k M -1 i=0 e jα i,k c(t -iT c )
where s n,k is the n th BPSK (Binary Phase Shift Keying) modulated data bit of the k th crawling vehicle, e jα i,k represents the i th component of the k th crawling vehicle's spreading sequence, c(t) is the chip shape and T c is the chip duration. Thus, the n th data transmitted by all the crawling vehicles working in association with the swimming vehicle in its zone is given by

x n (t) = Re K-1 k=0 s n,k M -1 i=0 e jα i,k c(t -iT c )
From the adopted channel model, the corresponding received signal at the swimming vehicle can be expressed as

y n (t) = Re L-1 p=0 h p K-1 k=0 s n,k M -1 i=0 e j(2πfct+φp+α i,k ) × c(t -iT c -τ p ) + η(t)
where L represents the total number of paths, f c the carrier frequency, φ p the phase offset of the p th path and η(t) the additive white Gaussian noise (AWGN) with power spectral density (PSD) N 0 . The path amplitude are assumed to be constant within one transmitted block i.e. h p (t) ≈ h p [START_REF] Mason | An OFDM Design for Underwater Acoustic Channels with Doppler Spread[END_REF].

IV. DIRECT SEQUENCE CODE DIVISION MULTIPLE ACCESS (DS-CDMA)

Code-Division is a multiple access technique in which the time and bandwidth are simultaneously utilized by various users. The transmitted signals of the users are modulate by orthogonal or semi-orthogonal spreading codes. The code structure is utilized by the receiver to separate out different users. One of the biggest advantages of spread spectrum technique is that dynamic co-ordination expected from the users with respect to time or frequency is little as they can be distinguished by the properties of the spreading codes alone. In case of TDMA and FDMA, the time and bandwidth is divided into different pieces which lays hard limits on the number of users which can simultaneously access the system. The same is true for CDMA systems using orthogonal codes. On the other hand, if semi-orthogonal codes are used, the interference lays a limit on the number of users that can be used.

A. Spreading Codes

In direct sequence (DS) CDMA, each user's information is generally spread by using orthogonal pseudo noise (PN) codes such as Walsh-Hadamard (WH) which helps distinguish between different users. In this report, the carrier interferometry (CI) codes have been used and compared with WH codes. A detailed description of the two spreading codes is given below:

1) Walsh-Hadamard (WH) codes: WH codes originate from the Walsh Hadamard transform which is well known in digital signal processing (DSP) literature [START_REF] Rappaport | Wireless Communications[END_REF] which is used to decompose a discrete-time signal onto a basis of 1's and -1's. The transform is an M ×M square matrix with individual elements given by

T W H (a, b) = (-1) f (a,b)
where T W H (a, b) is the element in the a th row and b th column and f (a, b) is a function chosen such that individual rows of the transform matrix T W H constitute independent vectors. The matrix can also be recursively generated by using the following procedure

T W H,M = T W H, M 2 T W H, M 2 T W H, M 2 -T W H, M 2
where T W H,2 is given by

T W H,2 = 1 1 1 -1
Each normalized row of this matrix constitutes a set of Walsh-Hadamard codes which is used for spreading data of individual users. As can be observed, WH transform matrix dimensions can only be powers of 2 and thus corresponding codes of other sizes are not possible.

2) Carrier Interferometry (CI) codes: The CI codes are based on the Discrete Fourier Transform (DFT) whose M ×M matrix elements are given by

T CI (a, b) = e j 2π M a,b
where T CI (a, b) is the element in the a th row and the b th column. Previously CI codes when applied to MC-CDMA systems [START_REF] Nassar | Introduction of carrier interference to spread spectrum multiple access[END_REF], [START_REF] Wu Zhiqiang | High-performance, high-capacity MC-CDMA via carrier interferometry[END_REF], [START_REF] Hui | MC-CDMA via carrier interferometry Codes in an Underwater Acoustic Channel[END_REF] produced better BER results as compared to those with Walsh Hadamard codes due to their phase offset independent SIR.

V. TRANSCEIVER DESIGN

The transmitter system for a DS-CDMA system is as shown in fig. 3. wherein each data symbol is spread by using a specific spreading sequence before transmission. The RAKE receiver coherently combines the components received on multiple paths by making use of the autocorrelation properties of the spreading codes. Independent multipath components that are separated by more than a chip time can be resolved. Thus, the RAKE receiver can be thought of as another form of diversity combining. The signal corresponding to the k th crawling vehicle on the j th finger of the RAKE receiver shown in fig. 4 is given by ŷk (t) =

(s n,k ) = +1 if Ŷ > 0 -1, otherwise VI. MULTIMEDIA TRANSMISSION IN UNDERWATER ACOUSTIC CHANNEL
With advancement in hardware equipments associated with Underwater Vehicles, the demand for deployment of underwater acoustic multimedia networks is on the rise. Such networks will play an important in accomplishing a variety of tasks related to monitoring and survey. A few examples of such activities are multimedia tactical surveillance of coastal areas to protect the national borders and enable safe passage of trade ships, studying marine life and habitat, exploration under the sea for various minerals and other aquatic resources, picture and video acquisition and classification and prevention of disasters. Due to the difference in requirements of the applications involved, it is necessary to provide different service support to traffic that is loss tolerant and delay tolerant as well as loss sensitive and loss tolerant.

A frequency dependent transmission loss, multipath effect, high propagation delay, battery and other resource constraints make the delivery of multimedia content a challenging task. Multirate transmission scheme are therefore important and variable spreading length (VSL) and multi-code schemes are two promising schemes to realize such transmissions.

A. Variable spreading length scheme

In variable spreading length scheme, a single data rate service is assigned a single spreading sequence. However, the length of the spreading sequence is different for different users. The processing gains are set to match the data rates of different services. Multipath interference becomes extremely significant as data rate increases due to the low processing gain assigned to the user.

B. Multi-code scheme

In the multi-code scheme, the high data-rate stream is spread into multiple low data rate sub-streams. Within each sub-stream, the data is spread by using a random spreading sequence and then transmitted in parallel multi-code channels. The spreading codes assigned to different sub-streams are orthogonal to each other. This is done to avoid interference between different channels when the time delay is zero. The high data rate is achieved by increasing the number of parallel sub-channels. While doing so, the processing gain is maintained constant. Similar to the VSL scheme, the multipath interference (MPI) becomes prominent at higher data rate due to the increased number of parallel sub-channels.

Multi-code schemes have a better performance in multipath environments as compared to variable spreading length schemes. Therefore, we use the multi-code scheme along with DS-CDMA/CI and study the performance of the system in an underwater environment.

VII. SIMULATION RESULTS AND DISCUSSION

The simulation has been performed in three stages. In the first stage, the performance of DS-CDMA/CI has been compared to that of DS-CDMA/WH on the basis of BER plot. Next, bit error rate plot for varying the number of crawling vehicles has been shown followed by the results obtained by transmitting an image in an underwater acoustic communication system using DS-CDMA/CI. Finally, a multicode DS-CDMA/CI transmission scheme has been used to study the use of carrier interferometry codes in multirate transmission systems. The block diagram of the general system is as given in fig. 5.

A. Comparison of CI codes with WH codes

For comparison between CI codes and WH codes, the system model described in Section III is used with the underwater channel model and a simulation is performed in MATLAB. We consider a system consisting of a single swimming vehicle which is accompanied by 4 crawling vehicles. The spreading length is taken as 32 for both CI codes and WH codes for making the comparison fair. The chip duration (T c ) is taken to be 2.5ms. The corresponding bit error rate (BER) plot is as shown in fig. 6. As can be seen from the plot, the performance of the system is far better with CI codes which demonstrate nearly 8 dB gain as compared to that with WH codes for a 6 crawling vehicle system. This is particularly important as underwater acoustic systems are expected to be power efficient. The reason for this high improvement in performance is the that the SIR in case of CI codes is independent of the phase offsets between different paths. On the other hand, in case of WH spreading codes, the SIR is dependent on the squared cosine of the phase offset which causes degradation in the BER performance [START_REF] Wu Zhiqiang | Novel orthogonal codes for DS-CDMA with improved crosscorrelation characteristics in multipath fading channels[END_REF]. As a result, CI is clearly an excellent candidate for use in such underwater vehicular networks.

B. UUV network operation using DS-CDMA/CI

For this simulation, the value of K is varied from 9 through 11. The spreading length is taken as 32 and the chip duration is 2.5ms. As expected, the BER shows degradation as the number of UUVs increase as is indicated by fig. 7. The system model and simulation parameters are the same as above. Following this, the system's performance is studied for image transmission. This simulation was intended to study image transmission at different SNR values. A sample image is transmitted through an underwater channel and is obtained at the receiver. DS-CDMA is used with CI spreading codes of the same spreading length and chip duration as taken previously. The number of crawling vehicles is kept fixed at 8. The results are as shown in fig. 8. As expected, the image quality is slightly diminished at lower SNR values and improves at higher SNR. It can be seen that at a low SNR value, the system performs satisfactorily.

C. Multirate transmission using multi-code DS-CDMA/CI

For this simulation, we consider an underwater vehicle deployed in an underwater environment. Such systems have been studied in [START_REF] Lin | Underwater Acoustic Multimedia Communication Based on MIMO-OFDM[END_REF] wherein the hardware design possesses an ability to transmit sensor data, image and audio. The transmission scheme used is Multi-code DS-CDMA/CI and the performance of the system has been evaluated. The number of bits transmitted is 10 6 and the spreading length is taken to be 32. The data rate of the higher data rate transmission is taken to be twice that of the low data rate transmission. In multimedia transmission schemes for underwater environment, the acceptable transmission BERs for data, audio and image packets are typically taken as 10 -6 , 10 -3 , 10 -4 respectively. The bit error rate (BER) plot for the above system is as shown in fig. 9. Thus, the multi-code DS-CDMA/CI has a satisfactory performance for such multirate communication in underwater acoustic channels.

VIII. CONCLUSION

This report studies the performance improvement in underwater acoustic communications in a multiple UUV network using DS-CDMA with carrier interferometry (CI) spreading codes. The performance has been compared to the case wherein the vehicles use Walsh-Hadamard (WH) spreading codes as the resultant crosscorrelation is independent of the phase offset between different paths. It has been found that use of CI spreading codes in place of WH codes greatly improve the system performance. Also, the performance of CI is thoroughly studied by varying the number of crawling vehicles in the network as well as by image transmission using DS-CDMA/CI. Furthermore, the performance of multi-code DS-CDMA/CI has also been shown and a satisfactory system performance has been achieved. Thus, CI can be considered as an excellent candidate for underwater UUV networks which use DS-CDMA technique for communication.

IX. FUTURE WORK

The network model used in this report considers multiple crawling vehicles operating under a single swimming vehicle. However, in a general case, a number of such groups might be operating in the vicinity of each other. Simulations for such cases in which the value of the number of swimming vehicles is greater than one will be a good future work.

The multi-rate transmission scheme in this report has considered the high data rate user to transmit at twice the rate of the low data rate user. Although, this simulation is satisfactory for the current work, it can be further used to study the multi-rate services in which the high data rate is an integral multiple of the low data rate-the integer being a positive number other than two.

While operating in a specific region, it is likely that the UUV network will come across another network operating in the same area. As a result, the UUV network might end up obstructing this network causing performance degradation in them. Therefore, designing a suitable scheme for minimizing this disturbance is necessary and this is another possibility for future work.

Fig. 2 .

 2 Fig. 2. Multiple UUV communication network

Ts 0 T

 0 k (t -τ j )y(t) dt these individual signals on RAKE receiver fingers are then combined together to generate a variable ( Ŷ ) which then is used to perform a hard decision. The diversity combining

Fig. 3 .Fig. 4 .

 34 Fig. 3. Multiple UUV communication network

Fig. 5 .Fig. 7 .

 57 Fig. 5. Block diagram

Fig. 8 .

 8 Fig. 8. (a) Transmitted image (b) Received image at SNR=15dB (c) Received image at SNR=20dB (d) Received image at SNR=25dB

For more information on WLANs, reader can refer to[START_REF] Nayak | Multi-user downlink with single-user uplink can starve TCP[END_REF],[START_REF] Nayak | AP-side WLAN Analytics[END_REF],[START_REF] Nayak | Modeling Multi-User WLANs Under Closed-Loop Traffic[END_REF],[START_REF] Nayak | Virtual Speed Test: an AP Tool for Passive Analysis of Wireless LANs[END_REF],[START_REF] Nayak | Performance Evaluation of MU-MIMO WLANs Under the Impact of Traffic Dynamics[END_REF],[START_REF] Nayak | Modeling Multi-User WLANs Under Closed-Loop Traffic[END_REF],[START_REF] Nayak | uScope: a Tool for Network Managers to Validate Delay-Based SLAs[END_REF],[START_REF] Nayak | A Novel Compact Dual-Band Antenna Design for WLAN Applications[END_REF],[START_REF] Nayak | A Novel Triple Band Antenna Design for Bluetooth, WLAN and WiMAX Applications[END_REF],[START_REF] Nayak | Design and Analysis of a Multiband Fractal Antenna for Applications in Cognitive Radio Technologies[END_REF],[START_REF] Nayak | Underwater Communications for UUV System[END_REF],[START_REF] Nayak | Passive AP-side Tool for WLAN Analysis: CAP Poster Abstract[END_REF],[START_REF] Nayak | Modeling and Performance Evaluation of MU-MIMO Under the Impact of Closed Loop Traffic Dynamics: CAP Poster Abstract[END_REF]

For further reading on compact antenna design, readers can refer to[START_REF] Nayak | Ultrawideband (UWB) antenna design for cognitive radio[END_REF],[START_REF] Nayak | Multiband Antenna Designs for Wireless Communication Systems[END_REF],[START_REF] Nayak | Multiband fractal antenna design for Cognitive radio applications[END_REF],[START_REF] Nayak | Compact dual-band antenna for WLAN applications[END_REF],[START_REF] Endluri | A Low Cost Dual Band Antenna for Bluetooth, 2.3 GHz WiMAX and 2.4/5.2/5.8 GHz WLAN[END_REF],[START_REF] Nayak | Novel Compact Dual-Band Antenna Design for WLAN Applications[END_REF],[START_REF] Nayak | A Study of Ultrawideband (UWB) Antenna Design for Cognitive Radio Applications[END_REF].


Peter K Kang

Marco Dentz

Tanguy Le Borgne

Ruben Juanes

Spatial Markov Model of Anomalous Transport Through Random Lattice Networks

Keywords: numbers: 05.60.Cd, 05.40.Fb, 05.10.Gg, 92.40.Kf

Flow through lattice networks with quenched disorder exhibits strong correlation in the velocity field, even if the link transmissivities are uncorrelated. This feature, which is a consequence of the divergence-free constraint, induces anomalous transport of passive particles carried by the flow. We propose a Lagrangian statistical model that takes the form of a continuous time random walk (CTRW) with correlated velocities derived from a genuinely multidimensional Markov process in space. The model captures the anomalous (non-Fickian) longitudinal and transverse spreading, and the tail of the mean first passage time observed in the Monte Carlo simulations of particle transport. We show that reproducing these fundamental aspects of transport in disordered systems requires honoring the correlation in the Lagrangian velocity.

Anomalous transport, understood as the nonlinear scaling with time of the mean square displacement of transported particles, is observed in many physical processes, including contaminant transport through porous and fractured geologic media [1], animal foraging patterns [2], freely diffusing molecules in tissue [3], tracer diffusion in suspensions of swimming micro-organisms [4], and biased transport in complex networks [5].

Anomalous transport often leads to a broad-ranged particle distribution density, both in space and time [6][7][8]. Understanding the origin of the slow-decaying tails in probability density is essential, because they determine the likelihood of high-impact, "low-probability" events and therefore exert a dominant control over the predictability of a system [9]. This becomes especially important when human health is at risk, such as in epidemic spreading through transportation systems [START_REF] Balcan | Proc. Natl. Acad. Sci. USA[END_REF] or radionuclide transport in the subsurface [START_REF] Molinero | [END_REF].

Past studies have shown that high variability in the flow properties leads to anomalous transport [1,7]. Depending on the nature of the underlying disorder distribution anomalous behavior can be transient or persist to asymptotic scales [12,13]. The continuous time random walk (CTRW) formalism [14,15] offers an attractive framework to understand and model anomalous transport through disordered media and networks [1,5,16]. The CTRW model is intrinsically an annealed model because the disorder configuration changes at each random-walk step. A particle that returns to the same position experiences different velocity properties. The validity of the CTRW approach for average transport in quenched random environments has been studied for purely diffusive transport [e.g., 7] and biased diffusion [e.g., 9,[17][18][19]. Most studies that employ the CTRW approach assume that transition times associated with particle displacements are independent of each other, therefore neglecting velocity correlation between successive jumps [20]. Indeed, a recent study of transport on a lattice network has shown that CTRW with independent transition times emerges as an exact macroscopic transport model when velocities are uncorrelated [9].

However, detailed analysis of particle transport simulations demonstrates conclusively that particle velocities in massconservative flow fields exhibit correlation along their spatial trajectory [17,21,22]. Mass conservation induces correlation in the Eulerian velocity field because fluxes must satisfy the divergence-free constraint at each intersection. This, in turn, induces correlation in the velocity sequence along a particle trajectory. To take into account velocity correlation, Lagrangian models based on temporal [22,23] and spatial [17,21] Markovian processes have recently been proposed. These models successfully capture many important aspects of the Lagrangian velocity statistics and the particle transport behavior. In particular, the study of Le Borgne et al. [17] shows that introducing correlation in the Lagrangian velocity through a Markov process in space yields an accurate representation of the first and second moments of the particle density. The model is restricted, however, to particle trajectories projected onto the direction of the mean flow, and the study leaves open the question of whether spatial Markov processes can describe multidimensional features of transport.

Here, we investigate average transport in divergence-free flow through a quenched random lattice from the CTRW point of view. We introduce a multidimensional spatial Markov model for particle velocity, and confirm that the model exhibits excellent agreement with Monte Carlo simulations. We show that accounting for the spatial correlation in the Lagrangian velocity is essential to capture the fundamental macroscopic transport behavior.

Random Lattice Network. We consider a lattice network consisting of two sets of parallel, equidistant links oriented at an angle of ±α with the x-axis. The distance between nodes is l [Fig. 1(a)]. Flow through the network is modeled by assuming Darcy's law [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF] for the fluid flux u ij between nodes i and j, u ij = -k ij (Φ j -Φ i )/l, where Φ i and Φ j are the flow potentials, and k ij > 0 is the conductivity of the link between the two nodes. Imposing mass conservation at each node i, j u ij = 0, leads to a linear system of equations, which is solved for the flow potentials at the nodes. A link from node i to j is incoming for u ij < 0 and outgoing for u ij > 0. We denote by e ij the unit vector in the direction of the link connecting i and j. A realization of the random lattice network is generated by assigning independent and identically distributed random conductivities to each link. Therefore, the k values in different links are uncorrelated. The set of all realizations of the quenched random network generated in this way form a statistical ensemble that is stationary and ergodic.

We study a simple flow setting of mean flow in the positive x-direction, by imposing a no-flow condition at the top and bottom boundaries of the network, and fixed values of the potential at the left (Φ = 1) and right (Φ = 0) boundaries.

Once the fluxes at the links are known, we simulate transport of a passive tracer by particle tracking. We neglect diffusion along links, and thus particles are advected with the flow velocity between nodes. We assume complete mixing at the nodes. Thus, the link through which the particle exits a node is chosen randomly with flux-weighted probability. The Langevin equations describing particle movements in space and time are The system of discrete Langevin equations (1) describes coarse-grained particle transport for a single realization of the quenched random lattice. Particle velocities and thus transition times depend on the particle position. The particle position at time t is x(t) = x nt , where n t denotes the number of steps needed to reach time t. The mean particle density is P (x, t) = δ(xx nt ) , where the angular brackets denote both the noise average over all particles in one realization and the ensemble average over all network realizations. We solve transport in a single disorder realization by particle tracking based on Eq. ( 1) with the initial conditions x 0 = 0 and t 0 = 0 [Fig. 1(b)]. From this, we obtain the mean particle density P (x, t) by ensemble averaging.

x n+1 = x n + l v(x n ) |v(x n )| , t n+1 = t n + l |v(x n )| . ( 1
To develop a transport model for the average particle density P (x, t), we study average particle movements from a CTRW point of view. This could be done, for example, by interpreting first-passage time distributions in the CTRW framework and inferring an optimal distribution of transition times [20]. Here we follow a different rationale and analyze the ensemble statistics of the Lagrangian velocities because the CTRW model is based on the assumption that particle velocities sampled at given spatial positions along an average trajectory form a Markov process.

Spatial Markov Property. To characterize average particle movement from a CTRW point of view, we study the ensemble statistics of the series of Lagrangian velocities experienced by particles along individual trajectories. We consider the transition probability density to encounter a velocity v after n + m steps given that the particle velocity was v ′ after n steps, which in the variables (ν, θ) reads

r m (ν, θ|ν ′ , θ ′ ) = δ [ν -ν(x n+m )] δ θ,θ(xn+m) ν(xn)=ν ′ ,θ(xn)=θ ′ . (2)
We study the statistical properties of the Lagrangian velocity v(x n ) by particle tracking simulations in 10 3 realizations of an ensemble of random lattices characterized by a lognormal k distribution with variance σ 2 ln k = 5. The use of a lognormal distribution is motivated by measurements of conductivity in many natural media [START_REF] Sanchez-Vila | [END_REF]. The lattice size is 500×500 nodes and, in each realization, we release 10 3 particles at the origin. To evaluate the transition probability numerically, the particle velocity ν is discretized into classes, ν ∈ N j=1 (ν j , ν j+1 ). To emphasize the role of low velocities, velocity classes are defined on a near-logarithmic scale. We define the transition probability matrix

T m (i, θ|j, θ ′ ) = νi+1 νi dν νj+1 νj dν ′ r m (ν, θ|ν ′ , θ ′ ). (3)
The aggregate transition matrix T m (i|j) = θ,θ ′ T m (i, θ|j, θ ′ ) shown in Fig. 2a for m = 1, clearly indicates that particle velocities are correlated. The relatively large probabilities in the upper-left and lower-right corners of the transition matrix reflect flow reversal.

The series of Lagrangian velocities v(x n ) ≡ v n along particle trajectories can be approximated as a Markov process, if the transition matrix satisfies the Chapman-Kolmogorov equation [e.g., 26], which in matrix form reads Specifically, for a Markov process, the m-step transition matrix T m is equal to the m-fold product of the 1-step transition matrix T 1 with itself as T m = T m . Figure 2c shows the transition probabilities for m = 5 steps conditional to a low (j = 5) and high (j = 90) velocity class given by T 5 , which is obtained by direct Monte Carlo simulations, and under the Markov assumption from T 5 . The Markov model predicts accurately the transition probabilities, as well as the return probability for any number of steps [Fig. 2b, inset]. Our analysis suggests that the Markov model captures the Lagrangian velocity statistics accurately. We repeated the analysis for truncated power-law and the absolute value of Cauchy distributions of conductivity, and found that the Markovianity assumption holds for these conductivity distributions too. Therefore, a CTRW characterized by a one step correlation in velocity is a good approximation for describing average transport.

T n (i, θ|j, θ ′ ) = i ′ ,θ ′′ T n-m (i, θ|i ′ , θ ′′ )T m (i ′ , θ ′′ |j, θ ′ ). (4)
(v n ) ( v n ) +1 (
Continuous Time Random Walk Model. Particle movements in the random lattice can, on average, be described by the following system of Langevin equations

x n+1 = x n + l v n |v n | , t n+1 = t n + l |v n | . ( 5
)
We have already shown that the series of Lagrangian velocities {v n } ∞ n=0 is well approximated by a Markov process and thus fully characterized by the one-point density p(v) = δ(vv n ) and the one-step transition probability density

r 1 (v|v ′ ) = δ(v -v n+1 ) | vn=v ′ . ( 6
)
The particle density can be written as

P (x, t) = dv δ(x -x nt )δ(v -v nt ) , (7)
in which n t = max(n|t n ≤ t), x is the position of the node at which the particle is at time t, and v is the velocity by which the particle emanates from this node. Equation ( 7) can be recast as

P (x, t) = dv t t-l/|v| dt ′ R(x, v, t ′ ), (8a)
in which we defined

R(x, v, t ′ ) = ∞ n=0 δ(x -x n )δ(v -v n )δ(t ′ -t n ) . (8b)
The latter satisfies the Kolmogorov type equation

R(x, v, t) = δ(x)p(v)δ(t) + dv ′ r 1 (v|v ′ ) × dx ′ δ(x -x ′ -lv ′ /|v ′ |)R(x ′ , v ′ , t -l/|v ′ |). (8c)
For independent successive velocities, i.e., r 1 (v|v ′ ) = p(v), one recovers the CTRW model [e.g., 14]

P (x, t) = t 0 dt ′ R(x, t ′ ) ∞ t-t ′ dτ dxψ(x, τ ), (9a)
where R(x, t) satisfies

R(x, t) = δ(x)δ(t) + dx ′ t 0 dt ′ R(x ′ , t ′ )ψ(x -x ′ , t -t ′ ) (9b)
and the joint transition length and time density is given by

ψ(x, t) = dv ′ p(v ′ )δ(x -lv ′ /|v ′ |)δ(t -l/|v ′ |). (9c)
In the following, we refer to system (8) as correlated CTRW because subsequent particle velocities are correlated, and to model (9) as uncorrelated CTRW because subsequent particle velocities are uncorrelated.

Average Transport Behavior. The average transport behavior is studied in terms of the spatial particle density P (x, t), its mean square displacements in longitudinal and transverse directions and the distribution of the first passage time, t f (x), at a control plane at a distance x from the inlet. We compare the results obtained from direct Monte Carlo simulations to correlated CTRW and uncorrelated CTRW. Correlated CTRW is parametrized by the one-step transition matrix T 1 determined from numerical Monte Carlo simulations. Uncorrelated CTRW is parametrized by the Lagrangian velocity distribution p(v), which is obtained from Monte Carlo simulations as well.

The particle distribution is non-Gaussian and characterized by a sharp leading edge and an elongated tail [Fig. 3]. The non-Gaussian features persist even after the center of mass has travelled a distance of about 100 links in the direction of the mean flow. Correlated CTRW captures the shape of the particle plume with remarkable accuracy, including its leading edge, peak, transverse spread, and low-probability tail near the origin. Ignoring the correlated structure of the Lagrangian velocity leads to predictions of longitudinal and transverse spreading that deviate from the direct Monte Carlo simulation [Fig. 3, insets].

Figure 4a shows the time evolution of the longitudinal and transverse spreading. The Monte Carlo simulation shows that the longitudinal mean square displacement (MSD) with respect to the center of mass evolves faster than linear with time (slope of 1.33). Both the scaling and the magnitude of the longitudinal spreading are captured accurately by correlated CTRW. The model also reproduces accurately the magnitude and time scaling of the transverse MSD. The uncorrelated model underpredicts the magnitude of longitudinal spreading.

Nonlocal theories of transport, including CTRW, are often invoked to explain the empirical observation that the first passage time (FPT) distribution is broad-ranged [1]. Early arrival and slow decay of the FPT is also observed in our model system, even when the conductivity distribution is lognormal and has zero spatial correlation [Fig. 4(b)]. The cumulative FPT distribution from the Monte Carlo simulation exhibits a significantly slower decay than uncorrelated CTRW. This behavior is accurately captured by correlated CTRW, suggesting that the velocity correlation along particle trajectories is responsible for the emergence of the observed asymptotic behavior.

In conclusion, we have shown that the divergence-free condition arising from mass conservation is the source of strong and nontrivial correlation in the Lagrangian velocity, even when the underlying conductivity field is completely uncorrelated. Accounting for such correlation in the velocity is important to obtain quantitative agreement for the mean particle density and the FPT distribution. Here, we have proposed and validated a spatial Markov model of transport on a lattice network that explicitly captures the multidimensional effects associated with changes in direction along the particle trajectory. This study opens the door to understanding the interplay between two sources of velocity correlation: the divergence-free condition and the spatial correlation in the permeability field. Finally, we suspect that correlation in the Lagrangian velocity exerts an even more dominant control over mixing (understood as the decay of the variance of the particle density [START_REF] Le Borgne | [END_REF][28][29]) than it does on spreading. This remains an exciting open question.

We gratefully acknowledge funding for this research, provided by the DOE Office of Science Graduate Fellowship Program (to PKK), the Spanish Ministry of Science and Innovation through the project HEART (CGL2010-18450) (to MD), the European

) 5 FIG. 1 .

 51 FIG. 1. (a) Schematic of the lattice network considered here, with two sets of links with orientation ±α = ±π/4 and spacing l = 1. (b) Particle distribution at nodes (represented by circles of different sizes) at t = 30 for a single realization after injection at the origin at t = 0.

FIG. 2 .

 2 FIG. 2. (a) Aggregate transition matrix for N = 100 velocity classes distributed with logarithmic scale. (b) Transition probabilities after m = 5 steps from direct Monte Carlo computation (blue solid line) and calculated from the Markov assumption (green symbols). Shown are probability densities for two initial velocity classes: a low velocity class (j = 5, •), and a high velocity class (j = 90, * ). Inset: probability of returning to the same initial velocity class as a function of the number of steps for a high initial velocity (class j = 90).

FIG. 3 .

 3 FIG. 3. Contour plot of the mean particle density at t = 5 × 10 2 , computed from direct Monte Carlo simulation (blue solid line), correlated CTRW model (green solid line), and uncorrelated CTRW model (red solid line).

FIG. 4 .

 4 FIG. 4. (a) Time evolution of the longitudinal MSD. Inset: Transverse MSD. (b) Cumulative FPT distribution.


Foreword by Denis Noble

During most of the twentieth century experimental and theoretical biologists lived separate lives. As the authors of this book express it, "there was a belief that experimental and theoretical thinking could be decoupled." This was a strange divorce. No other science has experienced such a separation. It is inconceivable that physical experiments could be done without extensive mathematical theory being used to give quantitative and conceptual expression to the ideas that motivate the questions that experimentalists try to answer. It would be impossible for the physicists at the large hadron collider, for example, to search for what we call the Higgs boson without the theoretical background that can make sense of what the Higgs boson could be. The gigantic masses of data that come out of such experimentation would be an un-interpretable mass without the theory. Similarly, modern cosmology and the interpretation of the huge amounts of data obtained through new forms of telescopes would be inconceivable without the theoretical structure provided by Einstein's general theory of relativity. The phenomenon of gravitational lensing, for example, would be impossible to understand or even to discover. The physics of the smallest scales of the universe would also be impossible to manage without the theoretical structure of quantum mechanics.

So, how did experimental biology apparently manage for so many years without such theoretical structures? Actually, it didn't. The divorce was only apparent.

First, there was a general theoretical structure provided by evolutionary biology. Very little in biology makes much sense without the theory of evolution. But this theory does not make specific predictions in the way in which the Higgs boson or gravitational lensing were predicted for physicists. The idea of evolution is more that of a general framework within which biology is interpreted.

Second, there was theory in biology. In fact there were many theories, and in many different forms. Moreover, these theories were used by experimental biologists. They were the ideas in the minds of experimental biologists. No science can be done without theoretical constructs. The so-called Central Dogma of Molecular Biology, for example, was an expression of the background of ideas that were circulating during the early heydays of molecular biology: that causation was one way (genes to phenotypes), and that inheritance was entirely attributable to DNA, by which an organism could be completely defined. This was a theory, except that it was not formulated as such. It was presented as fact, a fait accompli. Meanwhile the pages of journals of theoretical and mathematical biology continued to be filled with fascinating and difficult papers to which experimentalists, by and large, paid little or no attention.

We can call the theories that experimentalists had in mind implicit theories. Often they were not even recognised as theory. When Richard Dawkins wrote his persuasive book The Selfish Gene in 1976 he was not only giving expression to many of these implicit theories, he also misinterpreted them through failing to understand the role of metaphor in biology. Indeed, he originally stated "that was no metaphor"! As Poincaré pointed out in his lovely book Science and Hypothesis (La science et l'hypothèse) the worst mistakes in science are made by those who proudly proclaim that they are not philosophers, as though philosophy had already completed its task and had been completely replaced by empirical science. The truth is very different. The advance of science itself creates new philosophical questions. Those who tackle such questions are philosophers, even if they do not acknowledge that name. That is particularly true of the kind of theory that could be described as meta-theory: the creation of the framework within which new theory can be developed. I see creating that framework as one of the challenges to which this book responds.

Just as physicists would not know what to do with the gigantic data pouring out of their colliders and telescopes without a structure of interpretative theory, biology has hit up against exactly the same problem. We also are now generating gigantic amounts of genomic, proteomic, metabolomic and physiomic data. We are swimming in data. The problem is that the theoretical structures within which to interpret it are underdeveloped or have been ignored and forgotten. The cracks are appearing everywhere. Even the central theory of biology, evolution, is undergoing reassessment in the light of discoveries showing that what the modern synthesis said was impossible, such as the inheritance of acquired characters, does in fact occur. There is an essential incompleteness in biological theory that calls out to be filled.

That brings me to the question how to characterise this book. It is ambitious. It aims at nothing less than filling that gap. It openly aims at bringing the rigour of theory in physics to bear on the role of theory in biology. It is a highly welcome challenge to theorists and experimentalists alike. My belief is that, as we progressively make sense of the masses of experimental data we will find ourselves developing the conceptual foundations of biology in rigorous mathematical forms. One day (who knows when?), biology will become more like physics in this respect: theory and experimental work will be inextricably intertwined.

However, it is important that readers should appreciate that such intertwining does not mean that biology becomes, or could be, reducible to physics. As the au-thors say, even if we wanted such a reduction, to what physics should the reduction occur? Physics is not a static structure from which biologists can, as it were, take things 'off the shelf'. Physics has undergone revolutionary change during the last century or so. There is no sign that we are at the end of this process. Nor would it be safe to assume that, even if it did seem to be true. It seemed true to early and midnineteenth century biologists, such as Jean-Baptiste Lamarck, Claude Bernard, and many others. They could assume, with Laplace, that the fundamental laws of nature were strictly deterministic. Today, we know both that the fundamental laws do not work in that way, and that stochasticity is also important in biology. The lesson of the history of science is that surprises turn up just when we think we have achieved or are approaching completeness.

The claim made in this book is that there is no current theory of biological organisation. The authors also explain the reason for that. It lies in the multi-level nature of biological interactions, with lower level molecular processes just as dependent on higher-level organisation and processes, as they in their turn are dependent on the molecular processes. The error of twentieth century biology was to assume far too readily that causation is one-way. As the authors say, "the molecular level does not accommodate phenomena that occur typically at other levels of organisation." I encountered this insight in 1960 when I was interpreting experimental data on cardiac potassium channels using mathematical modelling to reconstruct heart rhythm. The rhythm simply does not exist at the molecular level. The process occurs only when the molecules are constrained by the whole cardiac cell to be controlled by causation running in the opposite direction: from the cell to the molecular components. This insight is general. Of course, cells form an extremely important level of organisation, without which organisms with tissues, organs and whole-body systems would be impossible. But the other levels are also important in their own ways. Ultimately, even the environment can influence gene expression levels. There is no a priori reason to privilege any one level in causation. This is the principle of biological relativity.

The principle does not mean that the various levels are in any sense equivalent. To quote the authors again: "In no way do we mean to negate that DNA and the molecular cascades that are related to it, play an important role, yet their investigations are far from complete regarding the description of life phenomena." Completeness is the key concept. That is true for biological inheritance as well as for phenotypegenotype relations. New experimental work is revealing that there is much more to inheritance than DNA.

The avoidance of engagement with theoretical work in biology was based largely on the assumption that analysis at the molecular level could be, and was in principle, complete. In contrast, the authros write, "these [molecular] cascades may causally depend on activities at different levels of analysis, which interact with them and also deserve proper insights." Those 'proper insights' must begin by identifying the entities and processes that can be said to exist at the higher levels: "finding ways to constitute theoretically biological objects and objectivise their behaviour." To achieve this we have to distance ourselves from the notion, prevalent in biology today, that the fundamental must be conceptually elementary. As the authors point out, this is not even true in physics. "Moreover, the proper elementary observable doesn't need to be "simple". "Elementary particles" are not conceptually/mathematically simple."

There is therefore a need for a general theory of biological objects and their dynamics. This book is a major step in achieving that aim. It points the way to some of the important principles, such as the principle of symmetry, that must form the basis of such a theory. It also treats biological time in an innovative way, it explores the concept of extended criticality and it introduces the idea of anti-entropy. If these terms are unfamiliar to you, this book will explain them and why they help us to conceptualize the results of experimental biology. They in turn will lead the way by which experimentalists can identify and characterize the new biological objects around which a fully theoretical biology could be constructed.

Oxford University, Denis Noble

June 2013