Sylvain Arlot
email: sylvain.arlot@u-psud.fr

Stefano Marmi
email: stefano.marmi@sns.it

Duccio Papini
email: duccio.papini@uniud.it

Coupling the Yoccoz-Birkeland population model with price dynamics: chaotic livestock commodities market cycles

come   L'archive ouverte pluridisciplinaire

Introduction

Twenty years ago a new model of population dynamics which exhibits endogenous chaotic behaviour has been proposed by J.-C. Yoccoz and H. Birkeland [START_REF] Yoccoz | Informal commentaries on the numerical investigation of the toy model[END_REF], (we refer to [START_REF] Sylvain Arlot | Dynamique de populations de petits mammifères, saisonnalité et attracteur de Hénon -comment une question d'écologie a pu intéresser Jean-Christophe[END_REF] for a historical perspective on this work). The model was prompted from the evidence of aperiodic large oscillations (2-3 orders of magnitude on a 3-5 years time span) in the time evolution of the population of a species of rodents, Microtus Epiroticus (sibling vole) on Svalbard Islands. This species has a high fertility rate which has a strong dependence on seasonal factors (due to harsh weather conditions in winter) and on the population density. Indeed, few good reproduction spots are available and their quality decreases as the population increases. Despite the absence of any significant predator and the relative abundance of food, one observes high oscillations of population.

The Yoccoz-Birkeland model was studied through a mathematical analysis and some simulation experiments in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF][START_REF] Nieto | Long-Term and Short-Term Dynamics of Microtus epiroticus: A Yoccoz-Birkeland Model[END_REF]. In short, it has been showed that such a deterministic model can produce complex dynamics with a high sensitivity to initial conditions, only by the combination of density-dependent fertility, the lag due to the maturation age and a periodic seasonality. A detailed account of these results is provided in Section 2.

Date: June 11, 2018 Preliminary version. This work was started during a visit of SA to the Scuola Normale Superiore in Pisa and to the Complex Systems Center of the University of Siena, the hospitality of which is gratefully acknowledged. The authors acknowledge the support of the Centro di Ricerca Matematica Ennio de Giorgi and of UniCredit Bank R&D group for financial support through the "Dynamics and Information Theory Institute" at the Scuola Normale Superiore. Most of this work was done while SA was a researcher at CNRS (Département d'Informatique / École Normale Supérieure, UMR 8548 CNRS/ENS/INRIA), 45, rue d'Ulm, F-75230 PARIS Cedex 05, France). DP was supported by the INdAM-GNAMPA project "Proprietà qualitative di alcuni problemi ai limiti". We are also grateful to Paolo Nistri for contributing to the initial formulation of the model and to Nicola Zanda for useful discussions on livestocks farming (in particular Cinta Senese and pig).

In this paper, we introduce a new model coupling the population dynamics in the Yoccoz-Birkeland model with an equation modelling the price dynamics of a livestock commodity market inspired by [START_REF] Bélair | Consumer memory and price fluctuations in commodity markets: an integrodifferential model[END_REF]. A cattle population is split at the birth into reproducing females and cattle for butchery. The relative amount of the two is determined by the spot price of the meat whereas the logarithmic derivative of the price is determined by the unbalance between the demand and supply of the meat. On the population side, seasons (or artificial synchronization of births) and maturation lags are also taken explicitely into account.

The importance of the livestock commodities market in economics is related to the fact that it provides some of the oldest and best documented examples of business cycles. Approximately periodic fluctuations of supply and prices were first observed at the beginning of last century in hog markets [START_REF] Hanau | Die Prognose der Schweinepreise[END_REF], [START_REF] Ezekiel | The cobweb theorem[END_REF] and since then they have been the object of many studies. Beef cattle stocks "are among the most periodic time series in economics" [START_REF] Rosen | Cattle cycles[END_REF], a fact also related to the double role played by cattle as capital as well as consumption goods.

In many respects, the continuing presence of any price cycle is disturbing: if a predictable price cycle exists, then producers responding in a countercyclical fashion could earn larger than "normal" profits over time [START_REF] Dermot | Hog cycles and countercyclical production response[END_REF]. Such profits could occur even with lags in the production process (substantial gestation and maturation times stretch it over lengthy intervals of time) because predictable price movements would still influence production decisions. Eventually, countercyclical production response would smooth out price fluctuations at the market level, causing the cycle to disappear.

An alternative explanation for the existence of a business cycle is that the cycle itself is not perfectly predictable: the law of motion may be a deterministic nonlinear relationship that generates unpredictable patterns [START_REF] Grandmont | Nonlinear economics dynamics: Introduction[END_REF]). Cobweb models show that complicated price dynamics may indeed occur due to nonlinearities [START_REF] Chiarella | The cobweb model: its instability and the onset of chaos[END_REF][START_REF] Hommes | Dynamics of the cobweb model with adaptive expectations and nonlinear suppy and demand[END_REF] and that simple expectation rules in a nonlinear environment may lead to chaotic price fluctuations [START_REF] Hommes | Behavioral Rationality and Heterogeneous Expectations in Complex Economic Systems[END_REF].

Our model shows how, under quite natural assumptions, simply connecting the percentage of reproducing females with the price equation gives rise to a chaotic time evolution of price. This is characterized by a series of "cycles" of booms and busts (i.e. rapid increase or decrease). Another important feature of the model is that it takes into account some specific characteristics of the production of livestock commodities often neglected in the literature: for example the existence of timelags between the producer decision on the reproduction strategy and the butchery of the calves, the synchronization of births and their seasonality. We point out that some reproductive constraints included in Yoccoz-Birkeland model are realized also in the production of livestock commodities, besides synchronization of births and seasonality. For example, the development of breeding facilities requires time and limits the reproductive capacity.

Here follows the plan of our paper. In Section 2 we review the model proposed by Yoccoz and Birkeland and the relative results that were shown in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF][START_REF] Nieto | Long-Term and Short-Term Dynamics of Microtus epiroticus: A Yoccoz-Birkeland Model[END_REF].

Our model is derived in Section 3 by coupling in a suitable way the Yoccoz-Birkeland model with a differential equation for the logarithmic derivative of the price which is inspired by [START_REF] Bélair | Consumer memory and price fluctuations in commodity markets: an integrodifferential model[END_REF].

In Section 4 our model is rigorously analyzed. Global existence and uniqueness for initial value problems are obtained along with some useful estimates on the solutions. In particular we show that global boundedness and persistence follows under suitable assumptions involving some relevant biological parameters. Finally we prove the existence of a global attractor containing at least a non-trivial periodic solution.

In Section 5 numerical experiments show this new model can produce complex dynamics in the population size and the price. The main setting of the parameters is chosen having in mind the attractor detected in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] for Yoccoz-Birkeland model, on the one hand, and the hog market, on the other. The attractor we found has both sensitive dependence on initial conditions and noninteger dimension. The relevance of the presence of the market dynamics, based on the unbalance of demand and supply, is outlined in a second numerical experiment in which the population dynamics is suitably decoupled from the price evolution and gives rise to a behavior which looks like by no means chaotic, but asymptotically periodic, in fact. Bifurcation diagrams shows that complex dynamics persists for realistic values of the maximal fertility and also for weaker levels of dependence of the fertility rate on the total population. Technical details of the numerics are given in the Appendix.

The Yoccoz-Birkeland model

This section recalls the model proposed in [START_REF] Yoccoz | Informal commentaries on the numerical investigation of the toy model[END_REF] and gives a brief summary of results obtained in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF][START_REF] Nieto | Long-Term and Short-Term Dynamics of Microtus epiroticus: A Yoccoz-Birkeland Model[END_REF].

The Yoccoz-Birkeland model aims at modelling the population of mature females of a single species with density-dependent reproduction rate and whose reproduction strategy may be influenced by seasons or other external factors.

More precisely, the model proposed by Yoccoz and Birkeland goes as follows:

• t is the time measured in years;

• A 0 is the age (in years) at which females reach sexual maturity;

• A 1 is the maximum age (in years) for females;

• N (t) is the number of sexually mature females at time t;

• m(N ) is the density-dependent female reproduction rate and measures the average number of female cubs that a single female can give birth to in a year in optimal weather conditions when the total number of mature females is N ; it is reasonably a decreasing function of N ; • m ρ (t) is the seasonal factor and gives the fraction of females actually reproducing at time t; typically it is 1-periodic; • S(a) is the fraction of newborn females still alive at age a (in years); then the number of females with age ranging in [a, a + da] is given by

N (t -a)m N (t -a) m ρ (t -a)S(a)da .
Therefore, N satisfies the following integral equation (2.1)

N (t) = A1 A0 N (t -a)m N (t -a) m ρ (t -a)S(a)da ,
which allows to uniquely determine N (t) for t ∈ [t 0 , t 0 + A 0 ] (and for every other

t > t 0 by recursion) if N (t) is known for t ∈ [t 0 -A 1 , t 0 ].
This model have been proposed to explain the behaviour of the population of Microtus Epiroticus in Svalbard Isles which has high average fertility, while its numbers show large fluctuations so that sometimes the species looks to be even close to extinction in spite of the absence of any significant predation [START_REF] Yoccoz | Demography of small mammals in cold regions: the importance of environmental variability[END_REF]. In the case of this species of small rodents the biological explanation of the observed behavior relies on the following remarks. First, the adverseness of the environment in which the population lives causes a shortage of good reproduction spots and, as a consequence, the larger is the number of sexually active females (above some threshold) the smaller becomes the average individual fertility because the spots become overcrowded. Secondly, the shortness of the reproduction season (summer) induces a decrease of the age at which females reach sexual maturity. In this way, the females that are born at the beginning of summer are able to give birth to their first cubs before the end of the same summer. These facts translate into a density-dependent reproduction rate that decreases as the population increases and into a sexual maturation age that is less than the average length of summer.

In [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] the following choices were made:

S(a) = 1 - a A 1 , for a ∈ [0, A 1 ] m ρ (t) = 0 if 0 ≤ t < ρ (mod 1) 1 if ρ ≤ t < 1 (mod 1) (2.2) m(N ) = m 0 if N ≤ 1 m 0 N -γ if N > 1 (2.3)
where ρ ∈ (0, 1) stands for the average length of winter, m 0 is the average yearly female fertility in optimal weather and environmental conditions and γ ≥ 1.

Equation (2.1) gives rise to a continuous semi-group as follows. For each

t 0 ∈ R/Z set Y t0 = N ∈ C [-A 1 , 0] : N (0) = A1 A0 N (-a)m N (-a) m ρ (t 0 -a)S(a)da
and consider the phase space

Y = {(t, N ) : t ∈ R/Z, N ∈ Y t }
which is a complete metric space with respect to the distance d((s, M ), (t, N )) = |s -t| R/Z + M -N ∞ . Then the semi-group (T s ) s≥0 generated by (2.1) is given by T s (t, N ) = (t + s(mod 1), N s t ) where:

N s t (-a) =      N (s -a) if 0 ≤ s ≤ a ≤ A 1 A1 A0 N (s -a -b)m N (s -a -b) m ρ (t + s -a -b)S(b)db otherwise.
The following result holds.

Theorem 2.1 ([Arl04]). Assume that N → N m(N ) is uniformly continuous on [0, +∞) and that (1) m 0 /2 ≤ m(N ) ≤ m 0 if N ≤ 1 and min{1/2, N -γ }m 0 ≤ m(N ) ≤ m 0 N -γ if N ≥ 1, γ ≥ 1; (2) 0 ≤ m ρ (t) ≤ 1 for all t and m ρ (t) = 1 on an interval of length 1 -ρ > 0, ρ > 0; (3) A 1 ≥ max{2A 0 , A 0 + 1} and c 0 m 0 > 2 where c 0 = A0+1 A0+ρ S(a)da. Moreover, let N max = m 0 A 1 2 1 - A 0 A 1 and L = m 0 3 - A 0 A 1
and consider the set

