Patrice Loisel
email: patrice.loisel@inra.fr

Stochastic perturbations and fisheries management

Keywords: Piecewise Deterministic Markov Process (PDMP), optimal control, value function. Recommendations

As most natural resources, fisheries are affected by random disturbances. The evolution of such resources may be modelled by a succession of deterministic process and random perturbations on biomass and/or growth rate at random times. We analyze the impact of the characteristics of the perturbations on the management of natural resources. We highlight the importance of using a dynamic programming approach in order to completely characterize the optimal solution, we also present the properties of the controlled model and give the behavior of the optimal harvest for specific jump kernels.

• In the context of updated biomass, for a centrally disturbed biomass and sufficiently high effort the optimal harvest increases with biomass jump rate • In the context of jointly updated biomass and growth rate, for a centrally disturbed biomass and sufficiently high effort the optimal harvest increases with the biomass jump rate • In the context of jointly updated biomass and growth rate, for sufficiently high effort the optimal harvest decreases with the growth jump rate

Introduction

The evolution of natural resources is most often disturbed by random events. These disturbances occur at times that are not necessarily at regular intervals. Hence, the management of natural resources must take into account the characteristics of these disturbances.

The inclusion of stochastic perturbations in resource management has been the subject of numerous articles in the literature, [START_REF] Clark | Mathematical bioeconomics: the optimal management of renewable resources[END_REF], [START_REF] Pindyck | Uncertainty in the Theory of Renewable Resource Markets[END_REF], [START_REF] Saphores | Harvesting a renewable resource under uncertainty[END_REF], [START_REF] Sethi | Fishery management under multiple uncertainty[END_REF]. Most of these works concerns growth processes subject to perturbation continuously or at regular intervals, whereas, as mentioned above, perturbations occur most often at discrete random times. Our study tries to take into account the random nature of not only the perturbations magnitude of the resource but also the occurrence times of these perturbations.

For systems with random perturbations the state variables are updated at random times and between these random times, the state variables are governed by deterministic processes. The most appropriate framework for the study of such systems seems to be that of the Piecewise Deterministic Markov Process (PDMP) [START_REF] Davis | Markov Models and Optimization[END_REF].

Using this framework, Hanson and Tuckwell [START_REF] Hanson | Population Growth with Randomly Distributed Jumps[END_REF] study the time to extinction of a population with some specific growth function and random perturbation structure. Applications in dynamic population are studied in [START_REF] Ryan | Optimal Harvesting with Exponential Growth in an Environment with Random Disasters and Bonanzas[END_REF], [START_REF] Ryan | Optimal Harvesting of a Logistic Population in an Environment with Stochastic Jumps[END_REF]. Hanson [4] gives a panorama of the models developed in various fields with this framework. We are interested by the optimal management of fisheries in this framework. The goal of this paper is to study the behavior of the control variable with respect to the jump rate of the perturbation process.

In Section 2 we first present the resource growth model with its deterministic and stochastics components and the corresponding PDMP framework. Secondly, we express the controlled problem with updated biomass, we highlight the importance of using a dynamic programming approach in order to completely characterize the optimal solution, we also present the properties of the controlled model and give the application for a specific jump kernel. Finally in Section 4, we consider the case of updated biomass and growth rate.

We assume that the evolution of the biomass in a fishery is mainly governed by a determinist continuous process while it is observed or perturbated only at discrete random times. In absence of update, the evolution of the resource biomass x(t) at time t is governed by a deterministic growth model:

dx(t) dt = G(x(t)) -h(x(t), e(t)),
with initial condition : x(0) = x 0 , with 0 < x 0 < K. The parameter K is the carrying capacity of the studied system and h is the harvest and e(t) is the harvest rate.

We assume that G is a differentiable concave growth function such that G(0

) = 0, G ′ (x) > 0 for x < K, G ′ (x) < 0 for x > K. The most common example is the logistic growth function G(x) = rx(1 -x/K)
with the growth rate r.

We consider a resource submitted to random updates of the biomass Y 1 , Y 2 , .. at random times τ 1 , τ 2 , ....

We assume that updates occur in a Poisson process i.e. that updates occur independently of one another and randomly in time. The distribution of times between successive updates is an exponential distribution with mean 1/λ:

F (x) =1 -e -λx ,
with the constant jump rate λ. For each random time τ i , the biomass is updated:

x(τ + i ) = Y i d ∼ L(.|x(τ i )), at time τ i , for i ≥ 1.
where L is a conditional distribution.

Hence the dynamics of the biomass can be described by a Piecewise Deterministic Markov Process (PDMP) ( [START_REF] Davis | Markov Models and Optimization[END_REF]). The random jump process is described by the jump kernel operator Q. To each function θ the operator Q associates the function Q[θ] which is defined by: Q

[θ](x) = Y θ(Y)dL(Y|x), Q[θ]
is assumed continuous and Lipschitz for all θ ∈ C(Ω).

For example:

Y i = Z i x(τ i ) with Z i d ∼ U (z, z), hence the associated jump kernel Q is: Q[θ](x) = 1 z -z z z θ(zx)dz = 1 (z -z)x zx zx θ(y)dy.

The biomass growth process

In order to describe the biomass growth process, we define the function X(t; x, τ ) at time t. If the biomass was x at time τ , the evolution of X(t; x, τ ) is given by (S x,τ ):

dX(t; x, τ ) dt = g(X(t; x, τ ), e(t)) ≡ G(X(t; x, τ )) -h(X(t; x, τ ), e(t)),
with initial condition: X(τ ; x, τ ) = x.

Knowing these characteristics, denoting τ 0 at x + 0 = x 0 , the biomass growth process {X t (x 0 ) : t ≥ 0} starting with biomass x 0 at initial time, may be expressed, for i ≥ 1:

X t = X(t; x + i-1 , τ i-1 ), τ i-1 < t ≤ τ i , while at time τ i , x + i d ∼ L(.|X τ i ).
The process {X t (x 0 ) : t ≥ 0} starting at x 0 is composed of the successive curves X(t; x 0 , 0), X(t; x + τ 1 , τ 1 ),.., X(t; x + τ i , τ i ), .... This process depends on the successive random updated time τ i and random jump at corresponding τ i .

Remark: by composite construction, for τ i-1 < s < t < τ i , we have: X(t; X(s; x, τ i-1 ), s) = X(t; x, τ i-1 ) and X(t; x, τ ) = X(t -τ ; x, 0).

The control problem

Given a biomass x and an effort e, the instantaneous gain of consumption l(x, e) is determined. Therefore, we assume a regulator maximizing expected discounted gain on an infinite horizon:

J(x 0 , e(.)) = E +∞ 0 l(X t (x 0 ), e(t))e -δt dt , (1)
with the X t (x 0 ) solution obtained with successive systems (S x 0 ,0 ), (S x + τ 1 ,τ 1 ), .... The expectation in Equation (1) is related to the successive random updated time τ i and random jump at corresponding τ i . The effort e is subject to the constraints: 0 ≤ e(t) ≤ e for all t > 0. The instantaneous gain is assumed proportional to the effort: l(x, e) = l 0 (x)e. Thus we consider the function value V defined by: V (x) = max e(.)∈[0,e] J(x, e(.)).

Assuming V ∈ C 1 ([0, K]), we can formally deduce (see Appendix A, with restrictive conditions [START_REF] Dempster | Necessary and sufficient optimality conditions for control of piecewise deterministic processes[END_REF], [START_REF] Davis | Markov Models and Optimization[END_REF] gives mathematical justification) the Bellman Hamilton Jacobi (BHJ) equation:

max e∈[0,e] [V ′ (x)g(x, e) -(δ + λ)V (x) + l(x, e) + λQ[V ](x)] = 0. (2)
The harvest is assumed proportional to the effort e i.e. h(x, e) = h 0 (x)e, the BHJ equation becomes:

max e∈[0,e] [l 0 (x) -h 0 (x)V ′ (x)]e + V ′ (x)G(x) -(δ + λ)V (x) + λQ[V ](x) = 0. (3)
In Equation (3), the effort e depends on the biomass X. This highlights the existence of a critical value x * solution of:

l 0 (x) -h 0 (x)V ′ (x) = 0. (4)
As 0 ≤ e ≤ e, the optimal effort is a feedback control e * (t) = E(X t ) where the function E is defined by:

E(x) =                0, if l 0 (x) -h 0 (x)V ′ (x) < 0, G(x) h 0 (x) , if l 0 (x) -h 0 (x)V ′ (x) = 0, e, if l 0 (x) -h 0 (x)V ′ (x) > 0,
where the function value V is defined by:

[l 0 (x) -h 0 (x)V ′ (x)] + e + V ′ (x)G(x) -(δ + λ)V (x) + λQ[V ](x) = 0. (5)
But Equation ( 4) is not sufficient to characterize the critical values. By using a dynamic programming equation, we obtain a complementary condition based on Euler-Lagrange condition.

Euler-Lagrange condition

The value function V (x 0 ) is the solution to the optimization problem:

V (x 0 ) = J(x 0 , e * (.)) = max e(.)∈[0,e] E +∞ 0 l(X t (x 0
), e(t))e -δt dt , with X t (x 0 ) solution obtained with successive systems (S x 0 ,0 ), (S x + τ 1 ,τ 1 ), ....

Using the strong Markov property, we may express [START_REF] Davis | Markov Models and Optimization[END_REF]:

V (x 0 ) = J(x 0 , e * (.)) =E τ τ 0 l(X(t; x 0 , 0), e * (t))e -δt dt + Q[V ](X(τ ; x 0 , 0))e -δτ .
We consider the first term in the right hand:

E τ τ 0 l(X(t; x 0 , 0), e * (t))e -δt dt = λ +∞ 0 τ 0 l(X(t; x 0 , 0), e * (t))e -δt dte -λτ dτ = +∞ 0 l(X(t; x 0 , 0), e * (t))e -(δ+λ)t dt,
when by inverting integration with respect to t and τ Finally we obtain the dynamic programming equation:

V (x 0 ) = max e(.)∈[0,e] +∞ 0 [l(X(t; x 0 , 0), e(t)) + λQ[V ](X(t; x 0 , 0))]e -(δ+λ)t dt.
To simplify the expressions, we denote X(t, x 0 ) ≡ X(t; x 0 , 0). We have the opportunity, as in the deter-ministic control case, to deduce the expression of the effort e(t) in terms of the biomass:

l(X(t, x 0 ), e(t)) = l 0 (X(t, x 0 ))e(t) = l 0 h 0 (X(t, x 0 )(G(X(t, x 0 )) -Ẋ(t, x 0 ))
and then we obtain the new form of the objective. Thus the optimization problem becomes:

V (x 0 ) = max X(.)∈Cx 0 +∞ 0 l 0 h 0 (X(t, x 0 ))(G(X(t, x 0 )) -Ẋ(t, x 0 )) + λQ[V ](X(t, x 0 )) e -(δ+λ)t dt.
C x 0 being the set of admissible curves:

C x 0 = {X(.) ∈ BC 1 ([0, K]), X(0) = x 0 , G(X(t)) -h 0 (X(t))e ≤ Ẋ(t) ≤ G(X(t))},
with BC 1 stands for the bounded with bounded derivative function defined on the interval [0, K]. We deduce:

Proposition 2.1 Assuming that V and Q[V ] ∈ BC 1 ([0, K]
), a critical value x * is solution of the system of equations:

l 0 (x) -h 0 (x)V ′ (x) = 0, (δ + λ -G ′ (x)) l 0 h 0 (x) = l 0 h 0 ′ (x)G(x) + λ[Q[V ]] ′ (x), (6)
where the value fonction V is solution of Equation [START_REF] Hanson | Population Growth with Randomly Distributed Jumps[END_REF].

Proof: We have an implicit problem of Calculus of Variations. X(.) stands for an interior solution, let L(., .)

the non actualized integrand:

I(X, Ẋ) = l 0 h 0 (X)(G(X) -Ẋ) + λQ[V ](X)
, then X(.) has to satisfy the Euler Lagrange condition:

I X (X(t), Ẋ(t)) = d dt I Ẋ (X(t), Ẋ(t)) -(δ + λ)I Ẋ (X(t), Ẋ(t)).
The Euler Lagrange condition enhances:

l 0 h 0 ′ (X(t))G(X(t) + l 0 h 0 (X(t))G ′ (X(t)) + λ[Q[V ]] ′ (X(t)) =(δ + λ) l 0 h 0 (X(t)).
The differential equation is reduced to an algebraic Equation ( 6).

Let x * (λ) the lower critical value, in order to avoid scaling of function V (V is defined by Equation ( 5) up to a multiplicative constant for x < x * (λ)), by using Equations ( 4) and ( 6) becomes:

(δ -G ′ (x) + λ(1 - [Q[V ]] ′ (x) V ′ (x) )) l 0 h 0 (x) = l 0 h 0 ′ (x)G(x). (7)
Similarly to the standard optimal control problem (without update) the optimal effort e is given by a function E of the biomass X. But this function E depends on the jump rate λ by the intermediate of x * (λ).

The value function and the optimal control

We now analyze the behavior of the value function

V at critical value x * . Let A(x) = l 0 /h 0 (x) -V ′ (x).
Then, for x such that A(x) = 0 we may define

A ′ (x) = l 0 h 0 ′ (x) -V ′′ (x)
and:

V ′′ (x) = l 0 h 0 ′ (x) -A ′ (x), (8)
so we deduce the regularity of the function value with respect to biomass x:

Proposition 2.2
The function value V is continuously twice differentiable and:

V ′′ (x * ) = l 0 h 0 ′ (x * ), (9)
and with respect to jump rate λ:

Proposition 2.
3 For a sufficiently small value of jump rate λ, assuming l 0 (x) = pqx -c and h 0 (x) = qx with price p, catchability q and cost c, the critical value x * is an increasing (respectively decreasing)

function with respect to λ if [Q[V ]] ′ (x * ) -V ′ (x * ) > 0 (respectively < 0)
. Moreover the function value V and V ′ x are continuously differentiable with respect to jump rate λ.

Proofs of Propositions are given in Appendix B. In the following section, we will illustrate for a concrete case, with a specific jump kernel, by a study of the sign of

[Q[V ]] ′ (x * ) -V ′ (x * ). From V ′ x continuously
differentiable with respect to jump rate λ, for a sufficiently small jump rate λ, Equations ( 5) and ( 6) has a unique solution x * so the function E is given by:

E(x) =                0, if x < x * , G(x) h 0 (x) , if x = x * , e, if x > x * ,
and finally:

Proposition 2.4 For a sufficiently small value of jump rate λ, assuming l 0 (x) = pqx -c and h 0 (x) = qx with price p, catchability q and cost c, the value function is not three times differentiable, more precisely at λ = 0:

h 0 (x * )[x 2 V ′′ (x)] ′ + (x * ) = - x * Σ(x * ) e -E(x * ) <0 < h 0 (x * )[x 2 V ′′ (x)] ′ -(x * ) = x * Σ(x * ) E(x * ) , where Σ(x) = -x 2 V ′′ (x) G(x) x ′ + (δ -G ′ (x))(V ′ (x) + V ′′ (x)x) -G ′′ (x)V ′ (x)
x and Σ(x * ) > 0.

We now consider the growth function: G(x) = rx(1 -x/K). We may deduce the behavior of the critical value with respect to growth rate r:

Proposition 2.5 For a sufficiently small value of jump rate λ x and λ r , assuming l 0 (x) = pqx -c and h 0 (x) = qx, the critical value x * (r) is an increasing function with respect to growth rate r.

In Figure 1, we give an example of optimal evolution of the biomass x. The dash line represents the level of the critical value x * .

x

Y i = x(τ i )(1 + Z i ǫ)
where Z i follows the distribution H (symmetric and centered in 0), for a sufficiently small value of jump rate λ x and x * not close to K.

(i) if E[Z] = 0 then the critical value is increasing (respectively decreasing) with respect to jump rate

λ if E[Z] > 0 (respectively < 0).
(ii) if E[Z] = 0 then the critical value is increasing (respectively decreasing) with respect to jump rate

λ if E(x * ) is smaller (respectively larger) than e 2 .
More precisely in the latter case, with the growth function G(x) = rx(1 -x/K), the critical value is increasing (respectively decreasing) with respect to jump rate λ if x * > K(1eq/(2r)) (respectively <).

Proofs are given in Appendix B. The given result for E[Z] = 0 is not surprising: for instance if E[Z] > 0, higher jump rate leads to higher biomass, hence higher possible harvest. If E[Z] = 0, the result is more difficult to explain: for a sufficiently large value of e, it is optimal to use a higher level of critical value for the biomass.

3 The model with updated biomass and growth rate

We now consider a resource submitted two types of random updates (biomass updates Y, growth rate updates R) at random times:

-the times between two biomass updates follows exponential distribution with mean 1/λ x -the times between two growth rate updates follows exponential distribution with mean 1/λ r .

Hence the random update time τ between two updates follows an exponential distribution with mean 1/(λ x + λ r ). For each random time τ , the biomass or the growth rate is updated:

x(τ + ) = Y d ∼ L x (.|x(τ )) with probability λ x λ x + λ r and r(τ + ) = R d ∼ L r (.|r(τ )) with probability λ r λ x + λ r ,
where L x and L r are conditional distributions.

We assume that between two random updates the growth rate does not change, i.e. the grow is piece-wise.

The dynamics of the growth rate is:

ṙ(t) = 0.
As in the previous case, the dynamics of the biomass can be described by a Piecewise Deterministic Markov Process (PDMP). The random jump process are described by the jump kernels Q x and Q r . To each function θ of biomass x and growth rate r, the functions Q x [θ] and Q r [θ] are defined by: Q

x [θ](x, r) = Y θ(Y, r)dL x (Y|x), Q r [θ](x, r) = R θ(x, R)dL r (R|r).

The biomass growth process

In order to describe the biomass growth process we now define the function X(r, t; x, τ ) at time t. If the biomass was x at time τ , the evolution of X(r, t; x, τ ) is given by (S x,r,τ ):

dX(r, t; x, τ ) dt = G(X(r, t; x, τ ), r) -h(X(r, t; x, τ ), e(t)),
with initial condition: X(r, τ ; x, τ ) = x.

The process of the biomass {X t (x 0 , r 0 ) : t ≥ 0} starting at time τ 0 = 0 with x + 0 = x 0 , may be expressed for i ≥ 1:

X t = X(r i-1 , t; x + i-1 , τ i-1 ), τ i-1 < t ≤ τ i ,
where at time τ i :

r i = r i-1 , x + i d ∼ L x (.|X τ i ) with probability λ x λ x + λ r and x + i = X τ i , r i d ∼ L r (.|r i-1
) with probability λ r λ x + λ r .

Given a biomass x and an effort e, we assume a regulator maximizing expected discounted gain on an infinite horizon:

J(x 0 , r 0 , e(.)) = E +∞ 0 l(X t (x 0 , r 0 ), e(t))e -δt dt , with X t (x 0 , r 0 ) solution obtained with successive systems (S x 0 ,r 0 ,0 ), (S x + τ 1 ,rτ 1 ,τ 1 ), .... Thus we consider the function value V defined by: V (x, r) = max e(.)∈[0,e] J(x, r, e(.)).

Using the same formalism than for the updated biomass model and assuming

V ∈ C 1 ([0, K] ⋆ [r, r]), denoting λ.Q[V ] = λ x Q x [V ] + λ r Q r [V ],
we can formally deduce the corresponding Bellman Hamilton Jacobi (BHJ) Equation:

max e∈[0,e] [V ′ x (x, r)g(x, r, e) -(δ + λ x + λ r )V (x, r) + l(x, e) + λ.Q[V ](x, r)] = 0. (10)
The BHJ Equation becomes:

max e∈[0,e] [l 0 (x) -h 0 (x)V ′ x (x, r)]e + V ′ x (x, r)G(x, r) -(δ + λ x + λ r )V (x, r) + λ.Q[V ](x, r) = 0.
In this Equation, the effort e depends on the biomass x and the growth rate r. As 0 ≤ e ≤ e, the optimal effort is a feedback control e * (t) = E(X t , R t ) where the function E is defined by:

E(x, r) =                0, if l 0 (x) -h 0 (x)V ′ x (x, r) < 0, G(x,r) h 0 (x) , if l 0 (x) -h 0 (x)V ′ x (x, r) = 0, e, if l 0 (x) -h 0 (x)V ′ x (x, r) > 0.
Hence:

[l 0 (x) -h 0 (x)V ′ x (x, r)] + e + V ′ x (x, r)G(x, r) -(δ + λ x + λ r )V (x, r) + λ.Q[V ](x, r) = 0. (11)
The critical value x * (r) is the solution of the equation:

l 0 (x) -h 0 (x)V ′ x (x, r) = 0. ( 12
)
But this equation is not sufficient to characterize the critical value. By using a dynamic programming equation, we obtain a complementary condition based on Euler-Lagrange condition.

The value function V (x 0 , r 0 ) is the solution to the optimization problem:

V (x 0 , r 0 ) = J(x 0 , r 0 , e * (.)) = max e(.)∈[0,e] E +∞ 0 l(X t (x 0 , r 0 ), e(t))e -δt dt , with X t (x 0 , r 0 ) solution of the system (S x 0 ,r 0 ,0 ).

Using the same reasoning than in the previous section, we obtain the dynamic programming equation:

V (x 0 , r 0 ) = max e(.)∈[0,e] +∞ 0
[l(X(r 0 , t; x 0 , 0), e(t)) + λ.Q[V ](X(r 0 , t; x 0 , 0), r 0 )]e -(δ+λx+λr)t dt.

Proposition 3.1 Assuming that V, Q x [V ] and Q r [V ] ∈ BC 1 ([0, K] ⋆ [r, r]
), a critical value x * (r) is solu-tion of the system of equations:

l 0 (x) -h 0 (x)V ′ x (x, r) = 0 and (δ + λ x + λ r -G ′ x (x, r)) l 0 h 0 (x) = l 0 h 0 ′ (x)G(x, r) + [λ.Q[V ]] ′ (x), (13)
where the value function V is solution of the Equation (11).

For fixed r, using the same reasoning than in the previous section:

(δ -G ′ x (x, r) + λ x (1 - [Q x [V ]] ′ x (x, r) V ′
x (x, r)

) + λ r (1 - [Q r [V ]] ′ x (x, r) V ′ x (x, r )
)) l 0 h 0 (x) = l 0 h 0 ′ (x)G(x, r)
and replacing λ by λ x (respectively λ r ) and λQ[V ] by

λ x Q x [V ] (respectively λ r Q r [V ]
) we can obtain the equivalent of the Propositions 2.2 and 2.3.

Proposition 3.2

The function value V is twice differentiable with respect to biomass x and:

V ′′ x 2 (x * (r), r) = l 0 h 0 ′ (x * (r)). ( 14
)
Assuming l 0 (x) = pqx -c and h 0 (x) = qx with price p, catchability q and cost c, the critical value x * (r) is an increasing (respectively decreasing) function with respect to biomass jump rate

λ x if [Q x [V ]] ′ x (x * (r), r)- V ′
x (x * (r), r) > 0 (respectively < 0) and is an increasing (respectively decreasing) function with respect to

growth jump rate λ r if [Q r [V ]] ′ x (x * (r), r) -V ′
x (x * (r), r) > 0 (respectively < 0). Moreover the function value V and V ′ x are continuously differentiable with respect to jump rate λ x , and jump growth rate λ r .

In the following section, we will illustrate for a concrete case, with a specific jump kernels for biomass (respectively growth rate), by a study of the sign of

[Q x [V ]] ′ x (x * , r)-V ′ x (x * , r) (respectively [Q r [V ]] ′ x (x * , r)- V ′ x (x * , r)).
For each growth rate r and for sufficiently small values of jump rate λ x , λ r , Equations ( 12) and ( 13) have a unique solution x * (r), the function E is given by:

E(x, r) =                0, if x < x * (r), G(x,r) h 0 (x) , if x = x * (r), e, if x > x * (r), Remark 3.1
The expression of E(x, r) has the same formulation than in the case without updates but due to the difference between the Propositions 2.1 and 3.1, the corresponding x * (r) differs.

Proposition 3.3

The value function V is twice differentiable but not third differentiable and at λ = 0:

V ′′ xr (x * (r), r) = 0, (15)
V ′′′

x 2 r -(x * (r), r) = - Σ 1 (x * (r), r) E(x * (r), r) < 0 < V ′′′ x 2 r + (x * (r), r) = Σ 1 (x * (r), r) e -E(x * (r), r) ( 16
)
and where

Σ 1 (x, r) = δ r l 0 h 2 0 (x).
The third derivatives of the function value V with respect to biomass x and growth rate r are linked by:

V ′′′ x 2 r ± (x * (r), r)x * ′ (r) + V ′′′ xr 2 ± (x * (r), r) = 0. ( 17
)
Proofs are given in Appendix C. In Figure 2, we give an example of optimal evolution of the biomass x with alternatively biomass and growth rate updating. The dash lines represent the successive levels of the critical value x * (r),

x * (r 0 )

x * (r 2 )

x * (r 4 )

x

Y i = x(τ x i )(1+Z i ǫ) (respectively, R i = r(τ r i )(1+S i ξ))
where Z i (respectively, S i ) follows the distribution H x (respectively, H r ).

For a sufficiently small value of biomass jump rate λ x and x * (r) not close to K:

(i) if E[Z] = 0 then the critical value x * (r) is increasing (respectively, decreasing) with respect to biomass jump rate λ x if E[Z] > 0 (respectively, < 0). (ii) if E[Z] = 0 then the critical value x * (r) is increasing (respectively, decreasing) with respect to biomass jump rate λ x if E(x * (r), r)) is larger (respectively, smaller) than e 2 .
For a sufficiently small value of growth jump rate λ r and r not close to r and r:

(iii) the critical value x * (r) is increasing (respectively, decreasing) with respect to growth jump rate λ r if

E(x * (r), r) is larger (respectively, smaller) than e 2 .
We deduce the following properties:

Conclusion

In this article, we consider the evolution of a fishery following a continuous process and submitted to random updates at random times, we present the appropriate PDMP framework. We express the control problem with biomass updates, we highlight the importance of using a dynamic programming approach in order to completely characterize the critical value of the control. We give conditions which permit to deduce the behavior of the optimal control effort with respect to jump rate. An application to a specific jump kernel shows the possible variety of behavior of the optimal effort with respect to the random structure. For a centrally disturbed biomass and sufficiently high effort the optimal harvest increases with biomass jump rate Finally we study the more complex case for which biomass and growth rate in the dynamics are updated.

For a centrally disturbed biomass and sufficiently high effort the optimal harvest increases with the biomass jump rate and decreases with the growth jump rate.

where V (x) is independent of t, hence, formally differentiating with respect to t:

max e∈[0,e] [V ′ (x)g(x, e) -(δ + λ)V (x) + l(x, e) + λQ[V ](x)] = 0.

Appendix B

Proof of Proposition 2.2. In Ω -{x|A(x) = 0}, the value function is smooth and is solution of:

ηA(x)h 0 (x)e + V ′ (x)G(x) -(δ + λ)V (x) + λQ[V ](x) = 0, with η = 0 if A(x) < 0 and η = 1 if A(x) > 0.
By differentiation:

η(A ′ (x)h 0 (x) + A(x)h ′ 0 (x))e + V ′′ (x)G(x) -(δ + λ -G ′ (x))V ′ (x) + λ[Q[V ]] ′ (x) = 0. (18)
Let A ′ -(x * ) and A ′ + (x * ) (respectively, V ′′ -(x * ) and V ′′ + (x * )) the left and right derivatives of A (resp V ′ ) at the critical value x * (if A(x)(x -x * ) > 0 in the vicinity of x * and the reverse if not). From Equations ( 6) and (4), using A(x * ) = 0, at the critical value x * :

V ′′ -(x * )G(x * ) - l 0 h 0 ′ (x * )G(x * ) =0 and A ′ + (x * )h 0 (x * )e + V ′′ + (x * )G(x * ) - l 0 h 0 ′ (x * )G(x * ) =0,
hence respectively using the Equation ( 8):

-A ′ -(x * )G(x * ) = 0 and A ′ + (x * )(h 0 (x * )e -G(x * )) = 0, then A ′ -(x * ) = A ′ + (x * ) = 0, so A is differentiable and V is twice countinuously differentiable.
Proof of Proposition 2.3. In order to determine the behavior of the critical value x * in the vicinity of λ = 0, we differentiate the Equation ( 7) with respect to jump rate λ to obtain:

-G ′′ (x * ) l 0 h 0 (x * ) + (δ -2G ′ (x * )) l 0 h 0 ′ (x * ) - l 0 h 0 ′′ (x * )G(x * ) x * ′ λ (0) = l 0 h 0 (x * )( [Q[V ]] ′ x V ′ x (0, x * ) -1).
Using expression of l 0 (x) and h 0 (x) the second equation becomes:

-G ′′ (x * ) l 0 h 0 (x * ) + (δ + 2( G(x * ) x * -G ′ (x * ))) l 0 h 0 ′ (x * ) x * ′ λ (0) = l 0 h 0 (x * )( [Q[V ]] ′ x V ′ x (0, x * ) -1).
From G ′′ (x) < 0 and G(0) = 0, we deduce that G(x) -xG ′ (x) ≥ 0 and the left term in the brackets is

positive hence x * λ ′ (0) and [Q[V ]] ′ (x) -V ′ (x)
for λ = 0 has the same sign.

We derive Equations (5) for x < x * (λ) and ( 4) with respect to jump rate λ at λ = 0, by using Equation (18):

V ′′ xλ (0, x * )G(x * ) + Q[V ](0, x * ) -V (0, x * ) = δV ′ λ (0, x * ) and l 0 h 0 ′ (x * )x * ′ λ (0) = V ′′ xλ (0, x * ).
Hence: V and V ′ x are continuously differentiable with respect to jump rate λ.

Proof of Proposition 2.4. By differentiation of Equation (18) multiplied by x, for λ = 0:

η[x(A ′ (x)h 0 (x) + A(x)h ′ 0 (x))] ′ e + [x 2 V ′′ (x)] ′ G(x) x -Σ(x) = 0. Hence due to A(x * ) = A ′ (x * ) = 0: [x 2 V ′′ (x)] ′ -(x * )G(x * ) =x * Σ(x * ) and [x 2 V ′′ (x)] ′ + (x * )(G(x * ) -h 0 (x * )e) + (x * x l 0 h 0 ′ (x)h 0 (x) ′ (x * ) -x * 3 V ′′ (x * ) h 0 x ′ (x * ))e =x * Σ(x * ).
From expression of l 0 and h

0 , x l 0 h 0 ′ (x)h 0 (x) ′ (x) = 0 and h 0 x ′ (x) = 0 so: [x 2 V ′′ (x)] ′ + (x * )(G(x * ) -h 0 (x * )e) =x * Σ(x * ).
From expression of E(x * ):

h 0 (x * )[x 2 V ′′ (x)] ′ -(x * )E(x * ) =x * Σ(x * ), h 0 (x * )[x 2 V ′′ (x)] ′ + (x * )(E(x * ) -e) =x * Σ(x * ). From concavity of G, G(x) -xG ′ (x) ≥ 0 and G ′′ (x) < 0, hence Σ(x) > 0, so [x 2 V ′′ (x)] ′ -(x * ) > 0 and [x 2 V ′′ (x)] ′ + (x * ) < 0.
Proof of Proposition 2.5. In order to determine the behaviour of an optimal critical value x * , we differentiate the Equation ( 6) with respect to growth rate r at λ = 0 to obtain:

-G ′′ x 2 (x * , r) l 0 h 0 (x * ) + (δ -2G ′ x (x * , r)) l 0 h 0 ′ (x * ) - l 0 h 0 ′′ (x * )G(x * , r) x * ′ (r) = δ r l 0 h 0 (x), -G ′′ x 2 (x * , r) l 0 h 0 (x * ) + (δ + 2( G(x * , r) x * -G ′ x (x * , r))) l 0 h 0 ′ (x * ) x * ′ (r) = δ r l 0 h 0 (x).
From G ′′ x 2 (x, r) < 0 and G(0) = 0, we deduce G(x, r) -xG ′ x (x, r) ≥ 0 and the left term in the brackets is positive, hence from positive marginal gain.

Proof of Proposition 2.6. For biomass x not close to K: Q

[V ](x) = X V (x(1 + zǫ))dH(z), so: [Q[V ]] ′ (x) -V ′ (x) = X [(1 + zǫ)V ′ (x(1 + zǫ)) -V ′ (x)]dH(z), where V ′ (x(1 + zǫ))(1 + zǫ) -V ′ (x) = zǫ(xV ′′ (x) + V ′ (x)) + zǫ 0 (zǫ -t)[V ′′ ((1 + t)x)(1 + t)x + V ′ ((1 + t)x)] ′ t dt = zǫ(xV ′′ (x) + V ′ (x)) + zǫ 0 (zǫ -t)[V ′′′ ((1 + t)x)(1 + t)x 2 + 2V ′′ ((1 + t)x)x]dt.
For a sufficiently small ǫ, for all x = x * the integral can be approximated by

z 2 ǫ 2 2 (x 2 V ′′′ (x) + 2xV ′′ (x)) (i.e. z 2 ǫ 2 2 [x 2 V ′′ (x)] ′
) in the vicinity of critical biomass x * . Hence at critical biomass x * :

(i) lim ǫ→0 [Q[V ]] ′ (x * ) -V ′ (x * ) ǫ = E[Z](x * V ′′ (x * ) + V ′ (x * )) = E[Z] x l 0 h 0 ′ (x * ), (ii) lim ǫ→0 [Q[V ]] ′ (x * ) -V ′ (x * ) ǫ 2 = E[Z 2 ] 4 ([x 2 V ′′ (x)] ′ -(x * ) + [x 2 V ′′ (x)] ′ + (x * )). From Proposition 2.4, [x 2 V ′′ (x)] ′ -(x * ) + [x 2 V ′′ (x)] ′ + (x * ) = x * Σ(x * ) h 0 (x * ) e -2E(x * ) e -E(x * ) , hence using Propo- sition 2.3, x ′ (0) > 0 if E(x * ) < e 2 and x ′ (0) < 0 if E(x * ) > e 2 .

Appendix C

Proof of Proposition 3.3. From total differentiation of Equation ( 12) with respect to growth rate r:

(l ′ 0 (x * (r)) -h ′ 0 (x * (r))V ′′ xx (x * (r), r))x * ′ (r) = h 0 (x * (r))V ′′ xr (x * (r), r), and using Equation ( 14) we deduce Equation (15). Let A(x, r) = l 0 h 0 (x) -V ′ x (x, r).

In Ω -{x|A(x, r) = 0} the value function is smooth and is solution of:

ηA(x, r)h 0 (x)e + V ′ x (x, r)G(x, r) -

(δ + λ x + λ r )V (x, r) + λ.Q[V ](x, r) = 0, (19)
with η = 0 if A(x, r) < 0 (i.e. x < x * (r)) and η = 1 if A(x, r) > 0 (i.e. x > x * (r)).

Using Equation ( 14), from differentiation of Equation ( 19) with respect to biomass x and growth rate r at λ x = λ r = 0, at optimal critical value x * (r):

(G(x * (r), r) -ηh 0 (x * (r))e)V ′′′ x 2 r (x * (r), r) + V ′′ x 2 (x * (r), r)G ′ r (x * (r), r) + V ′ x (x * (r), r)G ′′ rx (x * (r), r) = 0, (G(x * (r), r) -ηh 0 (x * (r))e)V ′′′ From total differentiation of Equation ( 15) with respect to growth rate r we deduce (17).

Proof of Proposition 3.4 (i) and (ii) Using the same reasoning than in the previous section and replacing λ by λ x and λQ[V ] by λ x Q x [V ] we obtain the same result than in the Proposition 2.6.

(iii) For growth rate r not close to r and r:

Q r [V ](x, r) -V (x, r) = For a sufficiently small ξ, the integral can be approximated by s 2 ξ 2 2 r 2 (V ′′′ xr 2 ) {-sign(s)} (x, r) in the vicinity of the critical biomass x * (r) for all growth rate r. Hence, due to Equation ( 15) : + (x * , r))+O(ξ 3 ). From Equations ( 16) and ( 17), V ′′′ is proportional with the same sign to e -2E(x * (r), r).

