domain: math sde

text:
Rafael Bravo De La Parra
email: rafael.bravo@uah.es

Jean-Christophe Poggiale
email: jean-christophe.poggiale@univ-amu.fr

Pierre Auger
email: pierre.auger@ird.fr

The effect of connecting sites in the environment of a harvested population

Keywords: multisite fishery, Gordon-Schaefer model, time scales, maximum sustainable yield

This work presents a model of a harvested population in a multisite environment. Locally it has the shape of the Gordon-Schaefer model. This model gives rise, placing us in the case of a fishery, to an equilibrium of the stock and the fishing effort and, therefore, of the yield that is obtained per unit of time. Considering that the management of the fishery can act on the fishing costs, the yield is optimized as a function of the cost.

The objective of the work is to compare the maximum obtained yield in two extreme cases: unconnected sites and connected sites with rapid movements of both the stock and the fishing effort. The analysis of the model, first in an environment with two sites and later with any number of them, makes it possible to establish the conditions for one of the two cases to be more favorable from the point of view of the yield. In this way, it is proposed towards which of the two compared cases management should be directed.

Introduction

Multisite fisheries consist of a network of discrete fishing sites connected to each other by fish movements. Some sites can be created by man as areas where fishing is totally prohibited such as marine protected areas (MPA), [START_REF] Hilborn | Integrating marine protected areas with catch regulation[END_REF], or partially prohibited, especially at certain periods of the year. Other fishing areas are intentionally created as fish aggregating devices (FADs) because they attract fish that can be caught there more easily, [START_REF] Fonteneau | Observed changes in the species composition of tuna schools in the Gulf of Guinea between 1981 and 1999, in relation with the fish aggregrating device fishery[END_REF] and [START_REF] Dagorn | Behavior of yellowfin (Thunnus albacares) and bigeye (T. obesus) tuna in a network of fish aggregating devices (FADs)[END_REF]. The various fishing zones generally correspond to heterogeneous sites, in particular with regard to the environmental characteristics, such as the surface of the zone, the depth, the presence of a shelter zone, the underwater vegetation or even the resources available for the growth of the various species of commercial fish. All these local characteristics of the fishing sites affect the growth rate and carrying capacity of the fish species which is usually assumed to grow logistically. Heterogeneity can also result from human intervention, in particular by installing reefs or artificial habitats (AHs). AHs are known to increase the carrying capacity of the site where they are installed [START_REF] Polovina | Impacts of artificial reefs on fishery production in Shimamaki, Japan[END_REF] and attract fish [START_REF] Bortone | Resolving the attraction-production dilemma in artificial reef research: Some Yeas and Nays[END_REF]. In the case of a heterogeneous multisite fishery made up of sites that are more or less rich in the resources available for the fish, it is frequently assumed that the fish are distributed according to the resource available on each site. Fish will be more numerous in resource-rich areas than in poorer areas. For logistic sites, when the fish are distributed among them proportionally to the corresponding carrying capacities, they are said to adopt the Ideal Free Distribution (IFD), [START_REF] Bernstein | Predator Migration Decisions, the Ideal Free Distribution, and Predator-Prey Dynamics[END_REF] and [START_REF] Fretwell | On territorial behaviour and other factors influencing habitat distribution in birds. I. Theoretical development[END_REF].

In the case of a single fishing site, without adding economic assumptions, it is usual to consider the Schaefer model, [START_REF] Schaefer | Some considerations of population dynamics and economics in relation to the management of the commercial marine fisheries[END_REF][START_REF] Clark | Mathematical Bioeconomics. The Optimal Management of Renewable Resources[END_REF]:

N = rN 1 - N K -qEN (1)
It represents a population of size (stock) N (t) at time t, whose growth is depicted by the logistic model with intrinsic growth rate r > 0 and carrying capacity K. The exploitation term, qEN , follows the catch-per-unit-effort hypothesis and, thus, it is proportional to the 2

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d stock level, with E denoting the constant fishing effort and constant q being the catchability coefficient. If r -qE > 0, for any positive initial condition, N (t) tends to an equilibrium

N * = K(1 -qE/r).
It is also well known [START_REF] Clark | Mathematical Bioeconomics. The Optimal Management of Renewable Resources[END_REF] that the yield at equilibrium Y (E) = qK(1 -qE/r)E, as function of E, attains its maximum for E M SY = r/(2q). The corresponding yield, called maximum sustainable yield (MSY), is Y M SY = rK/4, and the equilibrium stock biomass is N M SY = K/2.

In [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF] a simple spatialized version of Schaefer's model is studied.

Only two sites are considered and the local dynamics is of the form (1). Furthermore, the fish stock move between the two sites on a faster time scale than that associated with local dynamics:

dN 1 dτ = m 2 N 2 -m 1 N 1 + ε r 1 N 1 (1 - N 1 K 1 ) -qE 1 N 1 dN 2 dτ = m 1 N 1 -m 2 N 2 + ε r 2 N 2 (1 - N 2 K 2 ) -qE 2 N 2 (2)
For i = 1, 2, N i is the population density on site i. Parameters r i , K i , and E i represent the intrinsic growth rate, the carrying capacity, and the fishing effort, respectively, on site i. The catchability is q. The movement rates for the stock, m i , are constant and site dependent. Time scales are included in the model by using the fast time variable τ together with the small positive parameter ε that represents their ratio.

As a consequence of the existence of two time scales, the two-dimensional model (2) can be analyzed through a one-dimensional Schaefer-type model. The MSY of (2) is easily approximated with the help of this one-dimensional model. From it, it is straightforward to compare the MSY for connected and isolated sites. In this latter case the MSY is calculated just by adding the two local MSY. An important result obtained in [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF] is that the MSY for a system of two connected fishing sites is always less than or equal to the MSY of the system with isolated sites. Thus, there is no interest in facilitating the connection between the sites because this would not increase the MSY of the system. This is the conclusion when we have not yet brought economic issues into the model.

Fishing boats also travel between the various fishing areas. The fishing boats are evenly distributed among the various fishing areas, avoiding having to meet all at the same time on the same site. The distribution of the boats in the fishing fleet on the different sites has an effect on the overall catch made on the multisite fishery. In a global way, the combined distribution of fish stock and fishing vessels among the various fishing areas has a significant effect on the overall catch of the multisite fishery, [START_REF] Brochier | Successful artificial reefs depend on getting the context right due to complex socio-bio-economic interactions[END_REF].

From the shipowner's point of view, it is important to distribute the fishing fleet over the different fishing areas to maximize his benefice which is constituted by the catch multiplied by the price of the resource on the market, which is usually called the landed value. The net income is obtained by subtracting from it the various operating costs of the fishery, such as the price of fuel necessary for the boats, the various taxes to be paid to the government, the wages of the fishermen, the maintenance of the vessels and again the minimum profit desired which corresponds to a threshold below which the shipowner considers that the fishing operation is not sufficiently profitable. The operating costs of the fishery therefore constitute an important parameter and it is legitimate to seek the cost per unit of fishing effort allowing the overall catch of a multisite fishery to be optimized.

To deal with this question it is necessary to introduce the economy in Schaefer's model.

A widespread way of doing this is to consider fishing effort as a system variable whose rate of change is proportional to the benefit, i.e., the difference between revenue and cost. We use Gordon's model of an open access fishery together with Schaefer's model to approach the analysis that we propose in this work. We want to study the conditions for a network of discrete fishing sites, connected to each other by fast displacements, to obtain a higher or lower fishing yield than if the sites were isolated.

We highlight the importance of the cost per unit of fishing effort by choosing it as the parameter to calculate the MSY. The Gordon-Schaefer model, if benefit is positive when the stock is at its carrying capacity, has a positive equilibrium to which all positive solutions tend. With the constant values that this equilibrium gives us for the stock and the fishing effort, the yield expression, that depends on all model parameters, is obtained.

The maximum of this expression as a function of cost c is what we consider to be the MSY.

We assume that the fishing sites are relatively close to each other allowing fast movements of the stock and the fishing fleet between these various fishing areas. When the sites 4

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d are close enough, it makes sense to consider that the operating costs are the same in all fishing areas, same wages, taxes, fuel and maintenance costs as well as minimum profit required. The government or the national fisheries management institution has no reason to impose different taxes on fishing sites close to each other.

Starting from the previous hypotheses, in this work, for a multisite fishery, we develop a comparison between the MSY with connected sites and with isolated sites. For it, in Section 2 we recall the basic Gordon-Schaefer model results including the optimization of the sustainable yield in terms of the cost. Section 3 considers a multisite fishery with two fishing sites and compares the MSY in the case where the sites are connected by movements of the stock and the fishing effort to the case where the sites are isolated in the absence of migrations. Section 4 extends the previous results to a multisite fishery with any number of sites. The manuscript ends with a discussion of the results, a conclusion and application perspectives.

A single-site harvested population model

As said in the introduction, the influence of economics on harvesting of renewable resources can be described with the help of the Gordon model of an open-access fishery, that is one in which anyone can harvest the resource. In its simplest form it is assumed a constant price p per unit of harvested biomass, so that the total revenue is obtained as pY (E), where Y (E) is the yield resulting from the effort E. The total cost is expressed as cE where c is a constant representing the cost per effort unit. Finally, the effort is considered variable with its rate of change being equal to the benefit pY (E)-cE. Equation

(1) together with the effort equation give the following simple model that it is at the base of all the models that we are presenting.

dN dt = rN 1 - N K -qEN dE dt = pqEN -cE (3)
We refer to (3) as the basic Gordon-Schaefer model. It has the same form as the Lotka-Volterra predator-prey model with prey logistic growth [START_REF] Iannelli | An Introduction to Mathematical Population Dynamics[END_REF]. The 5This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d asymptotic behaviour of its solutions is simply expressed. If c < pqK, i.e., the benefit is positive when the stock is at the carrying capacity, then all positive solutions of the system tend to the so-called bionomic equilibrium:

N * = c pq , E * = r q 1 - c pqK , (4)
with an associated yield