K = (t 0 , N ) ∈ Y : c 0 m 0 2 N 1-γ max ≤ N (s) ≤ N max ∀s ∈ [-A 1 , 0] and N (s 1 ) -N (s 2 ) ≤ L|s 1 -s 2 | ∀s 1 , s 2 ∈ [-A 1 , 0] .

Then

(1) K is compact and T s (t 0 , N ) ∈ K for all s ≥ 0 and all (t 0 , N ) ∈ K;

(2) for each (t 0 , N ) ∈ Y there exists s 0 ≥ 0 such that T s (t 0 , N ) ∈ K for all s ≥ s 0 ;

(3) the compact set Λ = n≥0 T n (K) is invariant: T s (Λ) = Λ for all s ≥ 0;

(4) for each neighbourhood U of Λ and each (t 0 , N ) ∈ Y there exists s 0 ≥ 0 such that T s (t 0 , N ) ∈ U for all s ≥ s 0 ; (5) Λ is a global attractor of (T s ) s≥0 , Y .

Simulations were also performed in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] on a discretization of (2.1) with a smoothed version of (2.2) and (2.3) for m ρ (t) and m(N ), and the choices m 0 = 50 and A 1 = 2. The parameter space was explored with respect to γ ∈ [2, 16], ρ ∈ [0, 0.5], A 0 ∈ [0, 0.4]. The computations showed the existence of periodic points, of possible Hopf bifurcations, the coexistence of different attractors and the presence of complex dynamics. In particular, a complex attractor is outlined in the cases ρ = 0.30, γ = 8.25 and A 0 = 0.15 or A 0 = 0.18; a detailed study of its dynamical features is done for A 0 = 0.15. In any case, the analysis and the simulations in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] show that the Yoccoz-Birkeland model recovers to some extent the general behavior and the biological characteristics of Microtus Epiroticus outlined above. Numerical solutions of (2.1) have large oscillations, with minima close to extinction, thanks to the interaction of the density-dependent fertility, the relatively quick sexual maturation of females and the average duration of winter.

Recently, the Yoccoz-Birkeland model has been the subject of the paper [START_REF] Nieto | Long-Term and Short-Term Dynamics of Microtus epiroticus: A Yoccoz-Birkeland Model[END_REF]. There the analogue of Theorem 2.1 and the existence of periodic points have been proved for the discrete version of the model and numerical simulations have been done with special emphasis to the case with A 0 = 0.18.

A model coupling market and population dynamics

This section presents the new model proposed in this paper for cattle population and price dynamics. The idea is to couple a population dynamics model similar to the Yoccoz-Birkeland with a market dynamics model. The (cattle) population is split into two parts: on the one hand, females for reproduction; on the other hand, cattle for butchery (all the males plus some of the females).

The mechanism is the following:

• At the birth of some babies, part of the newborn females are put in the reproduction line, and the remaining newborn females are put in the butchery line together with all newborn males. The fraction R of newborn females that will reproduce (chosen by the breeder) is only determined by the price of meat at birth time. • In the reproducing line, females have children between ages A 0 and A 1 .

Their fertility can be affected by seasons, or because births are synchronized by the breeder (through a function m ρ (t)). Reproducing females older than A 1 (hence, non fertile) are not taken into account anywhere in the model. • In the butchery line, cattle can be butchered between ages Ω 0 and Ω 1 . So, only the (alive) butchery population between ages Ω 0 and Ω 1 can count as a "supply" for the market. • The price evolution is a simple function of the supply (which comes from the butchery line population dynamics) and the demand (which depends only on the price). Note that contrary to the Yoccoz-Birkeland model, we assume no mortality before ages A 1 (resp. Ω 1 ).

Notation and parameters.

• t is the time measured in years.

• N r (t) is the total population of mature females that are in the reproducing line and can give birth to pups at time t.

• N b (t) is the total population of cattle that is suitable for butchery at time t (both males and non-reproducing females, old enough and in the butchery line). • R(P ) is the fraction of newborn females that are put in the reproducing line when the price of meat is P when they are born. • A 0 is the age from which females can have children (i.e., age of sexual maturity + length of the first gestation). • A 1 is the maximal age at which females can give birth to children (i.e., age of sexual unfertility + length of the last gestation).

• Ω 0 is the minimal age at which the cattle (male or female) can be butchered.

• Ω 1 is the maximal age at which the cattle (male or female) can be butchered.

Note that Ω 1 could possibly be enlarged (compared to its biological value) if the meat can be frozen after butchering. • m(N ) is the average annual female (resp. male) fertility of each mature female when the total population is N , i.e., the average number of female (resp. male) babies per year for a single mature female. Typically it is a decreasing function as in (2.3). We assume the sex ratio is 1/2, i.e., the average number of male babies is equal to the average number of female babies (hence, m(N ) is half of the average annual fertility). • m ρ (t) is the 1-periodic step function (with 1 0 m ρ (t)dt = 1) that accounts for a possible modulation of fertility during each year (births synchronization or seasonal effects).

• P (t) is the market price of meat at time t.

• D(P ) is the demand of the market (per time unit) when the price of meat is P (typically a decreasing function of P ). • S(t) is the supply to the market (per time unit) at time t (typically proportional to N b (t)). • λ is a "temperature" parameter of the meat market: higher values of λ correspond to bigger price variations in response to the same demand/supply imbalance. • F (D, S) is the function of demand and supply that rules the meat price dynamics.

Population dynamics.

The population dynamics model is strongly inspired from the Yoccoz-Birkeland model (2.1). In order to derive the equations satisfied by N r and N b , let us define the following additional notation:

• B f (t)
is the density of newborn female cattle at time t (i.e., B f (t)dt females are born between t and t + dt). • B m (t) is the density of newborn male cattle at time t. • B r (t) is the density of newborn (female) cattle that are put in the reproducing line at time t. • B b (t) is the density of newborn cattle that are put in the butchery line at time t. First, the male and female birth densities at time t are given by:

(3.1) B f (t) = B m (t) = m ρ (t)m N r (t) N r (t) .
The breeder decides at birth time t which fraction R(P (t)) of the newborn females is going into the reproducing line, which gives:

B r (t) = B f (t)R P (t) = m ρ (t)m N r (t) N r (t)R P (t) (3.2) and B b (t) = B m (t) + B f (t) 1 -R P (t) = m ρ (t)m N r (t) N r (t) 2 -R P (t) . (3.3)
We assume no mortality at all between birth and the end of reproduction time for females, so the number of mature females is given by (3.4)

N r (t) = A1 A0 B r (t -a)da = A1 A0 m ρ (t -a)m N r (t -a) N r (t -a)R P (t -a) da .
Similarly, we assume no mortality at all between birth and butchering time for males and females in the butchery line, so the size of the cattle population suitable for butchery is (without any butchering before age Ω 1 ) (3.5)

N b (t) = Ω1 Ω0 B b (t-a)da = Ω1 Ω0 m ρ (t-a)m N r (t-a) N r (t-a) 2-R P (t-a) da .

Market dynamics.

Inspired by [START_REF] Bélair | Consumer memory and price fluctuations in commodity markets: an integrodifferential model[END_REF] (see also [START_REF] Michael | Commodity price fluctuations: price dependent delays and nonlinearities as explanatory factors[END_REF]), we consider the following differential equation satisfied by the price as a function of the demand D(P ) and the supply S(t):

(3.6) P (t) P (t) = λF D P (t) , S(t) where F (D, S) = D -S D + S .
Other functions F could be considered such as F (D, S) = (D -S)/S. The parameter λ > 0 measures the "temperature" of the market, i.e., how fast can the price goes up or down. The function P → D(P ) is a decreasing function of the price P .

In order to define the supply function S, we assume that all the cattle in the butchery line is butchered exactly at age Ω 1 (and never before), while the market takes into account all the population N b (t) for determining the price in equation (3.6). This leads to choosing (3.7)

S(t) = N b (t) ∆Ω = 1 ∆Ω Ω1 Ω0 m ρ (t -a)m N r (t -a) N r (t -a) 2 -R P (t -a) da , for ∆Ω = Ω 1 -Ω 0 .
3.4. The population/market model. Our model then consists in coupling equations (3.4), (3.6) and (3.7):

N r (t) = A1 A0 m ρ (t -a)m N r (t -a) N r (t -a)R P (t -a) da (3.8a) P (t) P (t) = λF D P (t) , S(t) (3.8b) S(t) = 1 ∆Ω Ω1 Ω0 m ρ (t -a)m N r (t -a) N r (t -a) 2 -R P (t -a) da (3.8c) where F (D, S) = (D -S)/(D + S), and m : [0, +∞) → [0, +∞), m ρ : R → [0, +∞), R : [0, +∞) → [0, 1], D : [0, +∞) → [0, +∞), A 1 > A 0 > 0, Ω 1 > Ω 0 > 0, ∆Ω = Ω 1 -Ω 0
and λ > 0 have to be chosen. Possible choices for all these parameters are discussed in Section 5.1.

Analysis of the model

The section analyses mathematically the model defined by equations (3.8a)-(3.8c), under the following assumptions:

• m ρ : R → R is a non-negative, bounded, 1-periodic function such that 1 0 m ρ = 1 and we let m ρ (t) ≤ µ max and

(4.1) 0 < c 0 ≤ A1 A0 m ρ (t -a)da ≤ c 1 ∀t; • m : [0, +∞) → R is a continuous function that satisfies m 0 2 min{1, N -γ } ≤ m(N ) ≤ m 0 min{1, N -γ } ∀N > 0 with m 0 > 0 and γ ≥ 1; • R : [0, +∞) → R is a continuous function such that R 0 ≤ R(P ) ≤ R 1 for all P ≥ 0 and some constants R 1 , R 0 > 0 with R 1 ≤ 1; • D : [0, +∞) → R
is a strictly decreasing and locally Lipschitz continuous function such that D(+∞) = 0 and we set D 0 = D(0). We begin by setting up a phase space and a notion of solution suitable for our model. Let

T 0 = min{A 0 , Ω 0 }, T 1 = max{A 1 , Ω 1 } and X = L ∞ ([-T 1 , 0]; [0, +∞))× C 0 ([-T 1 , 0]; [0, +∞))
which is a complete metric space with respect to the distance induced by the norm

(N, P ) X := N ∞ + P ∞ := ess sup s∈[-T1,0] N (s) + sup s∈[-T1,0] P (s) .
In particular, when we consider N ∈ L ∞ ([-T 1 , 0]; [0, +∞)), we actually mean that ess inf N ≥ 0. Definition 4.1. Let (N 0 , P 0 ) ∈ X , and t 0 , T ∈ R with t 0 < T . A solution of (3.8a)-(3.8c) with initial data (N 0 , P 0 ) is a couple (N r , P ) :

[t 0 -T 1 , T ) → R 2 such that N r | [t0,T ) is continuous, P | [t0,T ) is differentiable, N r , P satisfy (3.8a)-(3.8c) for t ∈ [t 0 , T ), while N r (t 0 + a) = N 0 (a) and P (t 0 + a) = P 0 (a) for a ∈ [-T 1 , 0).
Our first result shows that a unique solution exists, is globally defined and satisfies some estimates: N r and S are globally bounded and the component N r turns out to be Lipschitz continuous on [t 0 , +∞). In particular all the obtained estimates are uniform with respect to the initial condition.

Proposition 4.2. Let (N 0 , P 0 ) ∈ X and t 0 ∈ R be given. Then there exists a unique solution pair (N r , P ) : [-T 1 + t 0 , +∞) → R 2 of (3.8a)-(3.8c) with initial data (N 0 , P 0 ). Moreover, N r , P are non-negative and