[Q r [V ]] ′ x (x, r)-V ′ x (x, r) = E[R 2 ]ξ 2 4 r 2 (V ′′′
5

Figure 2 :

 2 Figure 2: Optimal evolution of biomass x with biomass and growth rate updated

3. 2

 2 Application to specific jump kernels Proposition 3.4 Assuming the updated biomass Y i (respectively, the updated growth rate R i ) given by:

Corollary 3 . 1

 31 (i) The behavior of the critical value with respect to biomass jump rate λ x is of the same type that in the previous case with update of the biomass. (ii) If E[Z] = 0, the behavior of the critical value x * (r) is reversed with respect the biomass jump rate λ x and the growth jump rate λ r . (iii) The behavior of the critical value x * (r) with respect to growth jump rate λ r is independent of the expectation E[ξ].

x 2 r

 2 (x * (r), r) (r))G(x * (r), r) + l 0 h 0 (x * (r))G ′ x (x * (r), r)) = 0, then using Equation (13) (G(x * (r), r) -ηh 0 (x * (r))e)V ′′′ x 2 r (x * (r), r) (r)) = 0, h 0 (x * (r))(E(x * (r), r) -ηe)V ′′′ x 2 r (x * (r),r) positive marginal gain we deduce the expression of V ′′′ x 2 r -(x * (r)) and V ′′′ x 2 r + (x * (r)).

  x, r(1 + sξ)) -V (x, r))dH r (s), hence: [Q r [V ]] ′ x (x, r) -V ′ x (x, r) = r r [V ′ x (x, r(1 + sξ)) -V ′ x (x, r)]dH r (s)and:V ′ x (x, r(1 + sξ)) -V ′ x (x, r) = sξrV ′′ xr (x, r) + sξ 0 (sξ -t)r 2 V ′′′ xr 2 (x, r(1 + t))dt.