Y * = rc pq 1 - c pqK . (5)
In the open-access fishery, Gordon-Schaefer model, the fishery effort and the stock tend to reach an equilibrium at the level at which the benefit vanishes [START_REF] Clark | Mathematical Bioeconomics. The Optimal Management of Renewable Resources[END_REF], as the E equation imposes.

In the models that we are presenting, we consider the cost per unit of effort as a parameter that the manager can change, by providing the fishers with subventions to decrease it or by setting up some taxes to increase it. We thus look for a maximum sustainable yield considering the yield so far obtained as a function of the cost c.

In model (3), the value of c that maximizes the equilibrium yield ( 5) is

c opt = 1 2 pqK, (6)
obtaining

Y opt = 1 4 rK, with N opt = 1 2 K and E opt = r 2q .
So, we recuperate the situation of MSY (1) of the constant effort case.

Notice that c opt is one half of the maximum cost yielding positive benefit. It leads the stock N opt to one half of the environment carrying capacity K, that is the population size with the largest growth rate. Both, the optimal cost c opt and the stock N opt , as well as the optimal yield Y opt are proportional to K. Nevertheless, the corresponding effort E opt depends on the ratio intrinsic growth rate to catchability and it is independent of K. A consequence of this can be considered through the following example. A zone with carrying capacity K is managed in two different forms, as a whole and divided into two equal and 6

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d independent zones with carrying capacities K/2. In both cases the management follows

the dynamics of open-access fishery with chosen optimal cost. The global result is that at equilibrium the yield and the stock are equal in the two strategies, whereas the optimal cost is double in the first one and the effort is double in the second one.

A two-site harvested population model

We consider a harvested population in a multisite environment. The objective is showing that the maximum sustainable yield can be larger in a heterogeneous environment where stock and effort can move than it would be in the same environment without movements.

We begin by limiting the multisite environment to two sites. We will extend this for an arbitrary number of sites later on. In each site, the stock-effort dynamics is ruled by a local basic Gordon-Schaefer model ( 3). The movements between sites are described by constant rates and in a first approach, are considered fast in comparison to the local dynamics.

Thus, let us consider the following model:

dN 1 dτ = m 2 N 2 -m 1 N 1 + ε r 1 N 1 (1 - N 1 K 1 ) -qE 1 N 1 dN 2 dτ = m 1 N 1 -m 2 N 2 + ε r 2 N 2 (1 - N 2 K 2 ) -qE 2 N 2 dE 1 dτ = µ 2 E 2 -µ 1 E 1 + ε pqE 1 N 1 -cE 1 dE 2 dτ = µ 1 E 1 -µ 2 E 2 + ε pqE 2 N 2 -cE 2 (7)
where, for i = 1, 2, N i and E i are the population density and the fishing effort on site i. Parameters r i and K i represent the intrinsic growth rate and the carrying capacity, respectively, on site i. The catchability q, the price per harvest unit p and the cost per unit effort c are assumed equal in both sites. The movement rates for the stock, m i , and the effort, µ i , are constant and site dependent. Time scales are included in the model by using the fast time variable τ together with the positive parameter ε that represents their ratio.

We proceed now to compare the maximum sustainable yield in two different cases. In the first one, that we call connected sites case, stock and effort movements are allowed so that the dynamics is defined through system (7). In the second one, unconnected sites case, there is no stock and effort movements and, therefore, the dynamics is represented by the following independent systems, both with the form of (3), for i = 1, 2:

dN i dt = r i N i 1 - N i K i -qE i N i dE i dt = pqE i N i -cE i (8)
System ( 7) is not of the form of system (3) but it can be reduced to one of this form.

The fact that it includes two time scales allows us to apply the reduction method developed in (Auger et al. , 2008a,b). It consists in a sort of decoupling the fast and the slow parts of the system. The fast part of the system corresponds to the movements dynamics and let invariant the total stock N = N 1 + N 2 and the total fishing effort

E = E 1 + E 2 . Moreover,
it makes the proportions of stock and effort in each site rapidly tend to an equilibrium.

The equilibrium proportions for the stock are

u 1 = m 2 m 1 + m 2 and u 2 = m 1 m 1 + m 2 ,
and for the effort

v 1 = µ 2 µ 1 + µ 2 and v 2 = µ 1 µ 1 + µ 2 .
Assuming that movements have attained their equilibrium proportions, we can write a reduced system for the total stock and effort in terms of the slow time variable t = ετ . We add the two stock equations and the two effort equations, and substitute each of the state variables by the corresponding total variable times the associated site proportion, i.e., for

i = 1, 2, N i = N u i and E i = Ev i .
The result is the following basic Gordon-Schaefer model:

dN dt = rN (1 - N K ) -qEN dE dt = pqEN -cE (9) whose parameters, r = r 1 u 1 +r 2 u 2 , K = (r 1 u 1 + r 2 u 2 )K 1 K 2 r 2 u 2 2 K 1 + r 1 u 2 1 K 2 and q = q(u 1 v 1 +u 2 v 2 ), include
the local parameters and the equilibrium stock and effort proportions on sites. The so-8

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d called aggregated model ( 9) possesses a positive globally asymptotically stable equilibrium

provided that K > c pq . It is given by ( N , Ē) = c pq , r q 1 -N K .
The yield at equilibrium is thus

Y = rc pq (1 - c pq K ).
Thanks to the aggregated model, all these calculations are very simple and provide us with results concerning the complete model. Figure ( 1) illustrates that the above mentioned results work perfectly.

To compare the maximum sustainable yield in the connected and unconnected sites cases, we now consider the yield at equilibrium maximized as a function of a single cost.

In the first case using system (9) and in the second case with the help of systems (8).

Connected sites

The maximum sustainable yield as a function of the cost, in the case with fast movements of the stock and the efforts can be obtained by using (6) in system (9). For the cost

c opt = 1 2 pq K = 1 2 pq(u 1 v 1 + u 2 v 2 )(r 1 u 1 + r 2 u 2 )K 1 K 2 r 2 u 2 2 K 1 + r 1 u 2 1 K 2 the maximum yield is Y opt = 1 4 r K = 1 4 (r 1 u 1 + r 2 u 2 ) 2 K 1 K 2 r 2 u 2 2 K 1 + r 1 u 2 1 K 2 . ( 10
)
This expression does not depend on either the price, the catchability or the efforts movements rates. On the other hand, c opt does depend on all of these quantities. The stock and effort at equilibrium are

N opt = 1 2 K = 1 2 (r 1 u 1 + r 2 u 2 )K 1 K 2 K 1 r 2 u 2 2 + K 2 r 1 u 2 1 , E opt = r 2q = (r 1 u 1 + r 2 u 2 ) 2q(u 1 v 1 + u 2 v 2 )
.

9

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

m 1 = 2, m 2 = 1, µ 1 = 2, µ 2 = 3, r 1 = 1, r 2 = 2, K 1 = 5, K 2 = 5, q = 1, p = 1, and c = 1 . ε = 0.1 is chosen.
Top panel shows the total population density simulated with the complete (blue) and aggregated (dashed red) models. The centre panel shows the effort obtained with the complete (blue) and aggregated (dashed red) models. The bottom panel corresponds to the simulation of the yield. In blue and dashed red respectively are the yields obtained with the complete and aggregated models with the same parameter values as previously. In orange is the yield obtained with the complete model for a cost chosen in a such way that the yield becomes maximal at equilibrium. The dashed black horizontal line is the yield level obtained from the relation Y opt = rK/4.

The maximum yield (10) is strongly dependent on the equilibrium distribution of the stock to which the movements give rise. Let us show this with the help of a simple case. If we assume two biologically identical sites, r 1 = r 2 = r and K 1 = K 2 = K/2, then we have

Y opt = rK 4 • 1 2(u 2 1 + u 2 2 ) If the stock is homogeneously distributed, u 1 = 1/2 = u 2 , then Y opt = rK/4, that corre- 10
This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d sponds to the MSY considering a single stock with biological parameters r and K. Nevertheless, if u 1 approaches 1 and, thus, u 2 approaches 0 (or vice versa), i.e., the stock is mainly established in the first site (or the second one), then Y opt tends to rK/8. Uneven stock distribution could almost halve the maximum yield.

Let us see which is the equilibrium stock distribution that maximizes Y opt in a general situation. For this, we consider Y opt as a function of u 1 , making u 2 = 1 -u 1 :

Y opt (u 1 ) = 1 4 (r 1 u 1 + r 2 (1 -u 1 )) 2 K 1 K 2 r 2 (1 -u 1 ) 2 K 1 + r 1 u 2 1 K 2 .
Its derivative can be expressed as

Y opt (u 1 ) = r 1 r 2 K 1 K 2 r 1 u 1 + r 2 (1 -u 1 ) 2 r 2 (1 -u 1 ) 2 K 1 + r 1 u 2 1 K 2 2 K 1 -(K 1 + K 2 )u 1 ,
and, therefore, its sign just depends on expression

K 1 -(K 1 + K 2 )u 1 . The maximum is attained for u 1 = K 1 K 1 + K 2 , that implies u 2 = K 2 K 1 + K 2
. This stock distribution, proportional to the sites carrying capacities, is known in ecological literature as the Ideal Free Distribution (IFD) [START_REF] Fretwell | On territorial behaviour and other factors influencing habitat distribution in birds. I. Theoretical development[END_REF]. The corresponding yield is

Y opt K 1 K 1 + K 2 = 1 4 r 1 K 1 + r 2 K 2 , that it is associated to a cost c opt = 1 2 pq(K 1 v 1 + K 2 v 2 ) at an equilibrium N opt = 1 2 (K 1 + K 2 ) and E opt = 1 2 r 1 K 1 + r 2 K 2 q(K 1 v 1 + K 2 v 2 )
.

The yield Y opt (u 1 ) grows from 1 4 r 2 K 2 , for u 1 = 0 (the whole stock living in site 2), to its maximum, at the IFD, and then decreases to 1 4 r 1 K 1 , for u 1 = 0 (the whole stock living in site 1).

The closer the stock distribution is to the IFD, the larger is the MSY obtained.