N r (t) ≤ N max ∀t ≥ t 0 N r (t) -N r (s) ≤ L 1 |t -s| ∀t, s ≥ t 0 0 ≤ S(t) ≤ S max ∀t ≥ t 0 ,
where:

(4.2) N max := m 0 R 1 c 1 , L 1 := 2m 0 R 1 µ max S max := m 0 2 -R 0 ∆Ω sup s∈[0,1] Ω1 Ω0 m ρ (s -a)da.
Proof. We set N r (t) = N 0 (t -t 0 ) and P (t) = P 0 (t -t 0 ) for t ∈ [t 0 -T 1 , t 0 ) and remark that the equations (3.8a)-(3.8c) allow to extend N r and S on [t 0 -T 1 , t 0 + T 0 ) in a unique and continuous way. In particular we have that N r , S ≥ 0 and N r (t) ≤ m 0 c 1 R 1 for t ∈ [t 0 , t 0 + T 0 ) and S(t) ≤ S max for t ∈ [t 0 , t 0 + T 0 ) since N m(N ) ≤ m 0 for all N . Then (3.8b) can be uniquely solved in [t 0 , t 0 + T 0 ) with respect to P with S and P (t 0 ) given. Indeed, no blow-up can occur at or before t 0 + T 0 since F (D, S) ≤ 1 for all D > 0 and S ≥ 0. An inductive argument shows that the same properties hold true on the interval [t 0 + (k -1)T 0 , t 0 + kT 0 ) for all k ∈ N.

Now, let us fix any s, t ≥ t 0 such that s ≤ t and compute

N r (t) -N r (s) = t-A0 t-A1 N r (α)m N r (α) m ρ (α)R P (α) dα - s-A0 s-A1 N r (α)m N r (α) m ρ (α)R P (α) dα = t-A0 s-A0 N r (α)m N r (α) m ρ (α)R P (α) dα - t-A1 s-A1 N r (α)m N r (α) m ρ (α)R P (α) dα ≤m 0 R 1 t-A1 s-A1 m ρ (α)dα + t-A0 s-A0 m ρ (α)dα ≤2m 0 R 1 µ max |t -s|.
Remark 4.3. In fact, equation (3.8a) prescribes the value N r (t 0 ) which may be different from N 0 (t - 0 ). Hence the solution component N r may have a jump discontinuity at t 0 even if N 0 is continuous. However, N r is bounded and Lipschitz continuous on [t 0 , +∞) with constants that do not depend on initial data. On the other hand, it is clear from the proof of Proposition 4.2 that P (t) > 0 for all t > t 0 if and only if P (t 0 ) > 0. In particular, if P 0 is not identically zero but satisfies P 0 (t 0 ) = 0, then we have P (t) = 0 for all t ≥ t 0 , no matter what is N 0 .

The next results show that, under suitable assumptions, all the components of the solution eventually are uniformly bounded away from zero. In particular, the conditions in statements (2) and (3) of Proposition 4.4 require that the breeding strategy has to be suitably tuned to the maximal fertility rate. Moreover, the obtained estimates will allow to define a compact invariant set which absorbs in finite time all the relevant solutions.

Proposition 4.4. Let (N r , P ) be the solution of (3.8a)-(3.8c) with initial data (N 0 , P 0 ) ∈ X at time t 0 .

(

) If N r (t) ≤ N max for a.a. t ∈ [ t -A 1 , t] for some t ≥ t 0 , then N r (t) ≥ m 0 R 0 c 0 2 min ess inf [ t-A1, t] N r , N 1-γ max ∀t ∈ t, t + A 0 , 1
In particular this inequality holds for all

t ≥ t 0 + A 1 by Proposition 4.2. (2) If m 0 R 0 c 0 > 2 and inf [ t-A1, t] N r ≥ N 1-γ max for some t ≥ t 0 + A 1 , then N r (t) ≥ N min and S(t) ≥ S min ∀t ≥ t,
where: (4.3)

N min := m 0 R 0 c 0 2 N 1-γ max and S min := m 0 2 -R 1 2∆Ω N 1-γ max inf s∈[0,1] Ω1 Ω0 m ρ (s -a)da.
(3) If m 0 R 0 c 0 > 2 and N 0 (a) > 0 for almost all a ∈ [-A 1 , 0], then there exists t * ≥ t 0 such that N r (t) ≥ N min and S(t) ≥ S min for all t ≥ t * .

Proof.

(1) If N r (t) ≤ N max for a.a. t ∈ [ t -A 1 , t] then (4.4) N r (t)m(N r (t)) ≥ m 0 2 min ess inf [ t-T1, t] N r , N 1-γ max for a.a. t ∈ [ t -A 1 , t]
by our assumption on m, and Statement 1 follows immediately from (3.8a).

(2) If, moreover, m 0 R 0 c 0 > 2 and inf [ t-A1, t] N r ≥ N 1-γ max for some t ≥ t 0 + A 1 , then by Statement 1 we deduce that

N r (t) ≥ m 0 R 0 c 0 2 N 1-γ max = N min ≥ N 1-γ max ∀t ∈ [ t, t + A 0 ]
and, by induction, we obtain that N r (t) ≥ N min for all t ≥ t. The inequality for S just follows from (3.8c), (4.4) and the inequality just proved for N r and Statement 2 is proved.

(3) First note that, even if ess inf N 0 = 0, we have that N r (t) > 0 for all t ∈ [t 0 , t 0 + A 0 ] by (4.1) and (3.8a). An iteration of this argument shows that N r (t) > 0 for all t ≥ t 0 .

By Statement 2, if some t ≥ t 0 + A 1 exists such that inf [ t-A1, t] N r ≥ N 1-γ max , then Statement 3 holds true with t * = t. Let us assume this does not happen and show this implies a contradiction, which will end the proof of Statement 3. In other words, we now assume that (4.5)

∀t ≥ t 0 + A 1 , 0 < inf [t-A1,t] N r < N 1-γ max .
In particular, by Statement 1 with t =

t 0 + A 1 , ∀t ∈ [t 0 + A 1 , t 0 + A 1 + A 0 ] , N r (t) ≥ m 0 R 0 c 0 2 inf [t0,t0+A1] N r > inf [t0,t0+A1] N r
and applying the same reasoning k ≥ 1 times (since equation (4.5) is assumed to hold for every t ≥ t 1 ), we get that

∀t ∈ [t 0 + A 1 + kA 0 , t 0 + A 1 + (k + 1)A 0 ] , N r (t) ≥ m 0 R 0 c 0 2 inf [t0,t0+A1] N r . Taking k ≥ A 1 /A 0 we get inf [t0+A1,t0+2A1] N r ≥ m 0 R 0 c 0 2 inf [t0,t0+A1] N r .
Since equation (4.5) is assumed to hold for every t ≥ t 0 + A 1 , we can continue applying similar estimates and show that for every ≥ 1, inf

[t0+ A1,t0+( +1)A1] N r ≥ m 0 R 0 c 0 2 inf [t0,t0+A1] N r ,
which implies that the left-hand side tends to infinity as tends to infinity since m0R0c0 2

> 1 and inf [t0,t0+A1] N r > 0 . This is in contradiction with the boundedness of N r (see Proposition 4.2).

Proposition 4.5. Let (N r , P ) be a solution of (3.8a)-(3.8c) and assume that some t * ≥ t 0 exists such that 0 < S min ≤ S(t) ≤ S max for all t ≥ t * , where we recall that S min is defined in Proposition 4.4 and S max is defined in Proposition 4.2.

(1) Let P * ≥ 0 be such that D(P * ) < S min . If P ( t) > P * for some t ≥ t * , then we have

P (t) < P * ∀t > t + P ( t) -P * λP * • D 0 + S max S min -D(P * )
.

(2) Let P * ≥ 0 be such that D(P * ) > S max . If 0 < P ( t) < P * for some t ≥ t * , then we have (1) The first statement directly follows from equation (3.8b) since F (D(P ), S) < 0 for all S ∈ [S min , S max ] and P > P max , and F (D(P ), S) > 0 for all S ∈ [S min , S max ] and P < P min (recall that S(t) ≤ S max for all t ≥ t 0 by Proposition 4.2).