xr 2 -

 2 (x * , r)+V ′′′ xr 2

xr 2 - 2 +

 22 (x * , r)+V ′′′ xr (x * , r)

Figure 1 :Figure 2 :

 12 Figure 1: Optimal evolution of biomass x with biomass updated

Appendix A

Let the value function defined by: V (x) = max e(.)∈[0,e] J(x, e(.)). Consider a time t > 0, by the strong Markov property, the criteria satisfies:

), e * (u))e -δu du + V (X(t, x))e -δt )I t<τ , E τ τ 0 l(X(u, x), e * (u))e -δu duI τ <t = λ t 0 τ 0 l(X(u, x), e * (u))e -δu du e -λτ dτ then by inverting integration with respect to t and τ = t 0 l(X(u, x), e * (u))(e -(δ+λ)u -e -λt-δu )du, and hence, rearranging the terms:

[l(X(τ, x), e(τ )) + λQ[V ](X(τ, x))]e -(δ+λ)τ dτ + e -(δ+λ)t V (X(t, x)) ,


Gaétan Bahl

Lionel Daniel

Matthieu Moretti
email: matthieu.moretti@irt-saintexupery.com

Florent Lafarge
email: florent.lafarge@inria.fr

Low-power neural networks for semantic segmentation of satellite images

Semantic segmentation methods have made impressive progress with deep learning. However, while achieving higher and higher accuracy, state-of-the-art neural networks overlook the complexity of architectures, which typically feature dozens of millions of trainable parameters. Consequently, these networks requires high computational ressources and are mostly not suited to perform on edge devices with tight resource constraints, such as phones, drones, or satellites. In this work, we propose two highlycompact neural network architectures for semantic segmentation of images, which are up to 100 000 times less complex than state-of-the-art architectures while approaching their accuracy. To decrease the complexity of existing networks, our main ideas consist in exploiting lightweight encoders and decoders with depth-wise separable convolutions and decreasing memory usage with the removal of skip connections between encoder and decoder. Our architectures are designed to be implemented on a basic FPGA such as the one featured on the Intel Altera Cyclone V family of SoCs. We demonstrate the potential of our solutions in the case of binary segmentation of remote sensing images, in particular for extracting clouds and trees from RGB satellite images.

Introduction

The semantic segmentation of images, which consists in assigning a semantic label to each pixel of an image, is a long standing problem in Computer Vision. The popularization of Convolutional Neural Networks (CNN) has led to a lot of progress in this field. Fully Convolutional Networks (FCN) [23] and related architectures, such as the popular U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], outperformed traditional methods in terms of accuracy while being easy to use.

In the quest towards accuracy, state-of-the-art neural networks often disregard the complexity of their architectures, which typically feature millions of trainable parameters, requiring days of training and gigabytes hard-drive space.

These architectures also require a lot of computing power, such as GPUs consuming hundreds of Watts. Consequently, these architectures cannot be easily used on edge devices such as phones, drones or satellites.

Our goal is to design neural network architectures operating on low-power edge devices with the following objectives:

• High accuracy. The architecture should deliver accurate semantic segmentation results. • Low complexity. The architecture should have a low number of parameters to learn. • Adaptability. The architecture should be easily implementable on low-power devices. Some efficient architectures, such as ESPNet [START_REF] Mehta | Espnet: Efficient spatial pyramid of dilated convolutions for semantic segmentation[END_REF], and convolution modules [START_REF] Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF] have been proposed for phones and autonomous vehicles. These architectures are however too complex to be exploited on devices with lower power budget and computational resources. For instance, this is the case of FPGA cards embedded into system-on-chips (SoCs) on satellites, which have limited floating-point capabilities and are restricted by their amount of hardware logic.

In this paper, we present two highly-compact neural network architectures for semantic segmentation, which are up to 100 000 times less complex than state-of-the-art architectures while approaching their accuracy. We propose two fully-convolutional neural networks, C-FCN and C-UNet, which are compact versions of the popular FCN [23] and U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] architectures. To decrease the complexity of existing networks, our main ideas consist in exploiting lightweight encoders and decoders with depth-wise separable convolutions and decreasing memory usage with the removal of skip connections between encoder and decoder. Our architectures are designed to be implemented on a basic FPGA such as the one featured on the Intel Altera Cyclone V family of SoCs. We demonstrate the potential of our solutions in the case of binary segmentation of remote sensing images, in particular for extracting clouds and trees from RGB satellite images. Our methods reach 95% accuracy in cloud segmentation on 38-Cloud [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF] and CloudPeru2 [START_REF] Ghassemi | Convolutional Neural Networks for On-Board Cloud Screening[END_REF] datasets, and 83% in forest segmentation on EOLearn Slovenia 2017 dataset.

Related works

We first review prior work in two related aspects of our goal: neural networks for semantic segmentation and neural network inference on edge devices.

Neural networks for semantic segmentation

Many deep learning methods have been proposed to improve performance of segmentation networks, in terms of accuracy on popular benchmarks such as [6, [START_REF] Lin | Microsoft COCO: Common objects in context[END_REF][START_REF] Mottaghi | The role of context for object detection and semantic segmentation in the wild[END_REF] and speed, with the goal of real-time semantic segmentation [START_REF] Siam | RTSeg: Real-Time Semantic Segmentation Comparative Study[END_REF]. These methods are compared in recent surveys [START_REF] Guo | A review of semantic segmentation using deep neural networks[END_REF]7,[START_REF] Siam | RTSeg: Real-Time Semantic Segmentation Comparative Study[END_REF] without studying the networks from a complexity point of view.

Encoder-decoders. Most of the popular segmentation neural networks, such as U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], SegNet [2], or DeepLabv3+ [START_REF] Chen | Encoder-decoder with atrous separable convolution for semantic image segmentation[END_REF] adopt an encoder-decoder architecture. The encoder part of the network, also called backbone, is typically a Deep Convolutional Neural Network composed of multiple stages of convolution operations, separated by pooling operations. This allows the encoder to capture highlevel features at different scales. It is common to use an existing CNN architecture such as ResNet-101 [START_REF] He | Deep residual learning for image recognition[END_REF], VGG16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF], or MobileNet [START_REF] Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF] as a backbone, since pre-trained weights on databases such as ImageNet [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF] are available for transfer learning. However, while achieving great accuracy, these encoders are complex and not suited for inference on edge devices. Decoders are in charge of upsampling the result to get an output that has the same size as the original image. Deconvolutions, also called transposed convolutions or up-convolutions, introduced by H. Noh et al. [START_REF] Noh | Learning deconvolution network for semantic segmentation[END_REF] in the context of semantic segmentation, are used to learn how images should be upscaled. Skip connections. Skip connections are often used between convolution blocks, i.e. short skips, in architectures such as residual networks [START_REF] He | Deep residual learning for image recognition[END_REF], and/or between the encoder and the decoder, i.e. long skips, in architectures such as U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. These connections are performed by concatenating feature maps from different parts of the network in order to retain some information. In the case of segmentation networks, skip connections are used to keep high-frequency information (e.g. corners of buildings, borders of objects) to obtain a more precise upscaling by the decoder. While skip connections can yield good results in practice, they are often memory consuming, especially the long skips, as they require feature maps to be kept in memory for a long time. Consequently, they cannot easily be used on edge devices with limited memory. In particular, we show in Section 4 that they are not useful for our use cases.