In the unconnected case, this question does not apply because, as there are no stock movements, its distribution is constant.

11

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

Unconnected sites

Now, let us assume that there is no movement. The sites are independent of each other from the point of view of stock and effort, though they share the catchability q, the price p, and the cost c.

In this case, the yield of the system is the sum of the yields on each separate site. Using expression [START_REF] Auger | Aggregation methods in dynamical systems and applications in population and community dynamics[END_REF] in both systems (8), one gets:

Y S = r 1 c pq 1 - c pqK 1 + r 2 c pq 1 - c pqK 2
As a function of c, Y S has a maximum for the cost c S,opt = (r

1 + r 2 )pqK 1 K 2 2(r 1 K 2 + r 2 K 1 ) equal to Y S,opt = (r 1 + r 2 ) 2 K 1 K 2 4(r 1 K 2 + r 2 K 1 ) . (11)
The equilibria attained by the stock are the same in both sites:

N 1 S,opt = c S,opt pq = (r 1 + r 2 )K 1 K 2 2(r 1 K 2 + r 2 K 1 ) = N 2 S,opt .
Due to different biological parameters in the two sites, the equilibrium attained by the effort in each site is different

E 1 S,opt = r 1 q 1 - c S,opt pqK 1 = r 1 q 1 - (r 1 + r 2 )K 2 2(r 1 K 2 + r 2 K 1 ) , E 2 S,opt = r 2 q 1 - (r 1 + r 2 )K 1 2(r 1 K 2 + r 2 K 1 )
.

In the previous case of connected sites we showed that the MSY ( 10) is strongly dependent on equilibrium stock distribution. We found that to maximize MSY this distribution must be the IFD. In this case of unconnected sites the equilibrium stock distribution does not apply because they are independent. Nevertheless, we are showing that how the total carrying capacity of the system, K = K 1 + K 2 , is distributed between the sites does play a relevant role.

We start approaching the question in a simple case. Let us assume r 1 = r 2 = r, and represent the total carrying capacity distribution by means of the fraction α ∈ (0, 1) of K

12

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d in the first site, i.e., K 1 = αK and K 2 = (1 -α)K. Substituting these values in (11) we get

Y S,opt = α(1 -α)rK,
which attains the well known maximum 1 4 rK for α = 1 2 , that is, when there is half of the total carrying capacity in each site. On the other hand, if K tends to be mostly in one of the sites the MSY approaches zero.

In the general case, r 1 = r 2 , the MSY as a function of α is

Y S,opt (α) = (r 1 + r 2 ) 2 Kα(1 -α) 4(r 1 (1 -α) + r 2 α) .
This function also tends to zero when α tends to zero or one. Its maximum is attained for

α = √ r 1 √ r 1 + √ r 2 ,
that is, when the carrying capacities of the sites are proportional to the square roots of the corresponding intrinsic growth rates. If this is the case, the maximum yield is

Y S,opt √ r 1 √ r 1 + √ r 2 = 1 4 r 1 + r 2 √ r 1 + √ r 2 2 K.
where

r 1 + r 2 √ r 1 + √ r 2 2
is in between r 1 and r 2 .

Uniform costs in separated sites seem reasonable if the carrying capacities ratio is close to the ratio of the square roots of the intrinsic growth rates.

Comparison of the MSY for connected and unconnected sites

To compare the MSY obtained in the two cases, connected and unconnected sites, we define the ratio ρ Y of Y opt (10) to Y S,opt (11),

ρ Y = Y opt Y S,opt = (r 1 K 2 + r 2 K 1 )(r 1 u 1 + r 2 u 2 ) 2 (r 1 + r 2 ) 2 (r 2 u 2 2 K 1 + r 1 u 2 1 K 2 ) . ( 12
)
If ρ Y > 1 then we can conclude that the movements increase the MSY whereas ρ Y < 1 indicates that the movements reduce the MSY.

13

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482

P r e p r i n t n o t p e e r r e v i e w e d

We can remark that ρ Y depends only on the biological parameters, r i and K i , and the equilibrium stock distribution given by frequencies u 1 and u 2 = 1 -u 1 . On the other hand, it does not depend on the efforts frequencies v i , that is, changing the efforts movements rates will not change ρ Y .

As we pointed out in Section 3.1, the MSY in the case of connected sites attains its maximum when the equilibrium stock distribution coincides with IFD. The value of ρ Y in this case (u 1 = K 1 /(K 1 + K 2 )) can be expressed in the following form:

ρ IF D Y = (r 1 + r 2 ) 2 K 1 K 2 + r 1 r 2 (K 1 -K 2 ) 2 (r 1 + r 2 ) 2 K 1 K 2 , therefore, ρ IF D Y > 1 except in the case of equal carrying capacities, K 1 = K 2 , in which it is equal to 1.
If the stock movements tend to IFD, the system with connected sites gives a larger MSY than the one with unconnected sites. If IFD is not the equilibrium of stock movements then this is not always the case.

Considering ρ Y as a function of u 1 it is straightforward to find the conditions for ρ Y > 1 to hold. Solving the equation ρ Y (u 1 ) = 1 we obtain two solutions 1 2 and ū1 :=

(r 1 + 2r 2 )K 1 -r 2 K 2 2(r 2 K 1 + r 1 K 2 ) .
Depending on K 1 > K 2 , or K 1 < K 2 , the second root satisfies ū1 > 1/2, or ū1 < 1/2, and we have then ρ Y > 1 if and only if one of the next two conditions is met

K 1 > K 2 and u 1 ∈ 1 2 , min{ū 1 , 1} (13)
K 1 < K 2 and u 1 ∈ max{ū 1 , 0}, 1 2 (14)
We illustrate the results obtained with Figure 2.

To do it, we first express ρ Y in terms of the ratio between K 1 and K 2 ,

ρ K = K 1 /K 2 ,
and the ratio between r 1 and r 2 , ρ r = r 1 /r 2 , getting If both sites have equal carrying capacities, K 1 = K 2 , i.e., ρ K = 1, we have

ρ Y = ρ r u 1 + u 2 ρ r + 1 2 ρ r + ρ K ρ r u 2 1 + ρ K u 2 2 . ( 15
) (a) ρ K = 1 (b) ρ K = 2 (c) ρ K = 5 (d) ρ K = 10
263 ρ Y = ρ r u 1 + u 2 ρ r + 1 2 ρ r + 1 ρ r u 2 1 + u 2 2 = (ρ r u 1 + u 2 ) 2 (ρ r u 1 + u 2 ) 2 + ρ r (u 1 -u 2 ) 2 ,
then ρ Y = 1 for u 1 = 1/2 = u 2 and is less than 1 for any other value of u 1 . Figure 2a shows 264 the value of ρ Y as a function of u 1 (u 2 = 1 -u 1 ) for ρ K = 1 and different values of ρ r .

265

The rest of the figures in Figure 2 make this same representation for three increasing 266 values of ρ K . In them you can clearly see the range of values of u 1 (13), the fraction of 267 the stock that stabilizes in site 1, for which ρ Y > 1, i.e., for which the optimal yield is higher if the sites are connected. It can be seen that in all cases the maximum of ρ Y is approximately reached when u 1 represents the IFD, that is, for u 1 = ρ K /(1 + ρ K ). We note that to obtain ρ Y > 1 it is necessary for the equilibrium stock distribution to have more weight in the site with the highest carrying capacity. In Figure 2 we see that for any u 1 < 1/2 the site connection does not improve the maximum yield, i.e., ρ Y < 1. Concerning the role of the intrinsic growth ratio, ρ r , we observe that the maximum value of ρ Y , for fixed ρ K , is obtained for ρ r = 1, that is, when r 1 = r 2 .

An L-site harvested population model (L > 2)

In this section we extend the analysis carried out in section 3 to a general multisite environment encompassing L different sites.

We continue to consider a harvested population and are interested in comparing the maximum sustainable yield between the case with connected sites and the case with isolated sites.

In each site, the stock-effort dynamics is ruled by a local basic Gordon-Schaefer model (3).

In the case of connected sites, we assume that movements between sites are fast compared to local dynamics. They are described by constant rates. Let us denote m ij and µ ij , i = j and i, j ∈ {1, . . . , L}, the transition rates from site j to site i, for stock and effort respectively. To lighten the model writing, let us define, for i ∈ {1, . . . , L},

m ii = - L j = 1 j = i m ji , µ ii = - L j = 1 j = i µ ji .
They represent the rate at which stock or effort leave site i to go to the rest of sites.

16

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482

P r e p r i n t n o t p e e r r e v i e w e d

The extension of model ( 7) to an L-site environment reads as follows:

dN i dτ = L j=1 m ij N j + ε r i N i 1 - N i K i -qE i N i , dE i dτ = L j=1 µ ij E j + ε pqE i N i -cE i , i = 1, . . . , L. (16)
where, N i and E i are the population density and the fishing effort on site i. Parameters r i and K i represent the intrinsic growth rate and the carrying capacity, respectively, on site i. The catchability q, the price per harvest unit p and the cost per unit effort c are still assumed to be equal in all sites. The fast time variable is τ and positive parameter ε 1 represents the time scales ratio.

In the case of unconnected sites, as there is no stock and effort movements, the dynamics is represented by the following independent systems, for i = 1, . . . , L:

dN i dt = r i N i 1 - N i K i -qE i N i dE i dt = pqE i N i -cE i (17)
Assuming that matrices M 1 = (m ij ) and M 2 = (µ ij ) are irreducible, the same method we used to reduce system (7), (Auger et al. , 2008a,b), apply to system ( 16). The interpretation of the irreducibility of a matrix in this context is that the connections between sites allow to establish a path joining any site to any other. The movements dynamics let invariant the total stock

N = N 1 + • • • + N L and the total fishing effort E = E 1 + • • • + E L .
Moreover, it makes the proportions of stock and effort in each site rapidly tend to an equilibrium. Let ū = (u 1 , . . . , u L ) and v = (v 1 , . . . , v L ) be these equilibrium proportions for stock and effort respectively. Vectors ū and v are the right eigenvectors of matrices M 1 and M 2 associated to eigenvalue 0 and whose entries sum up to 1, see Section A.5. in Smith and Thieme (2011) for details on this kind of matrices called quasipositive irreducible.