P (t
(2) The second statement is a straightforward consequence of Proposition 4.2, Statement 3 of Proposition 4.4, Proposition 4.5 with (P * , P * ) = (P min /C, CP max ) and the first statement.

(3) The third statement follows from equation (3.8b) since |F (D, S)| ≤ 1 for all D, S.

As is usual in delayed equations, we study the dynamical system produced on the space X of initial conditions by the solutions of (3.8a)-(3.8c). In particular, due to the natural periodicity of the seasonality function m ρ , we consider how the initial condition is transformed after 1 year. Namely, let Π : X → X be defined by Π(N 0 , P 0 ) = (N 1 , P 1 ) with N 1 (s) = N r (1 + s) and P 1 (s) = P (1 + s) for s ∈ [-T 1 , 0], where (N r , P ) is the unique solution of (3.8a)-(3.8c) with initial data (N 0 , P 0 ) at t 0 = 0. Proposition 4.7. Π is continuous.

Proof. Let (N r , P ), (M r , Q) be the solution couples of (3.8a)-(3.8c) with initial data (N 0 , P 0 ), (M 0 , Q 0 ) ∈ X , respectively, at time t 0 = 0 and let S (N0,P0) , S (M0,Q0) be the corresponding supply functions given by (3.8c). Thanks to the continuity of m and R it is straightforward to show that |N r (t) -M r (t)| and |S (N0,P0) (t) -S (M0,Q0) (t)| can be made arbitrarily and uniformly small on [0, T 0 ] provided that (N 0 -M 0 , P 0 -Q 0 ) X is small enough. The same holds also for |P (t) -Q(t)| on [0, T 0 ] by standard results on the theory of ordinary differential equations (see [Har64, Theorem II.3.2, p. 14]). It is sufficient to iterate this procedure a finite number k of times such that kT 0 ≥ 1 to complete the proof. Now we look for an attractor for Π. The first step is to identify a compact invariant set that absorbs (almost) all orbits of the dynamical system. Its definition is suggested by the propositions we proved and goes as follows.

Definition 4.8. Let K be the set of couples (N 0 , P 0 ) ∈ X such that they are Lipschitz continuous with constants 2m 0 R 1 µ max and 2λP max , respectively (see Proposition 4.2 and Corollary 4.6), they satisfy N min ≤ N 0 (s) ≤ N max and P min /2 ≤ P 0 (s) ≤ 2P max for almost all s ∈ [-T 1 , 0] and, moreover, We cannot expect the basin of attraction of K to be the whole X , since the equations (3.8a)-(3.8c) admit semi-trivial solutions (N r , 0) and (0, P ) besides the trivial one. Moreover, even if we provide initial data (N 0 , P 0 ) ∈ X such that N 0 > 0 and P 0 ≡ 0, we have that P (t) = 0 for all t ≥ 0 if P 0 (0) = 0. Therefore we consider the following subset of X : X * = (N 0 , P 0 ) ∈ X : N 0 (a) > 0 a.e. in [-A 1 , 0] and P 0 (0) > 0 .

We obtain the next result as an immediate consequence of Ascoli-Arzelà's Theorem, Propositions 4.2 and 4.4, Corollary 4.6 and the 1-periodicity of m ρ . Proposition 4.9. Assume that m 0 R 0 c 0 > 2, D 0 > S max and S min > 0. Then K is a compact subset of X , Π(K) ⊆ K and for all (N 0 , P 0 ) ∈ X * there exists k ∈ N such that Π k (N 0 , P 0 ) ∈ K.

Using the absorbing compact set K, we can now prove the existence of an attractor for all initial conditions of X * . Theorem 4.10. Assume that

m 0 R 0 c 0 > 2 (4.7) D 0 > S max (4.8) S min > 0 (4.9)
with S max and S min given respectively in (4.2) and (4.3), and let

Λ = k∈N Π k (K) .
Then:

(1) Λ is a non-empty compact subset of X ;

(2) Π k (Λ) = Λ for all k ∈ N;

(3) for each neighbourhood U of Λ and every (N 0 , P 0 ) ∈ X * there exists k * ∈ N such that Π k (N 0 , P 0 ) ∈ U for all k ≥ k * ; (4) Λ is an attractor of the dynamical system generated by the iterates of Π and its basin of attraction contains X * .

Proof. Statements 1 and 2 follow from the continuity of Π, the compactness of K and its Π-invariance. Let (N 0 , P 0 ) ∈ X * . There exists k * such that Π k (N 0 , P 0 ) ∈ K for all k ≥ k * by Proposition 4.9. Therefore, the sequence {Π k (N 0 , P 0 )} is relatively compact, all its limit points lie in Λ and Statement 3 follows.

Let V = (Π k1 ) -1 (U )
where k 1 > T 1 is a fixed integer and

U := (N 0 , P 0 ) ∈ X : ess inf [-T1,0] N 0 > N 1-γ max and P min 3 < P 0 (s) < 3P max for all s ∈ [-T 1 , 0] .
V is an open subset of X , since U is open and Π is continuous, and contains K since Π k1 (K) ⊆ K ⊂ U . Thus, V is an open neighborhood of Λ. We claim that there exists an integer k 2 ≥ k 1 such that Π k2 (V ) ⊆ K. Let us fix (N 0 , P 0 ) ∈ V and call (N r , P ) the unique solution of (3.8a)-(3.8c) with initial data (N 0 , P 0 ) at time t 0 = 0. Moreover we define N k (s) = N r (k -s) and P k (s) = P (k -s) for s ∈ [-T 1 , 0], which means that Π k (N 0 , P 0 ) = (N k , P k ). By construction (N k1 , P k1 ) ∈ U and, therefore inf [k1-T1,k1] N r = inf [-T1,0] N k1 > N 1-γ max and, by Statement 2 of Proposition 4.4, N r (t) ≥ N min and S(t) ≥ S min for all t ≥ k 1 . In particular inf [-T1,0] N k ≥ N min for all k ≥ 2k 1 . Next we observe that P (t) ∈ (P min /3, 3P max ) for all t ≥ k 1 by Corollary 4.6 with C = 3. The estimates of Proposition 4.5 grants that, if P (k 1 ) ∈ [P min /2, 2P max ], then P (t) enters (and thereafter remains in) the interval [P min /2, 2P max ] at a time which is bounded above by

t 2 = k 1 + max D 0 + S max 2λ[S min -D(2P max )] , D 0 + S max 2λ[D(P min /2) -S max ] ,
a number that is independent of (N 0 , P 0 ) ∈ V . The claim follows by the choice

k 2 ≥ max{2k 1 , t 2 + T 1 }. Therefore, we have that Π k2 (V ) ⊆ K ⊂ V , k∈N Π k (V ) = Λ
and Λ is an attractor for the dynamical system generated by Π.

Since the coefficient m ρ is a 1-periodic function, we argue now about the existence of k-periodic solutions of (3.8a)-(3.8c), i.e. fixed points of Π k , k ∈ N. The analysis carried on till now allows to prove it as a straightforward application of Schauder fixed point theorem at least for k ≥ T 1 . If we consider the set D ⊂ X made of the couples (N 0 , P 0 ) such that they satisfy the same conditions that define K with the only exception of (4.6), we observe that D is a compact and convex subset of X (convexity is the reason why in the definition of D we have removed condition (4.6)). Moreover, using the same notation introduced in the proof of Theorem 4.10, N k is a translation of the restriction of N r to the interval [k -T 1 , k] and, hence, it is still continuous since k -T 1 ≥ 0 (see Remark 4.3). Thus, D is Π k -invariant by Propositions 4.2 and 4.4 and Corollary 4.6 and Schauder fixed point theorem applies to Π k : D → D.

On the other hand, if k < T 1 , this argument fails: indeed, N k may now have a jump discontinuity at time -k, since (4.6) is no more guaranteed, and D is no more Π k -invariant. However, the points of discontinuity may only appear at negative integer times and may only be of jump type, while between two successive integer times the component N k remains Lipschitz continuous with the same known constant. This remark suggests in which way the definition of the set D should be modified in order to preserve its Π k -invariance when k < T 1 . We detail the argument in the following result when k = 1 < T 1 and we remark that we did not pursue the question of finding the minimal period of k-periodic solutions: the fixed point of Π k obtained above might actually be a fixed point of Π.

Theorem 4.11. Assume that (4.7)-(4.9) hold true. Then equations (3.8a)-(3.8c) have a non-trivial 1-periodic solution.

Proof. We will consider only the case T 1 > 1 since we already explained how to deal with the case T 1 ≤ 1. Let i 0 ≥ 1 be the largest integer such that i 0 < T 1 and let J i = (-i, -i + 1), for i = 1, . . . , i 0 , while we let J i0+1 = (-T 1 , -i 0 ). Moreover, let D be made by the couples (N 0 , P 0 ) ∈ X such that:

• N 0 is Lipschitz continuous on each open interval J i , for i = 1, . . . , i 0 + 1, with constant 2m 0 R 1 µ max ; • P 0 is Lipschitz continuous with constant 2λP max ;

• N 0 (s) ∈ [N min , N max ] and P 0 (s) ∈ [P min /2, 2P max ] for all s ∈ [-T 1 , 0]. We note that the assumptions on N 0 ensure that the one sided limits N 0 (-i + ) and N 0 (-i -) exist and are finite for all i = 1, . . . , i 0 .

The set D is clearly convex and let us consider a sequence {(N 0,n , P 0,n )} n∈N ⊂ D in order to prove that D is compact in X . By Ascoli-Arzelà's Theorem it is possible to extract a subsequence, that we still call {(N 0,n , P 0,n )} n∈N , such that P 0,n and N 0,n | Ji converge uniformly for i = 1, . . . , i 0 + 1. It is easy to show that the convergence is with respect to the norm of X and that the limit belongs to D, which, thus, is compact. Now, by the very definition of Π, we have that, if (N r , P ) is the solution of (3.8a)-(3.8c) with initial data (N 0 , P 0 ) ∈ D at time t 0 = 0 and (N 1 , P 1 ) := Π(N 0 , P 0 ), then N 1 (s) = N 0 (s + 1) for s ∈ [-T 1 , -1), while N 1 (s) = N r (s + 1) for s ∈ [-1, 0]. Therefore, N 1 may have jump discontinuities only at i = -i 0 , . . . , -1. Propositions 4.2 and 4.4 and Corollary 4.6 (with C = 2) grant that Π(N 0 , P 0 ) ∈ D and Schauder fixed point theorem concludes the proof.

Results of simulation experiments

This section shows numerical experiments for the model described in Section 3. 5.1. Parameters of the model. We discuss here in which way one could choose some of the model's parameters in order to take into consideration real-world cattle, mainly hog. However, we recall that in our numerical experiments our choice will be to remain close to the values used in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] for Microtus Epiroticus (see Subsection 5.1.5), since our primary goal here is to check the effects that the interaction with the market equation has on the dynamics of Yoccoz-Birkeland equation (and viceversa). 5.1.1. Biological parameters. For the common pork the following choices are supported from the literature [START_REF] Marinus | Animal Reproduction, Principles and Practices[END_REF] as well as from discussion with farmers [Zan09].

In our experiments, similarly to [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF], we take the fertility function defined by

m(N ) = m 0 max{N, 1} -γ ,
where γ and m 0 remain to be chosen. The parameter m 0 should count the largest number of female pups (approximately one half of the total number of pups) a fertile female can give birth to in a year. For common pork it may reasonably range in the interval [10, 13] [KH13]. On the other hand, when one considers races of pork, like Cinta Senese in Italy, which are bred in a semi-wild environment, a lower choice of m 0 should be chosen. Assuming a maximum of two litters per year [START_REF] Campodoni | I sistemi di allevamento. Il bosco e l'allevamento della razza Cinta Senese[END_REF] and the data provided for instance in [START_REF] Franci | Italian autochthonous pigs: progress report and research perspectives[END_REF], the value of m 0 for Cinta Senese should range in the interval [5, 10].

The following table is taken from [START_REF] Marinus | Animal Reproduction, Principles and Practices[END_REF] and could be used to choose values for A 0 depending on the species one is considering: Species age at puberty (months It is more difficult to suggest a choice for A 1 since the maximal age of reproducing cattle is a quantity that depends heavily on the breeding strategies of the breeder. In our setting, A 0 and A 1 are taken as in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF]. 5.1.2. Seasonality. We assume that the births are synchronized [Zan09] and we consider the following 1-periodic piecewise-constant seasonal factor:

)
m ρ (t) = 1 1 -ρ 1 t-t ∈[0,1-ρ) .
5.1.3. Demand function. In the numerical experiments we perform in this section, we choose the following demand function for some parameters D 0 , α D > 0. For simplicity we took the parameters Ω 0 and Ω 1 in the supply function to be equal to A 0 and A 1 , respectively. As we have already noted the normalization of the supply is ∆Ω = Ω 1 -Ω 0 .

5.1.4. Breeder strategy. The function R defines how the breeder takes into account the current meat price for deciding how to split newborn females among the reproduction line (for a long-term increase of the supply) and the butchery line (for a short-term increase of the supply, but still delayed). A short term strategy would be to take R close to 0 when the price is high. A long term strategy is to take R close to 1 when the price is high. In the numerical experiments we choose (5.2)

R(P ) = R logistic (P ) := R 0 + (R 1 -R 0 )f d (P/P 0 )
where

f d (x) = x d 2 if x ∈ [0, 1) 1 1+exp( -2d(x-1) )
otherwise and R 0 , R 1 ∈ [0, 1] and P 0 , d > 0 have to be chosen. 5.1.5. Parameter setting. The default set of parameters (called SP, because it is close to realistic values for pork, see Section 5.1.1) and functions is the following:

• population dynamics: A 0 = 0.18, A 1 = 2, m defined by (2.3) with m 0 = 5 and γ = 8.25, and the seasonality function is

m ρ (t) = 1 1 -ρ 1 t-t ∈[0,1-ρ)
with ρ = 0.79. • market dynamics: λ = 1 and the demand function is D = D exp is defined by equation (5.1) with D 0 = 5 and α D = 1. • interaction between population and market: Ω 0 = 0.18, Ω 1 = 2, and R = R logistic is defined by equation (5.2) with R 0 = 0 (minimal value), R 1 = 1 (maximal value), P 0 = 1 (price threshold) and d = 4 ("degree of R(P ) for small P "). The setting SP was chosen close to the main setting studied in [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF] (up to the changes in the model), with m 0 one order of magnitude smaller (i.e., 5 instead of 50, which is more realistic, and necessary to obtain a reasonable behaviour because of the change in the survival function) and ρ modified (the latter choice resulted from a rough exploration of the main parameters of the model). We remark that the choice of parameters in SP is not far from realistic ones for livestock production.

In our numerical experiments, only a couple of parameters will be varying: m 0 and γ, which are the ones that mostly influence the population dynamics. In particular, figures 8 and 9 show that the model has a chaotic behavior also for values of γ close to 2, i.e. a much weaker dependence of fertility on over-population.

Study of one parameter set (setting SP).