Neural networks on edge devices

Efficient convolution modules. Several low-complexity convolution modules have been developed, with the goal of reducing power consumption and increase inference speed on low-end or embedded devices, such as phones. This is the case of Inception [START_REF] Szegedy | Going deeper with convolutions[END_REF], Xception [START_REF] Chollet | Xception: Deep learning with depthwise separable convolutions[END_REF], ShuffleNet [START_REF] Zhang | ShuffleNet: An Extremely Efficient Convolutional Neural Network for Mobile Devices[END_REF], or ESP [START_REF] Mehta | Espnet: Efficient spatial pyramid of dilated convolutions for semantic segmentation[END_REF]. These modules are based on the principle of convolution factorization and turn classical convolution operations into multiple simpler convolutions. S. Mehta et al. provide a detailed analysis of these modules in [START_REF] Mehta | Espnet: Efficient spatial pyramid of dilated convolutions for semantic segmentation[END_REF]. However, we cannot use them as they still feature too many trainable parameters and often require to store the results of multiple convolution operations performed in parallel or short skip connections. In contrast, we use Depth-wise Separable Convolutions, introduced in [START_REF] Sifre | Rigid-motion scattering for image classification[END_REF] and used in MobileNets [START_REF] Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF].

Neural Network simplification Several techniques have been developed to simplify neural networks for systems with tight resource constraints. Training and inference using fixed-point or integer arithmetic and quantization is an active research topic [20, [START_REF] Jacob | Quantization and training of neural networks for efficient integer-arithmetic-only inference[END_REF][START_REF] Wu | Training and inference with integers in deep neural networks[END_REF], since floating-point operations are costly on many embedded systems such as FPGAs, and trade-offs have to be made between accuracy and speed. Neural network pruning, on the other hand, is the process of removing some weights [START_REF] Lecun | Optimal Brain Damage[END_REF] or entire convolution filters [START_REF] Li | Pruning filters for efficient convnets[END_REF]27] and their associated feature maps in a neural network, in order to extract a functional "sub-network" that has a lower computational complexity and similar accuracy. A lot of work has been done on enabling and accelerating NN inference on FPGA [START_REF] Abdelouahab | Accelerating CNN inference on FPGAs: A survey[END_REF], with tools being released to automatically generate HDL code for any neural network architecture [START_REF] Hamdan | VHDL auto-generation tool for optimized hardware acceleration of convolutional neural networks on FPGA[END_REF], with automated use of quantization and other simplification techniques. These approaches are orthogonal to our work as they could be applied to any neural network.

Proposed architectures

We now detail our compact neural networks, C-UNet and C-FCN. These architectures are Fully Convolutional. Thus the size of the networks does not depend on the size of the input images. In particular, the training steps and the testing steps can be performed on different image sizes.

C-UNet

As illustrated in Figure 1, C-UNet is a compact version of the popular U-Net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF].

To create a shallower network, three convolution stages have been removed with respect to the original architecture. Removing stages strongly reduces the computational cost as the number of feature maps is typically multiplied by a factor 2 at each stage. For instance, U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] has 1024 feature maps at the last encoder stage. A standard 3x3 convolution at this stage uses 3×3×1024×1024+1024 = 9438208 parameters. The number of stages to remove has been chosen after experimental validation, 3 being a good compromise between field of view and number of parameters. Indeed, we found experimentally that we could train a classifier to detect objects such as clouds and trees in 28x28 pixel images with good accuracy. The encoder in C-UNet has a receptive field of 50x50 pixels and is consequently more than capable of reaching the same accuracy.

Some of the standard convolution layers have also been replaced by depth-wise separable convolutions. Introduced in [START_REF] Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF], they allow a classical convolution to be split into a depth-wise convolution and a point-wise convolution, to use the fact that inter-channel features are often not related to spatial features. Since the convolution kernel is the same for all input feature maps in a depth-wise convolution, this allows us to reduce the number of learned parameters significantly. For example, a standard 3x3 convolution with 16 input and output feature maps has 2320 trainable weights, while a depth-wise separable convolution of the same size only has 416 trainable weights.

The number of filters per convolution has been chosen so that C-UNet has 51 113 parameters, ie around 500 times less than the architecture implemented in [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF]. C-UNet then has 8 filters per convolution at the first stage and this number is doubled at every stage in the encoder, and divided by two at every stage in the decoder, as in the original U-Net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF].

We also propose a variant of C-UNet, called C-UNet++, with an even more compact architecture. C-UNet++ contains 9 129 parameters, ie 3000 times less parameters than the architecture implemented in [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF]. As illustrated in Figure 1, the number of convolution layers has been reduced to one per stage, and only the first stage has standard convolutions. This gives C-UNet++ a receptive field of 22x22 pixels.

Note that skip connections between corresponding stages in encoder and decoder parts could be considered. However, we show in Section 4 that they are not necessarily useful, as adding them may not significantly increase performance and may require a lot of memory depending on the input image size.

C-FCN

The second proposed architecture, called C-FCN, relies upon a small CNN and a bilinear upsampler, as illustrated in Figure 2. The encoder has 3 stages, with a single convolution per stage. A 1x1 convolution then generates a class heatmap which is then upsampled by the bilinear upsampler. C-FCN uses depthwise separable convolutions to maximize the number of convolution filters we can use within our parameter budget, which is of 1 438 parameters for this architecture.

We also propose a variant of C-FCN with the same depth, called C-FCN++, which is designed to be the smallest viable architecture with only 273 parameters. The first layer is an atrous convolution layer with a dilation rate of 2, which broadens the receptive field of the network a little bit [START_REF] Yu | Multi-scale context aggregation by dilated convolutions[END_REF].

The number of filters of C-FCN and its variant C-FCN++ are showed in Table 1.

3 F 1 F 2 F 3 1 Figure 2: C-FCN

Experiments on cloud segmentation

In this section, we evaluate the performance of our architectures for cloud extraction in RGB remote sensing images. Performing real-time cloud segmentation directly onboard and removing useless cloud-covered images can allow significant bandwidth, storage and computation time savings on the ground. Accurate detection of clouds is not an easy task, especially when a limited number of spectral bands is available, as clouds share the same radiometric properties as snow, for example. Traditional cloud segmentation methods are either threshold-based or handcrafted, with the most popular ones being Fmask [START_REF] Zhu | Improvement and expansion of the fmask algorithm: Cloud, cloud shadow, and snow detection for landsats 4-7, 8, and sentinel 2 images[END_REF] and Haze Optimized Transform [START_REF] Zhang | An image transform to characterize and compensate for spatial variations in thin cloud contamination of landsat images[END_REF].

M. Hughes et al. [START_REF] Hughes | Automated Detection of Cloud and Cloud Shadow in Single-Date Landsat Imagery Using Neural Networks and Spatial Post-Processing[END_REF] were, to our knowledge, the first to use a neural network-based approach for cloud segmentation in Landsat-8 images. However, convolutional neural networks were not used. A lot of recent work has been done on cloud segmentation using CNN [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF][START_REF] Mohajerani | Cloud-net: An endto-end cloud detection algorithm for landsat 8 imagery[END_REF][START_REF] Ghassemi | Convolutional Neural Networks for On-Board Cloud Screening[END_REF][START_REF] Drönner | Fast cloud segmentation using convolutional neural networks[END_REF]22]. However, all these methods use computationnaly-heavy architectures on hyperspectral images and are designed to be used with powerful hardware on the ground, and thus, are not compatible with our use-case. Cloud segmentation using neural networks has already been integrated in an ARTSU CubeSat mission by Z. Zhang et al. [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF]. S. Ghassemi et al. [8] have proposed a small strided U-Net architecture for onboard cloud segmentation. Their architectures are still too big for our use cases and are designed to work on 4-band images (RGB + NIR). Nevertheless, we include a 3-band implementation of MobUNet [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF], MobDe-convNet [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF] and "Plain+" strided U-Net [8] architectures in our comparison.

Datasets. We use publicly available cloud segmentation datasets to train and test the models.