The associated reduced system for the total stock and effort in terms of the slow time variable t = ετ is the following basic Gordon-Schaefer model:

dN dt = L i=1 r i u i N 1 - u i N K i -qv i Eu i N = rN 1 - N K -qEN dE dt = L i=1 pqv i Eu i N -cv i E i = pqEN -cE (18)
where

r = L i=1 u i r i , K = r L i=1 (u 2 i r i K i ) , q = L i=1 (u i v i ) q. ( 19
)
To compare the maximum sustainable yield in the cases of connected and unconnected sites, we now consider the yield at equilibrium maximized as a function of a single cost. In the first case using system (18) and in the second case with the help of systems (17).

Connected sites

In the case of connected sites, we apply ( 6), the value of the cost that maximizes the equilibrium yield in the basic Gordon-Schaefer model, to system (18), obtaining

c opt = 1 2 pq K,
and the maximum equilibrium yield

Y opt = 1 4 r K = 1 4 L i=1 u i r i 2 L i=1 (u 2 i r i /K i ) . ( 20
)
It does not depend on either the price, the catchability or the efforts movements rates. On the other hand, it is strongly dependent on the equilibrium stock distribution. Let us show this through a simple extreme case as we did in the case of two sites. Assuming biologically identical sites, r i = r and K i = K/L, for i = 1, . . . , L, then we have

Y opt = rK 4 • 1 L L i=1 u 2 i .
18

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482

P r e p r i n t n o t p e e r r e v i e w e d

If the stock is homogeneously distributed, u i = 1/L, then Y opt = rK/4, that corresponds to the MSY considering a single stock with biological parameters r and K. Nevertheless, if the stock tends to concentrate in a single site, let us think that u 1 approaches 1 and u i approaches 0 for i = 2, . . . , L, then Y opt tends to rK/(4L), i.e., the previous maximum yield is divided by the number of sites.

If we now consider Y opt as a function of variables u i (i = 1, . . . , L) we can prove that the maximum of Y opt is obtained for a distribution that corresponds to IFD, i.e., u i =

K i / L j=1 K j . To see it, let us consider function Y opt (u 1 , . . . , u L ) with the constraint u 1 + • • • + u L = 1, u i > 0,
and use a Lagrange multiplier λ. For j = 1, . . . , L,

λ = ∂Y opt ∂u j = 1 4 r K = 1 2 r j K - 1 2 K2 u j r j K j ,
multiplying each equality by u j and summing up all of them yields

λ = L j=1 u j 1 2 r j K - 1 2 K2 u j r j K j = 1 2 r K - 1 2 K2 r K = 0,
substituting λ = 0 in the first equality gives

u j = K j K ,
and, as

u 1 + • • • + u L = 1, K = K 1 + • • • + K L .
So, as we wanted to prove, the stock distribution that maximizes Y opt is the IFD. This yield value is

Y IF D opt = Y opt (K 1 / K, . . . , K L / K) = 1 4 L i=1 r i K i ,
and it is associated to a cost c opt = 1 2 pq L i=1 r i K i at an equilibrium

N opt = 1 2 L i=1 K i and E opt = L i=1 r i K i / 2q L i=1 v i K i .
19

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

Unconnected sites

When the sites are not connected we must consider them independently. Thus, the yield of the system is just the sum of the yields on each separate site. Using expression [START_REF] Auger | Aggregation methods in dynamical systems and applications in population and community dynamics[END_REF] in all the L systems (17), one gets:

Y S = L i=1 r i c pq 1 - c pqK i = c pq L i=1 r i - c 2 p 2 q 2 L i=1 r i K i
We look for a common value of c that optimizes Y S . Considering Y S as a function of c, the maximum is attained at

c S,opt = 1 2 pq L i=1 r i / L i=1 r i /K i ,
and the maximum yield is equal to

Y S,opt = L i=1 r i 2 / 4 L i=1 r i K i . (21)
The equilibria attained by the stock is the same in every site j = 1, . . . , L

N j S,opt = 1 2 L i=1 r i / L i=1 r i /K i ,
but this need not be the case concerning the effort equilibria

E j S,opt = r j q 1 - L i=1 r i / 2K j L i=1 r i /K i .
In the previous case of connected sites we showed that the MSY ( 20) is strongly dependent on equilibrium stock distribution associated to stock fast movements. When this distribution is the IFD the maximum equilibrium yield attains its maximum. In the case of unconnected sites movements are not considered and, therefore, it does not make sense to speak about the equilibrium stock distribution they lead to. Nevertheless, we can consider the distribution among sites of the total carrying capacity of the system, K = L i=1 K i , and show that it does play a relevant role by treating a simple case.

20

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

Assuming, equal intrinsic growth rates in all sites, r i = r for i = 1, . . . , L, and setting

K i = α i K, with α i ∈ (0, 1), we have Y S,opt = L i=1 r 2 / 4 L i=1 r α i K = 1 4 rKL L L i=1 1/α i ,
where the term L/ L i=1 1/α i is the harmonic mean of the constants α i . As the harmonic mean is always less than or equal to the arithmetic mean we have that

Y S,opt ≤ 1 4 rK,
with the equality obtained if all sites share the same carrying capacity, i.e., α i = 1/L for all i = 1, . . . , L. On the other hand, it can be noted that if one of the constants α i tends to zero then the harmonic mean of all of them too, what entails the Y S,opt also tending to zero.

Uniform cost in separated sites with very different carrying capacities would lead to small Y S,opt .

If we treat the general case with possibly different local intrinsic growth rates r i , we obtain analogous results to those found in the model with two sites. Let us set, as above, K i = α i K, with α i ∈ (0, 1), and consider Y S,opt as a function of variables α i (i = 1, . . . , L).

To calculate its maximum we can use function Y S,opt (α 1 , . . . , α L ) with the constraint α 1 +

• • • + α L = 1, α i > 0, and use a Lagrange multiplier λ. For j = 1, . . . , L,

λ = ∂Y S,opt ∂α j = K 4 • L i=1 r i 2 L i=1 (r i /α i ) 2 • r j α 2 j , which implies that α j is proportional to √ r j α j = 1 2 K λ • L i=1 r i L i=1 (r i /α i ) √ r j
and, therefore, for j = 1, . . . , L,

α j = √ r j L i=1
√ r i .

If the carrying capacities K i of the sites are proportional to the square roots of the corresponding intrinsic growth rates √ r i then Y S,opt attains its maximum

Y S,opt √ r 1 L i=1 √ r i , • • • , √ r L L i=1 √ r i = 1 4 K L i=1 r i L i=1 √ r i 2 .

Comparison of Y opt and Y S,opt

To compare the maximum yield as obtained in the two studied cases, connected and unconnected sites, we define the ratio ρ Y of Y opt (20) to Y S,opt (21),

ρ Y = Y opt Y S,opt = L i=1 u i r i 2 L i=1 r i /K i L i=1 r i 2 L i=1 u 2 i r i /K i . ( 22
)
If ρ Y > 1, fast movements increase the maximum yield whereas ρ Y < 1 indicates that movements reduce it.

As it was remarked in the model with two sites, ρ Y depends only on the biological parameters, r i and K i , and on the stock movements through the equilibrium stock distribution u i , but it does not depend on the efforts frequencies v i . Changing the efforts movements rates will not change the value of ρ Y .

As we pointed out in Section 4.1, the maximum yield in the case of connected sites is obtained when the equilibrium stock distribution coincides with IFD. The value of ρ Y in this case (u i = K i / L i=1 K i , i = 1, . . . , L) is the following:

ρ IF D Y = L i=1 (r i K i ) L i=1 (r i /K i ) L i=1 r i 2 .
Using that function f (x) = x + 1/x satisfies f (x) > 2 for x ∈ (0, 1) ∪ (1, ∞) and f (1) = 2, we have, for any i, j = 1, . . . , L,

K i K j r i r j + K j K i r j r i > r i r j + r j r i if K i = K j , (23)
and the equality of both expressions holds if

K i = K j . Now,
K i K j r i r j + K j K i r j r i ) + L i=1 r 2 i L i=1 L j=i+1 (r i r j + r j r i ) + L i=1 r 2 i ,
we obtain that ρ IF D Y > 1 except in the case of equal carrying capacities, K i = K j for all i, j ∈ {1, . . . , L}, in which it is equal to 1.

If the stock movements tend to IFD, the system with connected sites gives a larger maximum yield than the one obtained with unconnected sites. If IFD is not the equilibrium of stock movements then this is not always the case. It does not seem straightforward to find necessary and sufficient conditions for ρ Y being bigger or less than 1. To show that unconnected sites can also lead to a larger maximum yield, let us calculate ρ Y in the simple case with equal intrinsic growth rates and equal carrying capacities (ECC) in all sites. Let r i = r and K i = K/L for i = 1, . . . , L, where K represents the total carrying capacity.

The ρ Y , in this case, is readily calculated to be

ρ ECC Y = 1 L L i=1 u 2 i .
As u i ≥ 0 and L i=1 u i = 1, it is straightforward to find that 1/L ≤ L i=1 u 2 i ≤ 1 and, thus,

1 L ≤ ρ ECC Y ≤ 1.
In this case, equal intrinsic growth rates and equal carrying capacities, the system with unconnected sites gives a larger maximum yield than the one obtained with connected sites. If stock movements tend to IFD, then ρ Y approaches 1, but it can be much less than 1 if stock movements equilibrium distribution gets away from IFD.

Conclusion