In this section, we focus on the setting SP, for which some interesting behaviour can be observed. 5.2.1. Continuous time dynamics. The continuous time dynamics of setting SP can be visualized on Figure 1 over 50 years. Note that on Figure 1 the mature reproducing population N r (t) goes through very small values, of order 10 -3 , which might seem unrealistic. This phenomenon can be interpreted as follows. There is no mortality, except for animals reaching the maximal ages A 1 and Ω 1 . So, if most of the mature reproducing population N r was born during a short period of time, N r (t) will drop down to (almost) zero A 1 years later. Here, we observe such drops, which do not endanger the whole population because it holds during the reproducing season and it doesn't occur too fast (the reproducing females were not all born during a too short time period). Therefore, once the drop has started, as soon as N r goes below 1, a large number of birth will happen, leading to mature reproducing females after a time delay A 0 . So, even if N r was going exactly to zero during this delay period, it would increase again as soon as the newborn females become mature.

A good way to visualize this phenomenon is given by Figure 2 where N r (t) is plotted (left) together with the totale female population in the reproducing line (right): the latter quantity never goes below 0.4.

A similar phenomenon holds with the butchery line (see Figure 3), where the "mature" butchery population (proportional to the supply S(t)) goes close to zero (left) but not the total butchery population (right).

Yearly dynamics.

A natural way to display the behavior of the dynamical system we are studying is to only look at the values it takes for t ∈ N, that is, at the very beginning of the birth period. We recall that the continuous time dynamics takes place in an infinite-dimensional space X corresponding to couples of functions (N (t), P (t)), where t varies in some interval. Our discretization (see Appendix A) leads to a discrete dynamical system in a phase space of dimension 2 × 201. Then, we can visualize the dynamics by plotting in R 3 the set { (N r (t), N r (t + 1), N r (t + 2)) , t ∈ N }, as shown by Figure 4, or the set { (P (t), P (t + 1), P (t + 2)) , t ∈ N }, as shown by Figure 5. Both figures are the projection of the same attractor on two different subspaces. We estimate the fractal dimension of the first set to 1.52 and of the second one to 1.84. It seems reasonable to conjecture the existence of a strange attractor of dimension d ∈ (1, 2). See Section A.4, Figure 12 and Figure 13, for details on how we estimate the fractal dimension. Note that Figure 4 shows an attractor similar to the one of the Yoccoz-Birkeland model with ρ = 0.30, γ = 8.25, A 0 = 0.18 and m 0 = 50 [Arl04, Figure 12]; only the center of the attractor here seems more complex. On the other hand, the interaction between the population and the price equation is crucial, as shown by Section 5.2.3 below.

5.2.3.

Comparison with setting SP 1 : SP with R constant. Let us call SP 1 the setting SP with the value of the fraction of reproducing females in (3.8a) freezed: R(P ) = R cst ≈ 0.955 for all P > 0. The constant value taken for R is equal to the empirical average of R(P (t)) in setting SP for t ∈ [290 000, 300 000]. So, comparing results obtained with SP 1 and SP shows the effect of the coupling between price and population. Note also that setting SP 1 then is an instance of the Yoccoz-Birkeland Figure 4. 3d plot of (N r (t), N r (t + 1), N r (t + 2)) with t ∈ N, 280 000 ≤ t ≤ 300 000 for setting SP. Figure 5. 3d plot of (P (t), P (t + 1), P (t + 2)), with t ∈ N, 280 000 ≤ t ≤ 300 000 for setting SP. model with A 0 = 0.18, A 1 = 2, γ = 8.25 and m 0 ≈ 4.78 (with slightly different functions m and m ρ , compared to [START_REF] Arlot | Étude d'un modèle de dynamique des populations[END_REF]).

The yearly dynamics of setting SP 1 can be visualized on Figure 6 for 95 000 ≤ t ≤ 100 000: it exihibits only a low-complexity seemingly nonchaotic orbit, very close to being periodic and totally different also from the original Yoccoz-Birkeland attractor. Figure 6. 3d plot of (N r (t), N r (t + 1), N r (t + 2)), with t ∈ N, 95 000 ≤ t ≤ 100 000 for setting SP 1 . 5.2.4. Analysis of chaos. In order to analyze the chaoticity of the dynamics of the price P , we followed a time series approach as in [START_REF] Kantz | Nonlinear Time Series Analysis[END_REF]: we computed the autocorrelation function and determined its first zero τ * . Then we sampled the price time series at time steps multiple of τ * and we computed the combinatorial entropy for the binary sequence obtained by looking at the sign of price returns (see Figure 7).

Given a times series (Y (t)) t≥0 , its autocorrelation function is defined by

∀τ > 0 , R Y (τ ) := Y -Y Y (• + τ ) -Y (• + τ ) Y -Y 2 Y (• + τ ) -Y (• + τ ) 2 ,
where • means an average over time. In other words, the autocorrelation R Y (τ ) is the correlation between Y (t) and Y (t + τ ) for a random t chosen uniformly in [0, T ], for some large time window T > 0. So, roughly speaking, R Y (τ ) close to 0 means that Y (t) does not provide information for predicting Y (t + τ ). The autocorrelation function of (P (t)) t≥0 is plotted on the top left of Figure 7. Its absolute value tends to be smaller for larger values of τ , and it crosses zero several times, first for τ = τ ≈ 1.37. Therefore, the discrete dynamical system (P (kτ )) k∈N is a good candidate for being unpredictable. So, we consider the corresponding "returns" r k := log 10 P (k + 1)τ -log 10 P (kτ ) (which are plotted on the top middle graph of Figure 7), and evaluate the combinatorial entropy H K of the binary sequence ((sign(r i )) k-K+1≤i≤k ) k∈N for various values of K 1 .

1 The procedure we follow provides a lower bound to the Kolmogorov-Sinai entropy of the flow. Indeed, this is defined as the supremum over all finite partitions of the rate of change of entropy due to the finer partitioning given by the flow at each time step. Here we have just fixed The top right plot of Figure 7 shows H K as a function of K: there is a clear linear trend, with a positive slope 0.611 (correlation coefficient 0.995), which suggests a positive entropy.

We also tried the same analysis for the discrete dynamical system (P (k)) k∈N , that is, the time series at integer times, for which the autocorrelation function is not estimated to be (almost) zero. The results are given on the bottom part of Figure 7, with the corresponding returns (middle) and entropy H K as a function of K (right): linear regression still suggests a linear trend with a positive slope 0.647 (correlation coefficient 0.983), even if the returns seem less unpredictable as with τ time steps. In particular, chaotic behavior seems to arise also for small values of γ (i.e. when the fertility rate m(N ) depends in a weaker way on the population N ) which are probably more reasonable in a human-controlled breeding facility. reproducing ages of females. So, for instance, by Eq. (A.1), the reproducing female population at time k A1 + 1 is the average of random variables uniform over [0, 2], so it should be close to 1, the threshold value for density-dependence to apply. A.3. Chaos analysis. In each setting and for each time series Y ∈ { N r , P }, we start from the "continuous" time series of Y (t) for t > T max -10 000, where T max is the total length of the simulation experiment (i.e., T max = 300 000 for setting SP, and T max = 100 000 for other settings). Then, we compute the empirical autocorrelation function of Y with lags τ ∈ [0, 100] (τ is expressed in years). We get its first "zero" τ as the first point where the empirical autocorrelation crosses zero, and we then check that the value of the empirical autocorrelation at τ = τ is smaller than an arbitrary threshold (here, 10 -2 ); this holds true in all the settings where we did such an analysis.

Then, we consider the sequence Y (t) only at times t equal to an integer multiple of τ , and compute the returns log 10 (Y ((i+1)τ )/Y (iτ )) and their respective signs ε i ∈ { -1, 1 }. From this (finite) sequence of signs, we compute for K = 1, . . . , 12 the combinatorial entropy H K (τ ) of ((ε i ) k-K+1≤i≤k ) k∈N :

H K (τ ) := - x∈{ -1,1 } K p x log 2 (p x )
where p x = P (ε i ) k-K+1≤i≤k = x the latter probability being with respect to k. Finally, we plot H K (τ ) as a function of K and perform a standard robust linear regression in order to estimate its slope.

We also do the same analysis when considering Y (t) at integer times, leading to a plot of H K (1) as a function of K, on which we estimate its slope by performing a standard robust linear regression.

A.4. Fractal dimension. For computing the fractal dimension of the attractors for setting SP, we started from their 3-dimensional visualization over 200 000 years,

m

  ρ (t -a)da > 0, ensures that D -1 ([S min , S max ]) is a compact interval.

N 0

 0 (-a)m N 0 (-a) m ρ (-a)D P 0 (-a) da .

  ) = D exp (P ) := D 0 e -α D P

Figure 1 .

 1 Figure 1. Continuous time dynamics for setting SP over 50 years: N r (top left), S and D (top right), B r (middle left), P (middle right), B b (bottom left) and R(P ) (bottom right).

Figure 2 .

 2 Figure 2. Setting SP, continuous time dynamics of reproducing female population over 20 years: mature (N r , left) and total population (right).

Figure 3 .

 3 Figure 3. Setting SP, continuous time dynamics of butchery population over 20 years: "mature" (∝ S, left) and total population (right).

Figure 7 .

 7 Figure 7. Analysis of chaos for P in setting SP. Top left: autocorrelation function. Its first (approximate) zero is at τ = 1.37. Middle: sequence of the returns r k (τ ) = log 10 P ((k +1)τ )-log 10 P (kτ ) for τ = τ (top) and τ = 1 (bottom). Right: estimated entropy H K as a function of K, with τ = τ (top; slope 0.611 and correlation coefficient 0.995) and τ = 1 (bottom; slope 0.647 and correlation coefficient 0.983). See Section 5.2.4 and Appendix A.3 for details.

5. 3 .

 3 Exploration of the parameter space. Starting from setting SP, we changed the parameters one by one between γ and m 0 . The resulting bifurcation diagrams are shown on Figures 8 and 10, for N r and Figures 9 and 11, for P . Appendix A.5 provides details about how the bifurcation diagrams have been obtained. 5.3.1. Setting SP with γ varying. A bifurcation diagram is plot on Figure 8, showing { N r (t), t ∈ N } as a function of γ ∈ [2, 10]. A corresponding bifurcation diagram for the price P is shown on Figure 9. The dynamics look like complex for the largest part of the interval, with only few and very small windows of periodicity.

  5.3.2. Setting SP with m 0 varying. A bifurcation diagram is plot on Figure 10, showing { N r (t), t ∈ N } as a function of m 0 ∈ [2, 8]. The corresponding bifurcation diagram for the price P is shown on Figure 11. Here period doubling cascades seem to arise in the windows of periodicity. Moreover, at a closer inspection of a particular type of partition (corresponding to the choice of sign of the returns). The use of τ * instead of 1 as a time unit only affects entropy by a multiplicative factor.

Figure 10 .Figure 11 .

 1011 Figure 10. Bifurcation diagram for N r (t) w.r.t. m 0 ∈ [2, 8] around setting SP.

Figure 12 .

 12 Figure 12. Estimation of the fractal dimension of the attractor of Figure 4 (setting SP, N r ).

  ) -S min D 0 + S max and, thus,P (t) ≤ -λP * S min -D(P * ) D 0 + S max (t -t) + P ( t) .

	Therefore, P (t) reaches the level P * before the time
	t +	P ( t) -P * λP *	•	D 0 + S max S min -D(P * )	,
	afterwards P (t) remains below P * since P (t) = P * implies P (t) ≤ 0 by equa-
	tion (3.8b), which proves Statement 1.
	(2) Statement 2 follows in a similar way once we observe that, as long as P (t) ≤
	P

) > P * ∀t > t + P * -P ( t) λP ( t) • D 0 + S max D(P * ) -S max . Proof. (1) As long as P (t) ≥ P * , D(P ) ≤ D(P * ) < S min so that by equation (3.8b) P decreases, P (t) ≤ λP * D(P * * , S max > D(P * ) ≥ D(P ) so that by equation (3.8b) P increases and P (t) ≥ λP ( t) D(P * ) -S max D 0 + S max . Corollary 4.6. Assume that m 0 R 0 c 0 > 2, D 0 > S max and S min > 0 and let C > 1 and P min , P max be such that [P min , P max ] = D -1 ([S min , S max ]). (1) If S(t) ≥ S min for all t ≥ t * and P ( t) ∈ [P min /C, CP max ]) for some t ≥ t * , then P (t) ∈ [P min /C, CP max ]) for all t ≥ t. (2) In any case, for every non trivial solution (N r , P ) of (3.8a)-(3.8c) there exists t such that P (t) ∈ [P min /C, CP max ]) for all t ≥ t. (3) Moreover, as long as P (t) stays in [P min /C, CP max ]), we have that |P (t)| ≤ CλP max . Proof. Remark that the assumption S min > 0, which is equivalent to requiring that R 1 < 2 and inf t∈[0,1]