We use the Cloud-38 dataset, introduced by S. Mohajerani et al. [START_REF] Mohajerani | Cloud-net: An endto-end cloud detection algorithm for landsat 8 imagery[END_REF]. The dataset is composed of 4-band Landsat-8 images with a 30m resolution. We only use RGB bands. Since some of the training scenes have black borders in the original dataset, we remove all the patches that have more than 50% of black pixels. The training and validation set is composed of 4 821 patches of size 384x384. The testing set is composed of 9 200 patches.

We also use the CloudPeru2 dataset, introduced by G. Morales et al. [START_REF] Ghassemi | Convolutional Neural Networks for On-Board Cloud Screening[END_REF]. This dataset is composed of 4-band PERUSAT-1 images with a 3m resolution. It features 22 400 images of size 512x512. We use 10% of these images for testing. Test images have been chosen randomly and are the same for all of our tests.

Training procedure. All models were implemented using the Keras framework using Tensorflow 1.13 as a backend in Python 3.6.

Models are trained using a round-based scheme. Training sets are randomly shuffled and split into training (85%), and validation (15%) sets at the beginning of each round. Within a round, all models are then trained on the same training and validation sets for a maximum of 150 epochs. At the end of all rounds, the best network for each architecture is selected using the testing set.

Data is augmented using random horizontal and vertical flips, as well as random rotations of 5 degrees maximum. The same random seed is used for all architectures. Models were trained for 4 rounds with batch-size 8 and the Adam optimizer with a learning rate of 0.001. Binary crossentropy loss is used for all networks. Early stopping is used with a patience of 15 and a maximum number of epochs of 150.

This training scheme takes around 10h per architecture, per round, per dataset on an Nvidia GTX 1080 Ti graphics card, which brings our total training time to around 1160h. Training could have been faster, as the card's VRAM was barely being used when training really small models, but we wanted to use the same training parameters (batch-size) for every model, to get the fairest possible comparison. The largest batch-size we could use to train the biggest network (U-Net) is 8.

Comparative architectures. We compare our networks to different semantic segmentation architectures from the litterature, that we re-implemented:

• U-Net for cloud segmentation, as presented by S. Mohajerani et al. [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF], which we use as a "baseline" for comparison, since U-Net, originally presented by Ronneberger et al. [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] for the segmentation of medical images is a well-known and easy to implement architecture that has proved useful for many segmentation tasks, • ESPNet A, introduced by S. Mehta et al. [START_REF] Mehta | Espnet: Efficient spatial pyramid of dilated convolutions for semantic segmentation[END_REF], which is a state-of-the-art segmentation network, in terms of efficiency and accuracy. We use K = 5, α 2 = 2, α 3 = 5. • MobUNet and MobDeconvNet, introduced by Z.

Zhang et al. [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF], are small versions of U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] and Deconv-Net [START_REF] Noh | Learning deconvolution network for semantic segmentation[END_REF] for onboard cloud segmentation, that use depth-wise separable convolutions [START_REF] Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF], • StridedUNet, introduced by Ghassemi et al.

[8] for onboard cloud segmentation, which uses strided convolutions for downsampling, • LeNetFCN, an FCN variation of the well-known LeNet-5 architecture [START_REF] Lecun | Gradientbased learning applied to document recognition[END_REF], which was originally created for written number classification. We replace the final Fully-Connected (Dense) layers by a 1x1 2D Convolution with a sigmoid activation, in order to output a segmentation map. LeNetFCN architecture looks like the "C-FCN" architecture (see Fig. 2) with F 1 = 3, F 2 = 5, except it uses 5x5 convolutions instead of 3x3 convolutions. We evaluate models using standard segmentation metrics: Total Accuracy, Precision, Recall, Specificity and Jaccard Index (IoU), as defined in [START_REF] Mohajerani | Cloud-net: An endto-end cloud detection algorithm for landsat 8 imagery[END_REF].

Qualitative results. Fig. 3 shows inference results of our networks on three image patches taken from the 38-Cloud test set, compared to the ground truth and U-Net. Our networks provide a visually good result, even on difficult terrains such as snow. The segmentation masks produced appear smooth, especially for C-FCN/C-FCN++ because of the bilinear upsampling. This is not an issue to estimate cloud coverage as a percentage.

Quantitative results. particular, this network matches Fmask's accuracy (within 1%), and achieve much better precision and IoU. This network also exceeds the performance of our implementation of mobDeconvNet [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF], using 4x less parameters. Our C-FCN++ matches the performance of mobDeconvNet while using only 273 parameters. Our C-UNet matches the performance of our implementation of mobUNet [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF], while using less memory thanks to the lack of skip connections. It also exceeds the performance of our implementation of StridedUNet [8], while being 20% smaller in terms of number of parameters. We see similar results on CloudPeru2 as on 38-Cloud. Our networks provide good results on this high-resolution dataset (3m). Fig. 4 shows the relation between the number of parameters of a network and its accuracy. Our networks match or exceed the performance or other methods at each complexity level.

Impact of skip connections. We tried adding skip connections in C-UNet and C-UNet++, but did not find a significant difference in performance metrics in the case of cloud segmentation, as shown in Table 3. Skip connections are useful to keep high frequency information and use it in the decoder part of the network, since the encoder part tends to act like a low-pass filter because of successive convolutions. Clouds do not have much high frequency information, which is maybe why skip connections were not particularly useful in our case. This means that a significant amount of memory can be saved during network inference by removing skip connections in use cases where they are not useful, as shown in Table 4. This can allow bigger networks to be used on systems with tight resource constraints. Performance. Many edge devices use ARM-based SoCs, which are known for their efficiency. This is why we chose the Raspberry Pi, an ARM computer the size of a credit card, for our testing. Our particular board is the Raspberry Pi 4 model B, which features a Quad core Cortex-A72 CPU clocked at 1.5GHz and 2GB of LPDDR4 SDRAM. We feel that this board is representative of the kind of power a typical low-power system could have, as this board consumes 3W when idle and 6W under load on average. We use Tensorflow 1.13.1 with Python 3.7.3 on Raspbian 10.0 to run our evaluation code on the board. The number of parameters, or trainable weights, of an architecture is an important metric, as it is not only linked to the storage space needed for these parameters but also (albeit loosely) correlated to computation time since each weight must be used in computation at some point. On systems such as FPGAs where networks can be "hardcoded" into a chip, the number of parameters is also linked to the amount of hardware logic that will be used by a network.

Table 5 shows the size used by the tested architectures in terms of number of parameters and storage space, as well as the number of FLOPs (Floating Point Operations) used for one forward pass on a 384x384 image. We can see that the number of FLOPs is indeed correlated to the number of parameters. Table 6 shows the execution time results we obtained. We can see that our networks do not run out of memory when trying to process relatively big images of 2048x2048 pixels at once. Even with a Python 32-bit float implementation, our architectures allow fast processing of images. This processing could be even faster by using Tensorflow Lite and quantized weights.

We notice a strange behaviour with stridedUNet, which performs significantly faster than networks that require less FLOPs. We suspect that Tensorflow, or at least the ARM version of it, is lacking some optimizations in the case of depthwise separable and atrous convolutions.

Adaptability on FPGA We successfully implemented our most compact network, C-FCN++, on an Altera Cyclone V 5CSXC6 FPGA with a dataflow/streaming architecture and 16-bit quantization based on VGT [START_REF] Hamdan | VHDL auto-generation tool for optimized hardware acceleration of convolutional neural networks on FPGA[END_REF], a tool that automatically generates HDL code for a given neural network architecture. This FPGA will be included aboard OPS-SAT, a CubeSat satellite that will be launched by the end of 2019 by the European Space Agency. Inference on a 2048x1944 image was done in less than 150ms at 100MHz. Thus, C-FCN++ allows real-time processing of images from OPS-SAT's camera as it captures 2000x2000 images at 7 fps. Comparison with the three other networks on FPGA were not possible as they are not compact enough to be implemented on this FPGA.

Experiments on forest segmentation

To provide a comparison on another use case, we also evaluate our models on a forest segmentation task, which is also important in remote sensing, for applications such as deforestation monitoring. Although vegetation detection is usually done using dedicated spectral bands on earth observation satellites, it would be useful to be able to perform it on devices that do not usually possess these bands, such as drones or nanosats. Neural networks have already been used successfully for land cover and vegetation segmentation [START_REF] Zhu | Deep learning in remote sensing: A comprehensive review and list of resources[END_REF][START_REF] Kamilaris | Deep learning in agriculture: A survey[END_REF]. However, to our knowledge, these tasks have never been done directly onboard.

For this experiment, we use the EOLearn Slovenia 2017 dataset [35], which is composed of multiple Sentinel acquisitions of Slovenia over 2017, with pixel-wise land cover ground truth for 10 classes, as well as cloud masks. We only use the class and the RGB bands. This dataset is originally split into 293 1000x1000 patches. We remove images that have more than 10% of clouds and split the remaining images into 500x500 pixel tiles. We use 24 patches as a test set. This means our training/validation set has 12 629 tiles, and our test set has 3 320 tiles. We use the same training scheme and settings as in our cloud segmentation experiments, which we detailed in 4.