In the early contribution [START_REF] Freedman | Mathematical models of population interactions with dispersal. I. Stability of two habitats with and without a predator[END_REF], it was shown that the global carrying capacity of two sites connected by fast migrations where each sub-population grows logistically could be globally greater than the sum of the carrying capacities of the two isolated sites. This result has received a lot of attention recently and has been discussed in several other articles, [START_REF] Holt | Population Dynamics in Two-Patch Environments: Some Anomalous Consequences of an Optimal Habitat[END_REF][START_REF] Poggiale | Global production increased by spatial heterogeneity in a population dynamics model[END_REF][START_REF] Deangelis | Effects of dispersal in a non-uniform environment on population dynamics and competition: a patch model approach[END_REF][START_REF] Deangelis | Effects of diffusion on total biomass in heterogeneous continuous and discrete-patch systems[END_REF][START_REF] Arditi | Is dispersal always beneficial to carrying capacity? New insights from the multi-patch logistic equation[END_REF][START_REF] Arditi | Asymmetric dispersal in the multi-patch logistic equation[END_REF][START_REF] Zhang | Carrying Capacity of Spatially Distributed Metapopulations[END_REF]. A condition for obtaining higher productivity from the system of two connected sites compared to the system of two isolated sites is to consider a heterogeneous environment, that is to say with two sites having different environmental properties, in particular having different growth rates and carrying capacities, [START_REF] Poggiale | Global production increased by spatial heterogeneity in a population dynamics model[END_REF].

In this work we have analyzed a similar problem. We have treated a population that grows according to the logistic model but that, in addition, is harvested. The actual system we have in mind is a fishery. The environment of the population can be considered divided into different sites, that is, we are dealing with the case of a multisite fishery. Finally, what we seek to compare between the system of connected and unconnected sites is not the total carrying capacity of the system, as in the works mentioned above, but the maximum sustainable yield.

This problem, as we mentioned in the introduction, has already been dealt with in [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF]. As can be seen in system (2), the population is harvested at a rate proportional to its size following the classical Schaefer's model. Yield optimization is performed with respect to the fishing effort parameter E. The analysis of the model offers a result that we could describe as negative. Regardless of the movement rates chosen and whether or not the sites are environmentally homogeneous, connecting the sites does not allow increasing the sum of the MSY of the unconnected sites. Equality occurs when the exploited species is distributed over the two sites according to the ideal free distribution.

This result is quite general as soon as the harvested population equation includes a logistic growth term from which the catch is subtracted. Indeed, at equilibrium, the capture is equal to a logistic growth term which is always maximum when the population is at half the carrying capacity. This always leads to the same value of the catch at MSY and, therefore, to the inability to increase optimal yield by connecting sites. A positive answer is obtained in the same work [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF] considering a prey-predator community of fish in the same environment, where only the predators are harvested. It is shown that the total MSY of the system can be greater with connected sites. The condition leading to this result consist in connecting a fishing site with a large prey carrying capacity and an average growth rate to a site with a small prey carrying capacity but a large growth rate.

The models analyzed in this paper includes economic factors, which were not considered in [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF]. The fishing effort E becomes a variable whose dynamics follow the hypothesis of the open access fishery. The non-spatialized model, (3), leads, under a certain condition, to an equilibrium of both the stock and the effort, and therefore also the yield. From this, we have chosen the cost c as a parameter with respect to which to optimize the yield. We have preferred c over other possible choices because it includes all the costs of operating a fishery, including the profit margin that the owner considers profitable. In addition, it is an element on which management can act through taxes or subsidies.

The model ( 3) is used as a local model to build a multisite fishery model. The management of the fishery is unique and the same cost and price are established in each of the sites. When unconnected sites are considered, the model for the fishery is made up of a collection of independent systems (3). When the sites are connected, with rapid movements of the stock and of the fishing effort, the model ( 7) is proposed for two sites and, in general, for any number L of sites, the model ( 16). Once the maximum sustainable yield (MSY) is obtained based on the cost c, the comparison between the fishery with connected and unconnected sites based on this criterion is developed. The ratio ρ Y between the MSY with connected sites, Y opt , and with unconnected sites, Y S,opt , is used to decide whether to prefer a connected fishery, if ρ Y > 1, or an unconnected fishery, whenever ρ Y < 1.

In the case of two sites the ratio ρ Y depends on two other ratios, that of the carrying capacities, ρ K , and that of the intrinsic growth rates, ρ r , and on the equilibrium distribution of the stock reflected in u 1 . The rest of the parameters do not play any role. A first general conclusion is that if the carrying capacities of the two sites are equal, ρ K = 1, then for any values of u 1 and ρ r we have that ρ Y ≤ 1, with the equality valid only if u 1 = 1/2. If the two sites have the same carrying capacities, it is better to keep the sites unconnected.

Another general result but in the opposite sense follows. If the equilibrium distribution of the stock coincides with the IFD then it is better to use connected sites. Figure 2 shows that whenever the carrying capacity of site 1 is greater than that of site 2, ρ K > 1, there is an interval of values of u 1 , whose left endpoint is 1/2 and, on the right, can be as high as 1, such that ρ Y > 1. If the fishery can be managed so that the equilibrium distribution of the stock is in this interval, connecting the sites will give a higher MSY. Note that the best results are obtained when the equilibrium distribution of the stock is around the IFD.

The values of ρ Y obtained grow with ρ K and are reached for values of ρ r close to 1, i.e. equal intrinsic growth rates. A symmetric analysis is valid if the largest carrying capacity is that of site 2. Similar general conclusions can be reached in the case of L sites.

This work shows the importance of maintaining corridors allowing individuals to move from one site to another in a network of sites. It also shows that the heterogeneity between the different sites is a favorable factor for the increase in the total optimal catch in the multisite fishery. This work makes it possible to specify the conditions favorable to this increase in the capture at the MSY. The manager of the multisite fishery can use certain devices such as artificial reefs to increase the total productivity of the multisite fishery.

In the case of two connected fishing areas, it would be possible to increase the carrying capacity in only one of the two zones by installing artificial reefs there, which are known to increase the productivity of the halieutic resource by increasing the local carrying capacity. Artificial reefs have a priori no effect on the growth rate of the exploited species.

Thus it appears possible to create an asymmetry between the fishing zones by increasing the ratio of local carrying capacities ρ K while maintaining the ratio of growth rates ρ r at one. Figure 2 shows that it is possible as soon as ρ K is equal to 2 or more to obtain ρ Y ≥ 1.

In [START_REF] Zhang | Carrying capacity in a heterogeneous environment with habitat connectivity[END_REF], the authors performed experiments with aquatic plants growing logistically in containers. They simulated in real world the fast migrations by transferring plants from one container to another with a periodicity of a few days. They considered a heterogeneous environment with different growth conditions from container to container of different sizes. Under these conditions, they showed that the total productivity, in terms of carrying capacities, could be greater in heterogeneous conditions rather than in homogeneous conditions. An important perspective of this work is to consider carrying out experiments of the same type making it possible to generalize the results obtained in [START_REF] Zhang | Carrying capacity in a heterogeneous environment with habitat connectivity[END_REF] to the case of exploited species, either a single commercial species as in this work with an identical operating cost in all sites or in the case of fish prey-predator systems with exploitation of the predator, [START_REF] Auger | Increase of maximum sustainable yield for fishery in two patches with fast migration[END_REF].

Figure 1 :

 1 Figure 1: Comparison between the complete model (7) and the aggregated model (9) for the following parameters values:m 1 = 2, m 2 = 1, µ 1 = 2, µ 2 = 3, r 1 = 1, r 2 = 2, K 1 = 5, K 2 = 5, q = 1, p = 1, and

Figure 2 :

 2 Figure 2: Graphs of ρ Y as a function of u 1 for different values of ρ K : 1, 2, 5, and 10; and different values of ρ r : 0.1, 0.5, 1, 2, and 10.

  as we can write ρ IF D Y in the following form 22This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=4176482 P r e p r i n t n o t p e e r r e v i e w e d

Acknowledgements

Rafael Bravo de la Parra is supported by Ministerio de Universidades (Spain) with a grant in the Programa Estatal de Promoción del Talento y su Empleabilidad en I+D+i, Subprograma Estatal de Movilidad, of the Plan Estatal de Investigación Científica y Técnica y de Innovación 2017-2020.


Mladen Čičić

Carlos Vivas

Carlos Canudas De Wit

Francisco Rubio

Francisco Rubio Optimal

Optimal Renewable Energy Curtailment Minimization Control Using a Combined Electromobility and Grid Model

Keywords: Electromobility, Optimal control, Renewable energy sources, Energy curtailment, Grid optimization

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Optimal Renewable Energy Curtailment Minimization Control Using a Combined Electromobility and Grid Model

Mladen Čičić * Carlos Vivas * * Carlos Canudas-de-Wit * Francisco R. Rubio * *

INTRODUCTION

As efforts towards decarbonization of all economical sectors become a major priority, Electric Vehicles (EVs) have started to emerge as one of the main components of sustainable transportation systems worldwide. Since EVs are projected to reach around 40% of the total fleet in the EU by 2030 [START_REF] Conway | A review of current and future powertrain technologies and trends in 2020[END_REF], it is becoming clear that their integration with the city infrastructure (charging stations), and the electrical power supply network (power grid) poses unsolved problems that will be critical in the coming years (EU 2019).

Although the constant increase of the electrification of the transportation systems (electromobility) could be interpreted as a serious potential strain on the power grid due to the large charging power demands [START_REF] Fernandez | Assessment of the impact of plug-in electric vehicles on distribution networks[END_REF], the massive adoption of EVs will not necessarily hinder the development of future electric power systems.

In fact, the potential to use their batteries for energy storage will help with the introduction of renewable energy sources (RES) [START_REF] Wenzel | Real-time charging strategies for an electric vehicle aggregator to provide ancillary services[END_REF] by increasing their dispatchability, making the EVs an important component of the power system operations. From the power system side, the grid can use this capacity to deal with network congestion management issues, as they can provide EVflexibility (ability to quickly modify the demand to react to changing conditions), which is especially important in a foreseen future energy scenario dominated by intermittent RES. As the natural interface between the mobility and ⋆ E-mail: mladen.cicic@gipsa-lab.fr, carlos.canudas-de-wit@cnrs.fr, vivas@us.es, rubio@us.es power networks, the charging stations will play an essential role in the electromobility ecosystem. Today's power electronics technology and new DC grid topologies, together with V2G-enabled EVs, allow charging stations to operate as dynamic load/supplies providing ancillary services to the power grid in the form of frequency stabilization and congestion relief [START_REF] López | V2g strategies for congestion management in microgrids with high penetration of electric vehicles[END_REF]. Note that the EV flows on the roads can also be seen as additional virtual links in the power grid, which not only store energy but also carry energy from one point to another (a 3-lane highway with 25% of the traffic being EVs whose batteries have a capacity of 40 kWh is approximately equivalent to a 110 kV power line in terms of transported power!).