This section provides all technical details necessary to reproduce our numerical experiments.

A.1. Discretization of the model. Throughout the paper, we consider the following discretization with q = 100 steps per year. This will lead to a deterministic dynamical system in a phase space of dimension 2 × 201.

A.1.1. Notation and parameters. We replace the continous time parameter t ∈ [0, +∞) by indices i ∈ N\ { 0 }.

Roughly, i ≥ 1 replaces the interval i-1 q , i q where q ≥ 1 is an integer (number of steps per year).

• q ∈ N\ { 0 } is the number of steps per year • N r,i approximates the average of N r (t) over t ∈ i-1 q , i q :

N r,i ≈ q i q i-1 q N r (t)dt

• Similarly, N b,i , m ρ,i , S i and P i are respectively defined from N b (t), m ρ (t), S(t) and P (t).

• B r,i approximates the number of newborn females put in the reproducing line in the interval t ∈ i-1 q , i q :

Note the difference with N r coming from the fact B r (t) is a density (number per unit of time) contrary to N r (t) that is a number of individuals.

We also need to discretize some of the model parameters (with a few technical modifications in order to avoid some troubles):

is the closest integer to x. We want k A0 > 0 so that N r,k only depends on the past of B r in (A.1), otherwise we would have a circular definition.

because animals are supposed to die exactly at age A 1 , so individual of age

in the following order:

A.2. Initial condition. The initial condition is given by the values of the number of newborns (B r,i and B b,i ) for time steps i = 1, . . . , k A1 before the start of the simulation, since the discretized model above (Eq. (A.1)-(A.6)) only needs the past values of B r and B b to compute the future dynamics of the model.

For each parameter set, we choose an initial condition randomly (but with a known random seed, so we can use repeatedly the same initial condition to obtain our bifurcation diagrams, for instance) as follows: B r,1 , . . . , B r,k A 1 , B b,1 , . . . , B b,k A 1 are chosen independent with common distribution U 0,

. The reason for this choice is that k A1 -k A0 +1 is the number of time steps corresponding to that is, for N r for instance: K = N r (t), N r (t + 1), N r (t + 2) , 100 000 ≤ t ≤ 300 000 .

Then, for various values of > 0, we compute the number N (K) of cubes C i,j,k = [i ; (i + 1) ] × [j ; (j + 1) ] × [k ; (k + 1) ] that contain at least one point of K. Figure 12 thus represents log 10 N (K) as a function of log 10 .

Theoretically, the fractal dimension is the opposite of the slope of this graph at infinity. Here, since K is finite, N (K) is constant equal to Card(K) for small . So, we estimated the slope of the graph only for a limited set of values of . One can check on Figure 12 that the linear fit seems reasonably close to the original curve.

Because of the numerous approximations made during this estimation, the precise value of the fractal dimension should not be taken into account too seriously, but its order of magnitude should be correct. Figure 13 shows the corresponding estimation with P instead of N r .

A.5. Bifurcation diagrams. We took the same initial condition for all settings in all diagrams (by choosing a fixed arbitrary seed to the random number generator). Then, for each value of the parameter we plot the sequence of values of N r (t) (resp. P (t)) for t ∈ { 1 500, . . . , 2 000 }. In all diagrams, the parameters are varying by step of 0.01.


Cheng-Wen Wu

SOC Testing Methodology and Practice

published or not. The documents may come

Introduction

It is generally agreed that, for an SOC, the design and test engineers usually have to test the cores with very limited knowledge of the core test information. The issues of core access and isolation are being addressed by the IEEE P1500 standard [START_REF] Marinissen | On using IEEE P1500 SECT for test plug-n-play[END_REF]. Although the standard provides unified core test access methods, the test controller, test architecture, test access mechanism (TAM), and test integration are left to the SOC integrator. In [START_REF] Cheng | An SOC test integration platform and its industrial realization[END_REF], we stress two major issues in practical SOC test integration: 1) session-based test scheduling considering not only the realistic test control architecture and TAM bus, but also test IO limit; and 2) the coexistence of scan test and functional test for logic cores. A commercial digital still camera (DSC) controller SOC has been developed using our test integration platform-SOC Test Aid Console (STEAC). Results from the DSC controller chip show that our approach is more effective than others that use non-session-based test scheduling. Also, both the scan and functional tests are supported. All the tens of embedded memories are tested by the built-in self-test (BIST) circuits generated by our memory BIST compiler, BRAINS [START_REF] Cheng | BRAINS: A BIST complier for embedded memories[END_REF].

SOC Test Aid Console (STEAC)

Figure 1 shows the SOC test integration system called STEAC-SOC Test Aid Console [START_REF] Wang | A test access control and test integration system for system-on-chip[END_REF], which consists of four modules: the STIL Parser, Core Test Scheduler, Test Insertion Tool, and Pattern Translator. The STIL Parser parses the test information of each IP. The test information is written in STIL and is generated by commercial ATPG tools. Therefore, STEAC can be integrated into a typical design flow easily. The test information includes the IO ports, scan structure (number of scan chains, length of each scan chain, etc.), and test vectors. With these core test information, Core Test Scheduler will schedule the core tests to reduce the overall test time. The Scheduler partitions core tests into several test sessions, and assigns the TAM wires to each core to meet the power and IO resource constraints. If the IP is a soft core, the scan chains can be reconfigured. The Core Test Scheduler will then rebalance scan chains for each assigned TAM width. The results can be fed back to the SOC integrator to reconfigure the scan chains to balance the chain length. The scheduling results are also used to generate the Test Controller, TAM bus, and Test Wrapper. Finally, the generated test circuitry is inserted into the original SOC netlist automatically. A new SOC design with DFT will be ready in minutes. The core test patterns are generated at the core level. After the cores are wrapped, the test patterns must be translated to the wrapper level and then to the chip level. The test patterns are cycle based, which can be applied by external ATE easily. In Fig. 2 [5] we show our memory BIST architecture, which supports testing of plural heterogeneous memory cores. The tester can access all the on-chip memories via a single shared BIST Controller, while one or more Sequencers can be used to generate March-based test algorithms. Each Test Pattern Generator (TPG) attached to the memory will translate the March-based test commands to the respective RAM signals. With our automatic memory BIST generation system, BRAINS [START_REF] Cheng | BRAINS: A BIST complier for embedded memories[END_REF], one can generate the BIST circuit using the GUI or command shell, and evaluate the memory test efficiency among different designs easily. Moreover, BRAINS can be integrated with a memory compiler to deliver BISTed memory cores.

Experimental Results

A DSC test chip has been implemented and fabricated to verify the proposed approach. This test chip is implemented with a standard 0.25µm CMOS technology. The major digital part of the chip includes a processor, JPEG codec, TV encoder, USB, external memory interface, and tens of singleport and two-port synchronous SRAMs with different sizes. Figure 3 [START_REF] Cheng | An SOC test integration platform and its industrial realization[END_REF] gives the block diagram of this test chip. The details are given in a companion paper [START_REF] Chen | Integration, verification and layout of a complex multimedia SOC[END_REF]. 1 [START_REF] Cheng | An SOC test integration platform and its industrial realization[END_REF]. In the scan-test mode, the USB core has 2 scan vectors that are enabled by SE during Pulse-Clock, while the TV encoder has only one. When the test IO resource constraint is considered, parallel testing may not be better than serial testing. This is because more test control IOs are needed for parallel testing, so fewer IO pins can be used as the test data IOs (i.e., TAM IOs). Since there are also cases when parallel testing leads to shorter test time than serial testing, it is important to take chip IO pins into consideration so far as test time evaluation is concerned. In the DSC case, we tried several scheduling approaches, and found that the session-based approach (with three test sessions) has the shortest total test time-

Conclusions

We have presented an SOC test integration platform, with complete solutions for real problems in test scheduling, test IO reduction, functional test, scan IO sharing, embedded memory BIST, etc. The fabricated test chip has been verified, and the test design has been successfully implemented on a commercial DSC chip. Test results justify that short test integration cost, short test time, and small area overhead can be achieved. This work has been the result of the effort of many people, including my students who have coauthored the papers with me, as listed in the references.

Figure 1 :

 1 Figure 1: Test integration flow of STEAC [4].

Figure 2 :

 2 Figure 2: BIST architecture for multiple memory cores [5].

Figure 3 :

 3 Figure 3: Block diagram of the test chip [2]. The IPs to be wrapped in this test chip include the USB, TV encoder, and JPEG cores. The USB core has 4 clock domains, 3 reset signals, 1 scan enable (SE) signal, and 6 test signals. There are 4 scan chains with dedicated scan input and output for each clock domain. The TV encoder has both scan and functional tests. The test pins include one clock, reset, SE, and test enable signals. There are two scan chains in the TV encoder, where one scan chain shares the output with a functional output. The legacy JPEG core has only functional patterns and one clock domain. The clock signals for the IPs are generated by an internal PLL. The detailed test information of the IPs are shown in Table1[START_REF] Cheng | An SOC test integration platform and its industrial realization[END_REF]. In the scan-test mode, the USB core has 2 scan vectors that are enabled by SE during Pulse-Clock, while the TV encoder has only one.

  4,371,194 clock cycles as opposed to 4,713,935 cycles by non-session-based approach. The total test IOs of the three large cores are 19, including 6 clock signals, 4 reset signals, 7 test enable signals, and 2 SE signals. With shared test IOs, the test control IO counts are reduced. With STEAC, the Test Wrappers, TAM, and Test Controller have been automatically generated and inserted into the original test chip design in 5 minutes, using a SUN Blade 1000 workstation with dual 750MHz processors and 2GB RAM. The area of the WBR cell is equivalent to 26 two-input NAND gates. The Test Controller and TAM multiplexer require about 371 and 132 gates, respectively-their hardware overhead is only about 0.3%. Again, all the embedded memories are tested with the BIST circuits generated by BRAINS, which has been integrated into STEAC, as shown in Fig. 4.

Figure 4 :

 4 Figure 4: Integration of BRAINS into STEAC.

Table 1 :

 1 Test information of the cores[START_REF] Cheng | An SOC test integration platform and its industrial realization[END_REF].

	Core TI TO PI PO Scan chains (Lengths) Patterns
			(Type)
	USB 18 4 221 104 4 (1,629, 78, 293, 45) 716 (Scan)
	TV	6 1 25 40 2 (577, 576)	229 (Scan) 202,673 (Func.)
	JPEG 1 0 165 104 No scan	235,696 (Func.)

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Ilia Polian

Alejandro Czutro

Bernd Becker

Evolutionary Optimization in Code-Based Test Compression

Keywords: Test compression, code-based compression, evolutionary algorithms 1. L MVs v (1), v (2)

We provide a general formulation for the code-based test compression problem with fixed-length input blocks and propose a solution approach based on Evolutionary Algorithms. In contrast to existing code-based methods, we allow unspecified values in matching vectors, which allows encoding of arbitrary test sets using a relatively small number of codewords. Experimental results for both stuck-at and path delay fault test sets for ISCAS circuits demonstrate an improvement compared to existing techniques.