In this test, the performance of our networks, as shown in Table 7, is comparable to the other networks, which proves their adaptability to different tasks. In particular, C-UNet matches the performance of U-Net [START_REF] Mohajerani | A cloud detection algorithm for remote sensing images using fully convolutional neural networks[END_REF], and C-UNet++ matches the performance of mobDeconvNet [START_REF] Zhang | Cubesat cloud detection based on JPEG2000 compression and deep learning[END_REF].

Conclusion and future works

In this paper, we introduced lightweight neural network architectures for onboard semantic segmentation, and compared them to state-of-the-art architectures in the context of 3-band cloud and forest segmentation. We showed that our networks match the performance of Fmask, a popular cloud segmentation method, while only using 3 bands and at a low computational cost that allow them to be implemented on systems with restrictive power and storage constraints. We also match the performance of network architectures that have previously been used for onboard cloud detection, while using less memory and FLOPs. We talked about the usefulness of depth-wise separable convolutions for implementation of neural networks on edge devices and discussed the memory footprint of skip connections. Our work has some limitations which leave room for future work. We only tested our architectures on binary image segmentation, and while this is enough for our test cases, many applications require multi-class image segmentation.

Figure 1 :

 1 Figure 1: U-Net [31] compared to C-UNet and C-UNet++. Gray: input, white: conv3x3 + ReLU, yellow: depthwise separable conv3x3 + ReLU, green: 2x2 max pool, orange: conv1x1 sigmoid, blue: 2x2 deconvolution, red: output, arrows: skip connections. Numbers indicate the number of feature maps at each stage.

Figure 3 :

 3 Figure 3: Qualitative results of C-UNet, C-FCN, and C-FCN++ on 38-Cloud dataset, compared to Ground Truth and U-Net [25]. Our networks in bold font. Our networks give good results and are not fooled by snow on the third patch (m). C-FCN and C-FCN++ results appear very smooth compared to others because the segmentation map output by the network is 4 times smaller than U-Net's and upsampled with a bilinear upsampling. Thus, C-UNet produces a more refined segmentation than C-FCN.

Figure 4 :

 4 Figure4: Test accuracy with respect to number of parameters on the two tested datasets. Our networks in green and bold font, others in blue. Number of parameters is on a log scale. Our networks offer good accuracy at varying levels of complexity. There is a point of diminishing return in this task: our C-UNet performs about the same as a 500 times bigger U-Net.

Table 1 :

 1 Characteristics of C-FCN and C-FCN++. DW = usage of depthwise convolutions. Fx is the number of convolution filters and feature maps at stage x.

		and C-FCN++ architectures. Gray: input, white:
	conv3x3 ReLU, green: 2x2 max pool, orange: conv1x1 sigmoid, blue:
	4x4 bilinear upscaling.
	Arch.	F 1	F 2	F 3	DW	Atrous	Nparam
	C-FCN	10	20	40	Yes	No	1438
	C-FCN++	5	2	2	No	1 st conv.	273

Table 2 :

 2  and CloudPeru2 (bottom) results. Our networks in bold font.

	Model	Acc.	Prec.	Rec.	Spec.	Jacc.
	LeNet FCN	91.10	86.00	84.34	94.03	74.16
	C-FCN++ 93.23	91.67	85.45	96.62	79.30
	mobDeconvNet [41] 93.44	92.04	85.75	96.78	79.83
	C-FCN 93.91	91.23	88.39	96.31	81.47
	C-UNet++ 94.85	94.18	88.48	97.63	83.89
	ESPNet A [24] 95.08	93.94	89.55	97.49	84.66
	StridedUNet [8] 95.39	94.21	90.36	97.58	85.60
	mobUNet [41]	95.58	95.37	89.78	98.11	86.03
	U-Net [25] 95.61	95.69	89.56	98.25	86.08
	C-UNet 95.78	96.53	89.27	98.61	86.50
	Fmask [43] 94.89	77.71	97.22	93.96	75.16
	LeNet FCN	89.23	93.11	83.27	94.52	78.44
	C-FCN++ 91.37	93.50	87.76	94.58	82.71
	mobDeconvNet [41] 93.40	94.98	90.78	95.73	86.63
	C-FCN 94.09	95.15	92.15	95.82	88.01
	C-UNet++ 94.59	96.03	92.31	96.61	88.92
	U-Net [25] 95.82	96.29	94.78	96.75	91.44
	mobUNet [41]	95.97	97.20	94.14	97.59	91.66
	C-UNet 96.42	96.54	95.83	96.94	92.65
	stridedUNet [8] 96.45	96.49	95.95	96.89	92.72
	ESPNet A [24] 96.59	96.45	96.30	96.85	93.01

Table

2

shows accuracy metrics on 38-Cloud and CloudPeru2. We include the Fmask accuracy results from S. Mohajerani et al.

[START_REF] Mohajerani | Cloud-net: An endto-end cloud detection algorithm for landsat 8 imagery[END_REF]

for 38-Cloud, since our test dataset and metrics are the same.

The best network under 1500 parameters is C-FCN. In

Table 3 :

 3 Impact of skip connections on performance of C-UNet(++) on 38-Cloud dataset. We do not see a significant difference in accuracy when using skip connections for this cloud segmentation task.

	C-UNet	Acc.	Prec.	Rec.	Spec.	Jac.
	with skips 95.90	96.22 90.01	98.46	86.93
	w/o skips 95.78	96.53 89.27	98.61	86.50
	C-UNet++	Acc.	Prec.	Rec.	Spec.	Jac.
	with skips 94.51	92.30 89.33	96.76	83.14
	w/o skips 94.85	94.18 88.48	97.63	83.89
	Architecture / Image size	384x384	2000x2000
	U-Net [25]		34.9 MB	236.5 MB
	StridedUNet [8]		4.22 MB	114.4 MB
	MobUNet [41]		1.69 MB	45.8 MB
	C-UNet with skips	1.97 MB	53.4 MB
	C-UNet++ with skips	1.97 MB	53.4 MB
	C-UNet without skips	0 MB	0 MB
	C-UNet++ without skips	0 MB	0 MB

Table 4 :

 4 Memory footprint of skip connections for 32-bit float inference (in Megabytes), computed for two images sizes.

Table 5 :

 5 Network parameters, FLOPs (FLoating-point OPerationS) and storage size. Our architectures in bold font. FLOPs computed for a 384x384 input using Tensorflow's profiler. We see that the number of FLOPs is correlated to the number of paramters. Storage size is the size of the Keras h5 file for each model in MegaBytes.

	Architecture	Nparams	FLOPs	Storage (MB)
	C-FCN++	273	4 654	0.047
	LeNet FCN	614	10 455	0.049
	C-FCN	1 438	24 233	0.066
	mobDeconvNet [41]	7 919	134 256	0.149
	C-UNet++	9 129	154 800	0.172
	C-UNet	51 113	867 712	0.735
	mobUNet [41]	52 657	893 882	0.724
	StridedUNet [8]	79 785	1 355 704	1.0
	ESPNet A [24]	200 193	3 399 310	2.7
	U-Net [25]	31 094 497 528 578 588	357
	Model	384x384 1024x1024	2048x2048
	C-FCN++	0.130	0.773	2.663
	stridedUNet [8]	0.133	0.752	2.666
	C-FCN	0.148	0.903	3.017
	C-UNet++	0.204	1.242	4.295
	C-UNet	0.404	2.867	10.355
	mobUNet [41]	0.688	5.507	OOM
	mobDeconvNet [41]	0.755	N/A	N/A
	ESPNet A [24]	1.878	15.940	OOM
	U-Net [25]	5.035	OOM	OOM

Table 6 :

 6 Execution time on images of different sizes in seconds on Raspberry Pi 4, averaged over 20 iterations. OOM = Out of Memory, N/A = mobDeconvNet is not an FCN and thus cannot be executed on images of different sizes than training images.

Table 7 :

 7 Slovenia 2017 forest segmentation results (10m resolution).

	Model	Acc.	Prec.	Rec.	Spec.	Jacc.
	LeNetFCN 77.10	67.93	63.57	84.22	48.89
	C-FCN++ 77.40	69.52	61.28	85.87	48.31
	C-FCN 78.87	68.95	70.38	83.33	53.44
	mobDeconvNet [41] 80.61	71.52	72.68	84.78	56.36
	C-UNet++ 80.62	72.92	69.61	86.41	55.31
	stridedUNet [8] 80.67	72.52	70.70	85.92	55.77
	C-UNet 83.33	76.79	73.99	88.24	60.47
	U-Net [25] 84.04	73.37	84.28	83.92	64.54
	mobUNet [41]	84.27	75.19	81.13	85.92	63.99

ACKNOWLEDGMENT

This work was funded by IRT Saint-Exupéry. We thank the European Space Agency for providing us the FPGA development kit and for supporting our experiment. We would like to thank S. Mohajerani et al., G. Morales et al. and Sinergise Ltd. for making the 38-Cloud, CloudPeru2 and Slovenia 2017 datasets publicly available.