One of the main potential barriers to fully exploit the EVs' potential is the lack of tools and methods for forecasting EV fleets' flexibility in both time and space. This entails forecasting when and where EVs will move, how their State of Charge (SoC) will evolve, and how they will interact with the infrastructure. Though some approaches based on historical data do exist [START_REF] Morlock | Forecasts of electric vehicle energy consumption based on characteristic speed profiles and real-time traffic data[END_REF], a modelbased framework is preferable for optimal control purposes. Therefore, combining electromobility models with power grid models will enable the use of EVs' flexibility potential to minimize the RES curtailment and improve the use of the existing power transmission network. In [START_REF] Henry | Gym-ANM: Reinforcement learning environments for active network manage-ment tasks in electricity distribution systems[END_REF], RES curtailment was minimized by means of an optimal control strategy, which was also used to learn a computationally efficient control law based on reinforcement learning. Nevertheless, in this work storage and charging of EVs was assumed to be situated at a single point in the power network, whereas in reality both the EVs and the charging stations are distributed in time (in case of EVs) and space, and connected to different power grid nodes. In [START_REF] Zhou | Analyzing power and dynamic traffic flows in coupled power and transportation networks[END_REF], coupled traffic and power grid dynamics were considered, but the traffic flows were only described on graph level.

In this paper we propose for the first time in the literature an integrated model, combining the power grid model with a macroscopic electromobility model distributed in space, including charging stations under V2G operation (in Section 2). The proposed multi-class electromobility model extends the simplified CTEC model previously proposed in [START_REF] Čičić | Coupled macroscopic modelling of electric vehicle traffic and energy flows for electromobility control[END_REF] to capture a richer electromobility dynamics by allowing for several different co-located SoC levels in the vehicular flow. The electromobility layer is interfaced with the power layer through the charging stations, which act as predictable time-varying energy storage. This novel holistic model is used as a basis for optimal control design seeking to minimize RES curtailment, while accounting for the structural limitations of the grid and other SoC constraints necessary for the best possible operation of the EVs (in Section 3). Finally, the control framework is tested in simulations (in Section 4), based on which we are able to draw conclusions (in Section 5).

COMBINED ELECTROMOBILITY AND POWER GRID MODEL

In this section we present the combined electromobility and power grid model. We first introduce the new multi-class aggregated electromobility model including the charging stations, then we present the grid model, and finally connect the two parts into a combined model. The scenario studied in this work is similar to the one used in [START_REF] Henry | Gym-ANM: Reinforcement learning environments for active network manage-ment tasks in electricity distribution systems[END_REF], abstracting the situation where people commute between home and work using EVs, potentially stopping on the way at a public charging station.

As shown in Fig. 1, we consider a power grid with three nodes (Home, Public charging station, and Work, in further text denoted by h, p, and w, respectively) connected both by power lines and road links. The two considered road links connect nodes h and w, with a pair of off-and onramps at the middle, exiting towards node p. We assume that there is intermittent RES power generation and EV charging ports at all three nodes, and that there is some time-varying load at nodes h and w. While the combined model is adapted to this specific setup it can accommodate other types of road and power networks as well.

Electromobility model

In this section we introduce a multi-class extension of the simplified CTEC model [START_REF] Čičić | Coupled macroscopic modelling of electric vehicle traffic and energy flows for electromobility control[END_REF]

∂η(ε, t) ∂t + ∂(c(ε, t)η(ε, t)) ∂ε = 0, ( 3
)
where x is the position along the road, t the time, ρ(x, t) the traffic density, v(x, t) the traffic speed, ε(x, t) the SoC, D(v) the battery discharge as a function of EV speed, η(ε, t) the distribution of EVs at a charging station according to their SoC ε, and c(ε, t) the charging rate of EVs, potentially different for different ε. We extend the discretized version of the model by splitting the traffic flows into different classes which have the same behaviour while driving on the road, but may have independent SoC dynamics or behaviour at on-and off-ramps. These classes can be used to distinguish e.g. between combustion engine vehicles and EVs, or between vehicles with different routes and/or SoC, in order to avoid problems with the situation when two traffic flows with very different SoC merge.

Aggregate traffic density equations. Each road l is split into N l x cells of length L x satisfying L x ≥ v ff T , where v ff is the free flow speed of the traffic, and T is the sampling period for electromobility model which need to be suitable selected to fulfil the numerical stability requirements. The aggregate macroscopic state of EV traffic on each road link l is given by the traffic density ρ k i and SoC ε k i in each cell i at each discrete time instant k. Since here we consider multiple classes of vehicles, the aggregate macroscopic SoC ε k i is not used, but instead we will track the SoC of each class. The update equations for the aggregate traffic density are:

ρ k+1 i = ρ k i + T L x q k i--q k i+ , (4)
q k i+ = min v ff ρ k i , ω(ρ jam -ρ k i+1 ) -r k on,i+1 1 -β k i , (5)
q k i-=q k i-1+ (1 -β k i-1 ) + r k on,i , (6)
Here, q k i-denotes the traffic flow entering cell i at its upstream end, and q k i+ the traffic flow exiting cell i at its downstream end. The on-ramp flow into cell i is denoted r k on,i and assumed to enter it at its upstream end, and the off-ramp flow from cell i is assumed to leave it at its downstream end, and denoted r k off,i = β k i q k i+ , where β k i is the splitting ratio of mainstream traffic flow towards the off-ramp. The traffic is assumed to follow a triangular fundamental diagram with critical density ρ cr and jam density ρ jam , yielding congestion wave speed ω = v ff ρcr ρjam-ρcr , and we denote ρ k i = min{ρ k i , ρ cr } and ρ k i = max{ρ k i , ρ cr }. Note that we omit stating to which road link l we are referring for better readability.

Multi-class traffic density equations. The aggregate traffic flow may be split into some number of traffic classes ξ ∈ Ξ, with individual traffic densities denoted ξ ρ k i . Since here we assume that all vehicle classes have the same behavior in traffic, the evolution of their traffic densities is given by

ξ ρ k+1 i = ξ ρ k i + T L x ξ q k i--ξ q k i+ , (7)
ξ q k i+ = ξ ρ k i ρ k i q k i+ , (8)
ξ

q k i-= ξ q k i-1+ (1 -ξ β k i-1 ) + ξ r k on,i , (9)
where the ramp flows defined through ξ β k i-1 and ξ r k on,i are defined independently for each class ξ, with

ξ r k off,i = ξ β k i-1 ξ q k i+ .
The per-class quantities relate to the aggregate quantities according to

ρ k i = ξ∈Ξ ξ ρ k i , r k on,i = ξ∈Ξ ξ r k on,i , β k i = ξ∈Ξ ξ ρ k i ρ k i ξ β k i . ( 10
)
Once the vehicles reach the end of the road link from h to w or from w to h, all of them exit the road and enter the respective charging station, β k Nx = 1 for both road links. The splitting ratios towards the public charging stations ξ β k i , where i = i in ph for the link h to w, and i = i in pw for w to h, depend on the SoC of the approaching vehicles in the same cell i,

ξ β k i = 1 -1 + e - ξ ε k i -ε low γ -1 , (11)
as a monotonically decreasing sigmoid function centered at some ε low , with slope calibrated by some parameter γ. This function models the assumption that EVs are more likely to stop to charge if their SoC is low.

Multi-class SoC equations. The macroscopic SoC of each class is denoted ξ ε k i , and its evolution is defined by

ξ ρ k+1 i ξ ε k+1 i = ξ ρ k i ξ ε k i + ξ d k i T + T L x ξ ϕ k i--ξ ϕ k i+ , (12)
ξ ϕ k i-=( ξ q k i--ξ r k on,i ) ξ ε k i-1 + ξ d k i-1 T + ξ r k on,i ξ ε k on,i , (13)
ξ ϕ k i+ = ξ q k i+ ξ ε k i + ξ d k i T . ( 14
)
where the battery discharge rate ξ d k i depends on the traffic speed in cell i,

ξ d k i = D ξ v k i , v k i = q k i+ ρ k i , (15)
and is defined by potentially different discharge functions for each class D ξ .

Charging stations equations. Finally, the electromobility model is completed by the charging station dynamics.

For each charging station ζ, we split the SoC space into

N ε bins of length L ε , N ε L ε = 1. Each bin j = 1, . . . , N ε corresponds to a range of SoC [(j -1)L ε , jL ε ],
η k+1 j = ζ η k j + T L ε ζ c k j-1 ζ η k j-1 -ζ c k j ζ η k j -ζ c k j+1 η k j+1 ... ...+T ζ µ k in,j -ζ µ k out,j , (16)
where

ζ c k j = max{0, ζ c k j }, and ζ c k j = min{0, ζ c k j }, and ζ µ k in,j
and ζ µ k out,j represent the flows of vehicles entering and exiting the charging station, respectively.

The flows between the road and the charging station ζ are

r k on,i out ζ = Nε j=1 ζ µ k out,j , r k on,i out ζ ε k on,i out ζ = Nε j=1 ζ µ k out,j (j-1)L ε ,(17) r k off,i in ζ = Nε j=1 ζ µ k in,j , r k off,i in ζ ε k i in ζ = Nε j=1 ζ µ k in,j (j-1)L ε , ( 18
)
where i in ζ and i out ζ are the cells where the off-and onramp connecting the road to the charging station ζ are, respectively, ensuring both the vehicles and the energy are conserved. In this work, we assume that multiple classes of vehicles can enter a single charging station ζ,