Introduction

Test data compression becomes increasingly popular for distributing test complexity between automatic test equipment and on-chip structures. By transferring compressed information from the tester and decompressing it on-chip, the requirements on IC's pin throughput can be significantly reduced while on-chip hardware cost is kept in check. Most test compression techniques proposed so far are based either on codes or on reseeding.

Code-based approaches assign a codeword C(v) to a sequence of test data v (an input block). For instance, v could be 10011101, and C(v) could be 110. Then, if 110 is transmitted from the tester, an on-chip decoder would generate the sequence 10011101 and apply it to the circuit under test. Schemes based on run-length codes [START_REF] Jas | Test vector decompression using cyclical scan chains and its application to testing core-based designs[END_REF], selective Huffman codes [START_REF] Jas | Scan vector compression/decompression using statistical coding[END_REF], Golomb codes [START_REF] Chandra | Test data compression for system-on-a-chip using golomb codes[END_REF], frequency directed codes [START_REF] Chandra | Frequency-directed run length codes (FDR) codes with application to system-on-achip test data compression[END_REF], VIHC codes [START_REF] Gonciari | Improving compression ratio, area overhead, and test application time for system-on-a-chip test data compression/decompression[END_REF], LZ77 [START_REF] Wolff | Multiscan-based test compression and hardware decompression using lz77[END_REF], mutation codes [START_REF] Reda | Reducing test application time through test data mutation encoding[END_REF], packet-based codes [START_REF] Khoche | Packet-based input test data compression techniques[END_REF] and non-linear combinational codes [START_REF] Reddy | On test data volume reduction for multiple scan chain designs[END_REF] have been proposed.

Reseeding is a technique that is applied to a pseudorandom pattern generator in order to omit those parts of its sequence that are not needed for testing. If the control information (on which patterns to omit) is fed by the tester, reseeding becomes a test compression method. Pseudorandom pattern generators used in reseeding architectures so far include LFSRs [START_REF] Krishna | Test vector encoding using partial LFSR reseeding[END_REF][START_REF] Zacharia | Decompression of test data using using variable-length seed LFSRs[END_REF], multiple-polynomial LFSRs [START_REF] Hellebrand | Built-in test for circuits with scan based on reseeding of multiple-polynomial linear feedback shift register[END_REF][START_REF] Hellebrand | Generation of vector patterns through reseeding of multiplepolynomial linear feedback shift registers[END_REF], twisted-ring counters [START_REF] Chakrabarty | Built-in test pattern generation for high performance circuits using twistedring counters[END_REF] and folding counters [START_REF] Hellebrand | A mixedmode BIST scheme based on reseeding of folding counters[END_REF].

Further test compression techniques include Embedded Deterministic Test [START_REF] Rajski | Embedded deterministic test[END_REF], which is based on a ring generator unfolding low-bandwidth test data from the tester and feeding them into the scan chains of the circuit under test. It has been proposed to use a reconfigurable switch [START_REF] Tang | On reducing test data volume and test application time for multiple scan chain designs[END_REF] or a reconfigurable interconnection network [START_REF] Li | Deterministic BIST based on a reconfigurable interconnection network[END_REF] for this purpose. A test compression architecture based on reconfiguring some of the scan flip flops into pseudo-random pattern generators is presented in [START_REF] Krishna | 3-stage variable length continuous-flow scan vector decompression scheme[END_REF].