ζ µ k in,j = ξ∈Ξ ζ,ξ µ k in,j , (19)
and that all vehicles exiting the charging station are of the same class ξ out ζ , therefore ξ r k on,i out

µ k in,j =                      j- ξ ε k i in ζ L ε   ξ r k off,i in ζ , j-1≤ ξ ε k i in ζ L ε <j,ξ∈Ξ in ζ ,   ξ ε k i in ζ L ε -j+2   ξ r k off,i in ζ , j-2≤ ξ ε k i in ζ L ε <j-1,ξ∈Ξ in ζ , 0, otherwise, (20)
ζ µ k out,j ∈ 0, 1 T - ζ c k j L ε ζ η k j + ζ c k j-1 L ε ζ η k j-1 - ζ c k j+1 L ε ζ η k j+1 , (21)
where Ξ in ζ denotes the set of vehicle classes that enter charging station ζ, and the exiting flow ζ µ k out,j depends on the particular behavioural logic of each charging station. We assume that each traffic class ξ enters at most one charging station ζ from each cell i, but it may enter multiple charging stations in case off-ramps leading to them are in different cells i in ζ1 ̸ = i in ζ2 . Charging stations control variables. Charging stations are of particular importance because they serve as the interface between the electromobility layer and the power grid, and can be used as actuators to improve the situation in the wider power system. At each time step k, the power system operator gets from each charging station ζ the range of power that they can consume or generate, denoted by U Here we adopt the convention that if u k ζ > 0, the charging station is a net power consumer from the perspective of the grid, and is using the grid power to charge the EVs. Otherwise, if u k ζ < 0, the charging station is a net power provider to the grid, and is using the energy stored in the batteries of some EVs to provide V2G services. The actual charging station powers, denoted U k ζ , represent the control input to the electromobility system, and are given by

U k ζ = u k ζ U k ζ , u k ζ ≥ 0, u k ζ U k ζ , u k ζ < 0, (22)
therefore we have

U k ζ ≤ U k ζ ≤ U k ζ .
This charging station power relates to its EV charging rates as

U k ζ = Nε j=1 ζ c k j ζ η k j B ( 23
)
where B is the EV battery average capacity. In order to best utilize the available power, we employ a hierarchical charging scheme, where vehicles with lower SoC have a higher priority for charging. We split the vehicles into three groups according to their SoC: low SoC EVs, with 0 ≤ ε < ε low , medium SoC EVs, with ε low ≤ ε ≤ ε high , and high SoC EVs, with ε high < ε ≤ 1. Additionally, in case there is not enough power to charge the low SoC EVs, as well as if the grid operators requests power to be provided to the grid, we utilize the energy stored in batteries of EVs with higher SoC, discharging their batteries to satisfy the power demands.

Given the boundaries between these groups of EVs ε low and ε high , we define j low as the highest j for which ε j < ε low and j high as the lowest for which ε j > ε high , j ∈ {1, . . . , N ε }. The number of EVs currently in each group is defined by

ζ η k low = j low j=1 ζ η k j , ζ η k med = j high -1 j=j low +1 ζ η k j , ζ η k high = Nε j=j high ζ η k j . ( 24
)
For each charging station we define the maximum charging rate C max

c k j =    ζ c k low , j = 1, . . . , j low , ζ c k med , j = j low + 1, . . . , j high -1, ζ c k high , j = j high , . . . , N ε . ( 25
)
The exception is charging at home, ζ = h, where we assume that all vehicles charge at the same rate,

h c k j = min P k h h η k low + h η k med + h η k high , C max h . ( 26
)
Given the current provided or demanded power U k ζ , these charging rates are given by

ζ c k low =min ζ P k low ζ η k low , C max ζ , (27)
ζ c k med =      C min ζ , ζ η k low C max ζ > ζ P k low , min ζ P k med ζ η k med , C max ζ , ζ η k low C max ζ ≤ ζ P k low , (28)
ζ c k high =      C min ζ , ζ η k med C max ζ > ζ P k med , min ζ P k high ζ η k med , C max ζ , ζ η k med C max ζ ≤ ζ P k med , (29)
where we define maximum available charging power to each group of vehicles, ζ P k low , ζ P k med , and ζ P k high , as

ζ P k low = U k ζ -ζ η k med + ζ η k high C min ζ , ( 30
) ζ P k med = U k ζ -ζ η k med C max ζ -ζ η k high C min ζ , ( 31
) ζ P k high = U k ζ -ζ η k high C max ζ . (32)
Finally, the presented charging scheme results in a range of power that charging station ζ can consume or produce

[U k ζ , U k ζ ] given by U k ζ = ζ η k low + ζ η k med + ζ η k high C max ζ B, (33)
U k ζ = max 0, ζ η k low C max ζ + ζ η k med + ζ η k high C min ζ ζ η k j B .( 34
)
Finally, we define the flows of EVs exiting each charging stations according to their scheduled departure demand

function ζ μk out , ζ µ k out,j =      min ζ μk out , 1 T - ζ c k j K ε ζ η k j , j = j k ζ , 0, j ̸ = j k ζ . (35)
For ζ = h and ζ = w, ζ μk out is given externally, to represent EVs commuting from home around t = 6 h, and from work to home around t = 16 h. For ζ = p, we assume that the EVs stay at the charging station for 1 h, so the departure demand will depend on the flow entering the charging station in the past. At each sampling time, we randomly select the SoC of the departing EVs for each charging station j k ζ , with probabilities of each SoC level proportional to the number of vehicles currently at charging station ζ.

Power Grid Model

This section describes the dynamics of the alternating current (AC) power grid that is employed to model the interaction between a conventional grid and the electromobility layer. To avoid confusion regarding notation conventions, and reduce abuse of notation, we indicate the parameters and variables of the power grid model by tilde.

Nowadays, the vast majority of AC transmission and distribution networks dispatch electricity using a three-phase system. Three parallel circuits are employed, each associated with its own phase, shifted 120 o from one another. Balanced three-phase system can be represented in phasor form as an equivalent single-phase representation, where only one of the phases is taken into account. In this work, power transmission lines are assumed to constitute such three-phase balanced networks, and we adopt its equivalent single-phase representation in the following derivations. Also, for mathematical convenience, the power system will be analyzed using the per-unit (p.u.) notation, in which all electrical quantities are normalized with respect to a set of base quantities. This normalization approach is widely used as it removes the need to include nominal voltage levels in derivations.

The distribution grid is modelled as a set of nodes Ñ linked by a set of directed edges Ẽ that represent the transmission lines. Each edge ẽij ∈ Ẽ links buses i and j and can contain transmission lines, power transformers and/or phase shifters. Here the conventional π transmission line model is employed where each branch can be defined by five parameters: a series resistance rij , a series reactance xij , a charging susceptance bij , a tap ratio magnitude τij and phase shift θij . These magnitudes allows to define the link series admittance ỹij = 1/(r ij + √ -1 xij ) and shunt admittance ỹsh ij = √ -1 bij 2 as well as the complex tap ratio of transformers, tij = τij e √ -1 θij . See [START_REF] Zimmerman | Matpower: Steady-state operations, planning, and analysis tools for power systems research and education[END_REF] for a more detailed explanation of the branch convention employed.

Power Network equations. For a network with Ñb = | Ñ | buses, and after applying Kirchhoff's current law at each bus i ∈ Ñ , the constant impedance elements of the model are incorporated into a complex Ñb × Ñb bus admittance matrix Ỹ that relates the complex nodal current injections, Ĩ to the complex node voltages, Ṽ , as Ĩ = Ỹ Ṽ .

The impedance matrix admits a systematic definition as

Ỹij =                      - ỹij t * ij , i ̸ = j, ẽij ∈ Ẽ, - ỹij tij , i ̸ = j, ẽji ∈ Ẽ, ẽik ∈ Ẽ ỹik + ỹsh ik tik t * ik + ẽki ∈ Ẽ ỹki + ỹsh ki , i = j, 0,
otherwise. Power flow solution. The standard power flow or load flow problem involves solving for the set of voltages and flows in a network corresponding to a specified pattern of load and generation. The definition of grid admittance matrix Ỹ allows to define the nodal equation in an equivalent form in terms of power injections and voltage level at buses as

Sbus i = Ṽi Ĩ * i = Ṽi ( Ỹi Ṽ ) * = Ṽi Ỹ * i Ṽ * , ∀i ∈ Ñ , (36)
where

Sbus i = P bus i + √ -1 Qbus i
is the complex power phasor accounting for active ( P bus i ) and reactive ( Qbus i ) powers injected at bus i, and Ỹi denotes the i-th row of the admittance matrix Ỹ . Equation ( 36) is a set of Ñb complex-valued equations, that can be posed as 2 Ñb quadratic real valued equations with 2 Ñb real unknowns. By convention, a single generator bus is typically chosen as a reference bus to serve the roles of both a voltage angle reference and a real power slack. The voltage angle at the reference bus has a known value, but the real power generation at the slack bus is taken as unknown to avoid overspecifying the problem. Unknowns for the problem will thus be active and reactive powers at the slack bus together with complex voltage level at the remaining Ñb -1 buses.

These equations are solved at every iterations step of the problem providing equilibrium powers and voltages at the bus level. The solver of choice employed in this work is based on a standard Newton's method [START_REF] Tinney | Power flow solution by newton's method[END_REF] using a polar form and a full Jacobian updated at each iteration. Each Newton step involves computing the mismatch, forming the Jacobian based on the sensitivities of these mismatches to changes in and solving for an updated value of the unknowns by factorizing this Jacobian.

Combined model

We denote the encapsulated full state of the mobility model at time k as E k ,

E k = ξ ρ k l,i |l ∈ {hw,wh}, i = {1, . . . , N x }, ξ ∈ {h/w, p}, . . . ξ ε k l,i |l ∈ {hw, wh}, i = 1, . . . , N x , ξ ∈ {h/w, p}, . . . ζ η k j |j ∈ {1, . . . , N ε }, ζ ∈ {h, w, phw, pwh} , (37)
capturing the traffic density ξ ρ k l,i and SoC ξ ε k l,i for all vehicle classes ξ and all cells i on both road links l, and the numbers of charging vehicles ζ η k j at each SoC level j at each charging station ζ. The vehicles are split into classes based on whether they enter the road from h or w (h/w), or from p. Physical charging station p is modelled by two virtual charging stations pwm and phw, based on what link it is connected to.