In general, reseeding-based solutions embed test vectors into a pseudo-random sequence, which makes them attractive for testing combinational circuits. However, their applicability to sequential (non-scan), delay or functional test sets is limited. Embedded deterministic test is also designed for full scan-based circuits. In contrast, code-based compression techniques precisely reproduce the original encoded test set (possibly specifying the don't-care bits); they do not reorder the test set or add new vectors to it. Hence, code-based methods have applications for which they are inherently better suited than other compression techniques, irrespective of the achieved compression rates.

Recently, Tehranipour et al. [START_REF] Tehranipour | Nine-coded compression technique with application to reduced pin-count testing and flexible on-chip decompression[END_REF] have proposed a fixedlength input block compression technique based on exactly nine codes (9C compression). The codes are defined with respect to an even number K which is the length of an input block. For instance, let K be 6. Then, the test set is partitioned into 6-bit long pieces, i.e. input blocks. Every input block is mapped to one of the following nine matching vectors (MVs): v (1) = 000 000; v (2) = 111 111; v (3) = 000 111; v (4) = 111 000; v (5) = 111 UUU; v (6) = UUU 111; v (7) = 000 UUU; v (8) = UUU 000; v (9) = UUU UUU. Every vector v (i) is assigned a codeword C(v (i) ). For example, if 111000 appears in the test sequence, it will be encoded by C(v (4) ).

The symbol U denotes an unspecified value. For instance, 111100 and 111011 both match v (5) (as do other 6 input blocks). In order to distinguish between these 8 possible input blocks, the fill values for the U positions are spec-ified explicitly after the codeword. Thus, the input block 111100 will be coded as C(v (5) )100, and 111011 will be coded as C(v (5) )011. It is possible to represent 111000 both as C(v (4) ) and as C(v (5) )000, although the first code will probably be shorter. It is also possible to encode 111000 as C(v (8) )111 or as C(v (9) )111000. Note that any arbitrary input block can be encoded using v (9) . However, the number of U s in an MV corresponds to the number of fill values that have to be applied after the codeword for that vector; hence, it is better to use MVs with as few U values as possible.

In this paper, we generalize the use of the unspecified value U in fixed-length input block code-based test compression. We allow MVs with 1, 0 and U values on arbitrary positions instead of the fixed MV set v (1) , . . . , v (9) used in 9C compression. We also do not require the number of MVs (and thus codewords) to be 9. So, if the only input blocks starting with 111 are 111100 and 111110, we can use 1111U 0 instead of 111UUU, saving two fill values for each such code. Furthermore, if the test set string includes many input blocks 110100 and 110000, an efficient set of MVs should contain the matching vector 110U 00 (in 9C compression, MVs UUUUUU and UUU000 would be required, respectively). In general, our technique is advantageous when there are input blocks that "almost" match.

Given a test set (that can contain don't-care values), we are looking for a number L of MVs of length K that ensure best compression. K and L are inputs of the algorithm. Since any of the test vectors can have a 0, a 1 or a U on any of its K positions, the space of all possible alternatives is huge. Hence, we employ Evolutionary Algorithms (EAs) [START_REF] Bäck | Evolutionary Algorithms in Theory and Practice[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization & Machine Learning[END_REF] for optimization. EAs have been employed in many areas of VLSI CAD and design automation [START_REF] Mazumder | Genetic Algorithms for VLSI Design, Layout and Test Automation[END_REF][START_REF] Drechsler | Evolutionary Algorithms for VLSI CAD[END_REF], including test pattern generation [START_REF] Girard | A diagnostic atpg for delay faults based on genetic algorithms[END_REF][START_REF] Keim | Combining GAs and Symbolic Methods for High Quality Tests of Sequential Circuits[END_REF] and BIST configuration [START_REF] Ökmen | On optimizing BIST architecture by using OBDD-based approaches and genetic algorithms[END_REF][START_REF] Polian | Evolutionary optimization of Markov sources for pseudo random scan BIST[END_REF], but not in test compression. We assign codewords to the determined MVs using Huffman coding [START_REF] Huffman | A method for the construction of minimumredundancy codes[END_REF]. Experimental results for ISCAS circuits show a substantial improvement compared to 9C compression.

The remainder of the paper is organized as follows: the next section formally defines the fixed-length input block compression with unspecified values. The compression method, including EA-based calculation of matching vectors, is discussed in Section 3. Experimental results are reported in Section 4. Section 5 concludes the paper.

Problem Formulation

Suppose that the circuit has n inputs and the number of test patterns is T . Let tp (i) denote the i th test pattern, and let tp (i) j denote the j th bit position of the i th test pattern (tp

(i) j ∈ {0, 1, X}; 1 ≤ i ≤ T ; 1 ≤ j ≤ n).
Note that X denotes a don't-care value that can be set to 0 or 1 without violating the fault coverage targets. The overall number of bit positions in a test set is T • n. We con-sider the complete test set as one large string of information

t 1 t 2 . . . t T •n := tp (1) 1 tp (1) 2 . . . tp (2) 1 . . . tp (T ) 1 . . . tp (T )
n . Our target is to efficiently encode this string of information.

Let K be a positive integer number. (Note that in contrast to 9C compression it is not required to be odd.) The test set string t 1 . . . t T n is partitioned into fixed-length input blocks of length K, as formalized in the next definition. K should be a divisor of the length T n of the test set string; if that is not the case , the test set string is filled up by adding (K • T n/K -T n) X values in the end.

Definition:

An input block is a subsequence t (j-1)K+1 t (j-1)K+2 . . . t jK of the test set string t 1 . . . t T n (which has the length K), where 1 ≤ j ≤ T n/K (remember that K is a divisor of T n). The k th input block is denoted as

ib (k) ∈ {0, 1, X} K .
2 The fixed-length input block encoding is based on matching vectors to which the input blocks are mapped if they match. Definition: A matching vector (MV) v consists of K bit positions v 1 v 2 . . . v K , where v j ∈ {0, 1, U}, 1 ≤ j ≤ K (where U stands for "unspecified"). An MV v matches an input block ib if there is no bit position j for which either ib j = 1 and v j = 0 or ib j = 0 and v j = 1. (I.e. 1 matches with 1, 0 matches with 0, and X and U match with arbitrary values).

2 Let L be the number of MVs, and let the set of MVs be v (1) , v (2) , . . . , v (L) . Note that as the MVs can contain the U value, different MVs can match the same input block.

Every MV v (i) is assigned a codeword C(v (i) ) ∈ {0, 1} * . We allow the codewords to have different length, but we require the whole code {C(v (1) ), C(v (2) ), . . . , C(v (L) )} to be a prefix code, i.e. no codeword is identical to a beginning of a different codeword. An input block that matches an MV is encoded by the codeword of that MV followed by the values in the input block corresponding to unspecified values (U s) in the MV: Definition: Let a matching vector v have N U (v) unspecified values U at the positions u 1 , . . . , u NU (v) (i.e. 2 Our goal is to efficiently encode a test set string t 1 t 2 . . . t T n (which aggregates the whole test set) divided into T n/K input blocks. We formulate the following subproblems:

v u1 = • • • = v u N U (v) = U ), 0 ≤ N U (v) ≤ K. Let

Problem Formulation

Given the test set (i.e. the sequence of input blocks ib 1 , ib 2 , . . . , ib T n/K ), the size K and the number of MVs L, determine:

2. a prefix code {C(v (1) ), . . . , C(v (L) )}, such that the sum of the encoding lengths of the input blocks is minimal.

2 Note that the problem solved by Tehranipour et al. [START_REF] Tehranipour | Nine-coded compression technique with application to reduced pin-count testing and flexible on-chip decompression[END_REF] is a special case of our generic formulation: they use L = 9, a fixed set of MVs (specified in the introduction) and also a fixed encoding for the MVs.

Solution Approach

We solve the problem formulated above by iteratively applying the following process to the sequence of input blocks corresponding to the test set:

1. Determine L matching vectors (MVs).

(Covering)

Assign an MV to each input block; calculate the frequency-of-use for each MV.

(Encoding)

Encode the input blocks using the collected frequency-of-use data.

Suppose that the L MVs are fixed. Covering is done by selecting an input block and identifying an MV that matches it. If multiple MVs match an input block, we select the one with the lowest number of U s, as the encoding is more compact for such MVs (Note that it is possible that none of the L MVs matches the input block; in this case, encoding is impossible with this set of L MVs. This can be ruled out by setting all positions of one of the MVs to U .) Encoding is done using Huffman's algorithm.

The solution space for determining L MVs of dimension K (i.e. the number of different sets of L K-dimensional vectors over {0, 1, U}) is 3 KL (or 3 KL /(L!), observing that the ordering of the test vectors is irrelevant), which is a huge number even for relatively small values of K and L. Hence, we apply an evolutionary algorithm for this sub-problem. Details are reported in the sections below.

Matching vector determination

We employ Evolutionary Optimization [START_REF] Bäck | Evolutionary Algorithms in Theory and Practice[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization & Machine Learning[END_REF] to determine L matching vectors (MVs). We define an individual as a set of L K-dimensional MVs, i.e. a string v C new individuals (children) are generated by applying evolutionary operators (namely crossover, mutation and inversion) to randomly selected individuals from the initial population. The fitness of each of the C children is determined in the same way as for the individuals in the initial population. From these S + C individuals, S having highest fitness (thus leading to the highest compression) are selected to form the new population (second generation).

(1) 1 v (1) 2 . . . v (1) K v (2) 1 . . . v (L) 1 . . . v (L) K of length KL over the alphabet {0, 1, U}.
The crossover operator takes two individuals (parents) and produces two children by exchanging bit positions (genes) of the parents. Note that for our encoding, an individual has KL genes. The resulting children have genes of one parent in several positions and the genes of the other parent in others. The mutation operator generates one child from one parent by replacing one randomly selected gene of a parent by a random value. The inversion operator produces a child by reverting the ordering of the genes between two random positions of a parent.

The evolutionary optimization process stops when a termination condition is satisfied. We used limits on the number of generated legal solutions and on the number of generations in which no fitness improvement was registered as termination conditions. Then, the MV set given by the fittest individual is returned.

Covering

Once the L matching vectors (MVs) are generated, they are sorted in the order of increasing number of U s:

N U (v (1) ) ≤ N U (v (2) ) ≤ • • • ≤ N U (v (L)
). When matching an input block ib, this sorted list is processed, and the first matching MV v (i) is taken. We will comment on the implications of this technique on the optimality of the overall algorithm in the next section. For every MV v (i) , we calculate its fre-quency F i , i.e. the number of input blocks encoded by v (i) .

Encoding

Given the L MVs v (1) , . . . , v (L) and their frequencies F 1 , . . . F L , the optimal prefix code is obtained by Huffman coding [START_REF] Huffman | A method for the construction of minimumredundancy codes[END_REF]. However, the problem formulation is slightly different here: first, an MV with a frequency of 0 can be simply left out without allocating a codeword to it. Second, due to the existence of the unspecified values (U s) one MV may be subsumed by another. The following example demonstrates a situation in which the Huffman coding algorithm leads to a suboptimal solution due to these specifics: Example: Suppose that the MVs are v (1) = 111U with frequency F 1 = 5, v (2) = 1110 with frequency 3 and v (3) = 0000 with frequency 2. Huffman coding yields C(v (1) ) = 0 , C(v (2) ) = 10 and C(v (3) ) = 11 . Since v (1) contains one U value, it is encoded by | 0 | + 1 = 2 bits; v (2) and v (3) However, 3 input blocks that are matched by v (2) = 1110 are also matched by v (1) = 111U . If we omit v (2) altogether, we can encode v (1) as 0 and v (3) as 1 . Since v (1) now subsumes v (2) , its frequency is now 5+3 = 8. Its encoding is still 2 bits long, while the encoding of v (3) is 1 bit long. The total size of compressed data is 8 • 2 + 2 • 1 = 18 bits, which is less than 20 bits calculated by the Huffman algorithm. 2

The suboptimal result of the example is a consequence of the covering algorithm from the previous section: if the three input blocks had been covered by v (1) rather than v (2) , the Huffman algorithm would have produced the optimal result. Handling such cases explicitly could improve the compression rate.

Experimental Results

We applied our method to uncompacted stuck-at test sets with don't-cares obtained by the method from [START_REF] Kajihara | On identifying don't care inputs of test patterns for combinational circuits[END_REF] for IS-CAS 85 and combinational parts of ISCAS 89 circuits, and to uncompacted path delay test sets with don't-cares for combinational parts of ISCAS-89 circuits generated by the tool TIP [START_REF] Henftling | Bit Parallel Test Pattern Generation for Path Delay Faults[END_REF][START_REF] Tafertshofer | A SAT-based implication engine for efficient ATPG, equivalence checking, and optimization of netlists[END_REF] (with 100% robust path delay coverage). We used the package GAME [START_REF] Göckel | GAME: A software environment for using genetic algorithms in circuit design[END_REF] for evolutionary optimization.

Table 1 summarizes the results for stuck-at test sets. The first two columns contain the name of the circuit and the size of the test set (note that the circuits are sorted with respect to increasing test set size). All other columns quote the compression rate 100% Column '9C' reports the compression rate obtained by our reimplementation of the 9C compression method from [START_REF] Tehranipour | Nine-coded compression technique with application to reduced pin-count testing and flexible on-chip decompression[END_REF] for K = 8 (which yielded best results in that paper). Note that our numbers do not match the figures in that paper, as we use different test sets. Since the compression rates are higher for the test sets that we use than the rates reported in [START_REF] Tehranipour | Nine-coded compression technique with application to reduced pin-count testing and flexible on-chip decompression[END_REF], it would be unfair to quote those results. We did not

  an input block ib be matched by v. The encoding C(ib, v) of ib by v is given by C(v)ib u1 ib u2 . . . ib u N U (v) . The length of the encoding |C(ib, v)| is the sum of the length of C(v) and the number N U (v) of the unspecified values in v. Note that |C(ib, v)| is independent from ib.

Figure 1 Figure 1 :

 11 Figure 1: Evolutionary Algorithm

  are encoded by | 10 | = 2 bits and | 11 | = 2 bits, respectively, resulting in the total size of compressed data of 5 • 2 + 3 • 2 + 2 • 2 = 20 bits.

Table 1 :

 1 • original test set sizesize of compressed data original test set size of various methods. The compression is better for higher values of this number. Experimental results for stuck-at test sets

	Circuit	Test set	Compression rate
		size	9C 9C+HC	EA EA-Best
	s349	624 23.0% 30.0% 54.2%	55.8%
	s344	624 25.0% 33.0% 51.8%	55.8%
	s298	629 19.0% 27.0% 45.2%	51.2%
	s208	722 26.0% 32.0% 47.8%	50.4%
	s400	984 29.0% 36.0% 54.4%	56.4%
	s382	1008 29.0% 36.0% 52.0%	54.2%
	s386	1157	0.0% 13.0% 30.4%	30.6%
	s444	1176 40.0% 43.0% 54.4%	57.8%
	c6288	1216	8.0% 19.0% 17.6%	20.4%
	s510	1850 42.0% 45.0% 57.6%	57.6%
	c432	1944 26.0% 36.0% 49.2%	50.4%
	s526	1944 25.0% 29.0% 46.4%	46.4%
	s1494	2324 -1.0% 11.0% 23.0%	28.9%
	s420	2380 53.0% 55.0% 54.4%	56.2%
	s1488	2436	2.0% 15.0% 25.6%	30.0%
	s832	3404 35.0% 38.0% 43.8%	43.8%
	s820	3496 31.0% 35.0% 42.8%	43.4%
	c499	3854 43.0% 51.0% 45.0%	51.6%
	s713	4104 51.0% 52.0% 61.4%	61.8%
	s641	4212 51.0% 52.0% 60.2%	62.2%
	c880	4680 40.0% 42.0% 47.8%	49.8%
	c1908	4950 -2.0% 10.0% 18.4%	19.0%
	s953	5220 51.0% 53.0% 61.6%	63.2%
	c1355	5289 38.0% 45.0% 40.8%	44.8%
	s1196	6016 34.0% 38.0% 46.2%	46.2%
	s1238	6240 34.0% 37.0% 44.0%	45.8%
	s1423	8463 59.0% 59.0% 61.0%	61.6%
	s838	8509 67.0% 68.0% 66.2%	68.6%
	c3540	10350 36.0% 39.0% 43.8%	44.2%
	c2670	33086 70.0% 70.0% 70.4%	70.6%
	c5315	33108 65.0% 65.0% 66.2%	67.0%
	c7552	60030 63.0% 64.0% 63.2%	63.2%
	s5378	71262 73.0% 73.0% 76.8%	76.8%
	s9234	118560 75.0% 75.0% 76.2%	76.4%
	s35932	133988 71.0% 71.0% 73.8%	73.8%
	s15850	305500 80.0% 80.0% 83.0%	83.0%
	s13207	410200 83.0% 83.0% 85.8%	85.9%
	s38584	1250256 82.0% 82.0% 86.2%	86.2%
	s38417	2068352 84.0% 84.0% 87.0%	87.9%
	Average		42.6% 46.8% 54.2%	55.9%

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Acknowledgment

We would like to thank Kohei Miyase and Seiji Kajihara for providing the stuck-at test sets for ISCAS circuits containing don't-care values.

 [START_REF] Tehranipour | Nine-coded compression technique with application to reduced pin-count testing and flexible on-chip decompression[END_REF]: '0' for 000000; 10 for 111111; 11000 for 000111; '11001' for 111000; '11010' for 111UUU; '11011' for UUU111; '11100' for 000UUU; '11101' for UUU000; and '1111' for UUUUUU (for K = 6; the same assignment applies for other K values). We ran the reimplemented 9C algorithm with this encoding replaced by Huffman coding and report the results in column '9C+HC'.

The results of our proposed approach are given in column 'EA'. We used L = 64 and K = 12; in the evolutionary algorithm, we used the population size S of 10, the number of children C of 5, the crossover probability of 30%, the mutation probability of 30% and the inversion probability of 10%. One of the MVs was set to all-U , such that there were no insolvable instances. We report the average over 5 algorithm runs.

We generated data for numerous values of K and L. While the space does not allow us to include all data, we report our best results in the last column. Additional improvements could be achieved by varying the parameters of the evolutionary algorithm, which we did not do.

Our proposed approach outperforms the original 9C algorithm for all circuits except s838 (but our best result is higher for this circuit). This could be ruled out by adding the 9C matching vector set to the initial population (which we did not). If 9C compression is combined with Huffman coding, it outperforms our method with default values for 6 circuits, and it is more efficient than our best results for two circuits (out of 39). This can be explained by the large size of the solution space; a larger population size or a reduced evolutionary pressure might be appropriate for such circuits.

The last line of the table quotes the average values. It can be seen that our results are considerably better on average than those obtained by 9C compression, both with and without Huffman coding. It can also be seen that the difference between our technique with default values and the best compression rate is relatively small. We conclude that our algorithm is relatively stable with respect to the parameters K and L.

The structure of Table 2 is similar to that of Table 1. We report the results for the proposed method for L = 9 and K = 8 (column 'EA1') and L = 64 and K = 12 (column 'EA2'). Both values are averages over 5 successful algorithm runs with termination after 500 populations without improvement. The All-U MV was included into the population. The proposed methods yields best results for all but two circuits (s420 and s838), and the second set of parameters results in an improvement of the compression from 48.7% for the original 9C compression to 58.6%, on average.

Conclusions

We provided an Evolutionary Optimization-based solution for test set encoding with fixed-length input blocks. For the first time, our problem formulation allows unspecified values (U s) at arbitrary positions of the matching vectors. This enables the employment of compact on-chip decoders for arbitrary test sets. No decoder re-design is required in case of a test set modification, if an all-U matching vector is used; however, the compression rate might suffer. A reconfigurable decoder, into which the codeword/matching vector information can be loaded, would solve this problem.

Experimental results show the superiority of our technique. However, further improvements are possible by fitting the parameters of the Evolutionary Optimization, such as population size and operator probabilities. Another direction for further research is the application of our method in a multiple scan chain environment.