The electromobility state is updated according to ( 7), ( 12), and ( 16), as defined in Section 2.1, and we write jointly

E k+1 = M(E k , U k ), (38)
where U k is the collection of U k ζ to all physical charging stations, ζ ∈ {h, w, p}, yielding charging rates (25) that directly determine the evolution of ζ η k j through ( 16). Power grid equations at discrete time step k is given by a static mapping, and can be expressed compactly as

F k = F(U k , Λk , Ṽ k ), (39)
P k c = C( F k , Ṽ k ), (40)
P k r = R(E k ), (41)
where F k are the power flows at time k as a function of the EV power injected/retrieved U k , loads Λk , and RES distributed generation Ṽ k . RES curtailment is is given by P k c , and mobility power reserves by P k r,ζ . Except for power flows F k , all other variables in (39)-( 41) are collection of variables related to different ports ζ ∈ {h, w, p}.

OPTIMAL CONTROL DESIGN

Power system optimization

Active network management (ANM) refers in the power system literature to the decision-making procedures to design control schemes that modulate loads, generators, and the distributed energy resources connected to the grid, in order to keep power grid operation within nominal conditions in terms of bus voltages, stable frequency and power lines capacities.

The modulation resulting from ANM often causes a necessary reduction in the output of generators from what they could otherwise have produced given available resources, often referred to as the process of curtailment. Such generation curtailment, along with storage and transmission losses, constitute a major source of energy management inefficiencies that this work aims to improve. In addition, the ANM scheme must ensure a safe and reliable operation of the distribution network, expressed as a set of technical and operational constraints that must be satisfied.

The proposed optimization problem focuses on the reduction of curtailments by providing near future scheduling of the optimal power flow policy to/from EVs plugged at the different charging poles of the system (h, w, and p). This separation is designed to collect the spatiotemporal operation of a large majority of EV users commuting daily from home to work, and occasionally making use of public charging stations facilities.

The optimization problem is formulated based on a prediction model and operated in a receding horizon manner using the latest system information. In this approach, load injections and generation profiles of RES are noncontrollable variables. Only the active power (and reactive, within operational limits) of EV demand is assumed to be manageable. The objective function, constraints, and optimization algorithm will be presented next.

Optimal control formulation

We can now formulate ANM taking electromobility into account as a receding horizon optimal control problem minimize

ũk ζ J( k, ũk ζ )
subject to Mobility and grid dynamics (38)-( 41), Grid operational constraints, 43), and limitations on the amount of reactive to active power ratios the distributed generation can manage.

| Sij | ≤ Sij , ẽij ∈ Ẽ, (42)
Ṽ i ≤ Ṽi ≤ Ṽ i , i ∈ Ñ -Ñslack , (
Note that in order to ease the numeric burden, we solve this problem at a time scale different than that of the electromobility model. We take the optimization sampling period T as a multiple of electromobility sampling periods T , and denote the optimization time step k. The resulting control ũk ζ is then applied to the electromobility model over a number of time steps,

u k ζ = ũ⌊k T T ⌋+1 ζ , (45)
where ⌊ • ⌋ denotes rounding down.

The first part of the cost function accounts for the overall energy that is lost in the system as a result of necessary curtailment imposed by the ANM algorithm to satisfy grid constraints for a predicted sequence of control actions ũζ (t k+h |t k ), ζ ∈ {h, p, w}, h =, 1..., N u . The control actions correspond to the normalized active power that the controller schedules to be retrieved or injected (V2G) at port ζ, and N u is the control horizon considered. In this scheme, as is usually employed in the classic formulation of predictive control, for N u > 1 only the first control action of the sequence computed is applied. The control action terms ũζ (t k+h |t k ) ∈ [-1, 1] are considered in a normalized fashion, allowing to accommodate the constraints on the time-varying mobility reserve powers in a more convenient form. The second term is designed to maximize the mobility reserve storage available within the prediction horizon, which is taken for simplicity as a proxy of the overall SoC of the EV fleet.

SIMULATION RESULTS

Finally, we assess the behaviour of the proposed ANM EV charging scheme by comparing two simulation scenarios.

The first one corresponds to uncoordinated charging where a fixed percentage of EVs at each charging station ζ is charging at constant maximum rate C max ζ . We set these percentages to 70% for ζ = h, and 40% for ζ = w and ζ = p. The second scenario corresponds to a scheduled charging operation as proposed in this work using a receding horizon control scheme that balances curtailment losses in relation to average EV reserves. Note that in both cases vehicles decide whether or not they enter the public charging station depending on their SoC, according to (11). Algorithm 1 outlines the main steps to perform the proposed simulations and control design approach. Dynamic equations have been discretized to be approximated by a conventional fixed-step iteration algorithm.

Both scenarios have been simulated on the test layout described in Fig. 1 where for the sake of simplicity, the power network is reduced to three representative ports that lump loads, distributed generation and EV operation for aggregated domestic, work and public EV charging stations. Simulations were performed assuming a grid configuration with balanced capacities on the three ports in the order of 30 MW. Power lines capacities were assigned Algorithm 1 Sketched simulation algorithm 1: procedure GridMobSimulation Due to the daily cyclic nature of power loads, RES production and mobility patterns, simulations for a period of 24h were considered representative of the system. For the computation of the coordinated charging schedule, we used a prediction horizon of H T = 5 h, and a control horizon of N u T = 30 min. To speed up computations, simulations were performed with different sample times, with T = 1 min for the electromobility model, 5 min for the power grid model and T = 30 min as sample time for the receding horizon control.

As a primary objective of the optimization problem proposed, the performance index in (44) has been balanced to maximize the use of RES (reduce curtailment losses), while preserving a reasonable average state of charge of the EVs. Figure 2 plots the compared curtailment induced by the uncoordinated and coordinated schemes, the shaded area corresponding to curtailment energy saving by the later approach. More precisely, in this particular simulation, the uncoordinated charging scenario induces 249.7 MWh losses in curtailment, while the coordinated approach provides 124.1 MWh, a significant reduction of 50.3%.

Figure 3 shows the evolution of the power used by EVs.

For the uncoordinated charging scenario, the power used for charging is directly proportional to the number of EVs currently present at the charging station. Coordinated charging allows for a more nuanced use of power depending on the power grid load and RES production at different times of the day.

In order to illustrate the state of the electromobility model, we show the state of charging stations over time in Figure 4. As expected, all EVs start at home, then travel to work, and finally return home, with some of them stopping at the public charging station along the way. In this case EVs end the day with somewhat lower SoC, which will be rectified through overnight charging at home.

CONCLUSIONS

In this work, we present a control framework combining electromobility and the power grid. Using an extended multi-class electromobility model, we are able to capture the SoC dynamics of EVs on their daily commute. The electromobility layer is interfaced with the power grid through charging stations, providing us with a way to store energy and control the power flows. We use this framework to design a receding horizon control law for for active network management, in order to reduce RES production curtailment, which is shown in simulation to significantly improve the situation.

The models and results presented herein could serve as a basis for a number of future research directions. Though the network studied here was adopted for simplicity, the presented control approach should be extended to a more general case and more complex networks, both in terms of roads and power lines. In this work, we did not utilize the ratio of EVs that enter the public charging station as a control input, but this behavior could be influenced by e.g. introducing dynamic pricing of charging, which opens an entirely new area of control problems.

Fig. 1 .

 1 Fig. 1. Layout of the combined electromobility and power grid model used for this study.

=0

  for ξ̸ =ξ out ζ . The flows entering and exiting the charging stations are further determined by ζ,ξ

  k ζ and U k ζ , respectively. Based on these limits, the operator can choose the normalized current power flow to (or from) each charging station ζ, denoted u k ζ ∈ [-1, 1].

ζ≥ 0

 0 and maximum discharging rate C min ζ ≤ 0 per EV, yielding C min ζ ≤ ζ c k j ≤ C max ζ , respecting |C max ζ | ≤ C and |C min ζ | ≤ C. We assign a different charging rate to each group of EVs, denoted ζ c k low , ζ c k med , and ζ c k high , and given by ζ

  Optimize for horizon H ( Pr ) 8: Ẽk +1 ← Update Mobility Model(ũ k) 9: Solve Instant Power Flow Balance (ũ k) 10: k ← k + 1

Fig. 2 .

 2 Fig. 2. Curtailment cut. Uncoordinated vs coordinated charging scenario

  and we describe the state of the charging station through the number of vehicles in each bin ζ η k j . Different classes of vehicles at a single physical charging station can be represented by using multiple virtual charging stations ζ to represent charging of each vehicle class ξ. We allow the rate of charging ζ c k j for each charging station ζ and for each SoC level j, to vary in time within some range ζ c k j ≤ C, with L ε ≥ CT required for numeric stability. The charging station state update is then given by

	ζ

  Pc,ζ (t k+h |t k)-w r Pr,ζ (t k+h |t k) (44) where Pc,ζ (t k+h |t k) represent the curtailed power at time t k+h inferred from the information available at time t k, at port ζ. Similarly Pr,ζ (t k+h |t k) accounts for the mobility reserve power available at time t k+h . Constants w c and w r weigh the different optimization terms. The optimization problem is subject to constraints on the maximum apparent power | Sij | (42), power line voltage Ṽi deviation at buses (

		43)
		Reactive to active power ratio limitations.
	The optimization problem is articulated as constrained
	minimization of the cost function
	H
	J( k,ũ k ζ )=	w c
	h=1 ζ∈{h,p,w}

ACKNOWLEDGEMENTS

This work has received support from the Scale-FreeBack project, funding by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement N 694209), and by projects PID2020-115561RB-C32 granted by the Spanish Ministry of Science and Innovation.