Askan Can

Mavi
email: camavi@univ-paris1.fr

Antoine D'autume

Katheline Schubert

Katrin Hippolyte D'albis

Mouez Millock

Bertrand Fodha

Fabien Wigniolle

Lionel Prieur

Karine Ragot

Constant

What Can Abrupt Events Tell Us About Sustainability ?

HAL is

Introduction

The fact that uncertain abrupt climate events could cause large scale damages is widely recognized [START_REF] Alley | Abrupt climate change[END_REF], [START_REF] Field | Managing the risks of extreme events and disasters to advance climate change adaptation[END_REF]). A very considerable number of studies focus on decision making concerning exploitation policy of natural resources under uncertainty [START_REF] Bretschger | Best policy response to environmental shocks : Applying a stochastic framework[END_REF], [START_REF] Tsur | Pollution control in an uncertain environment[END_REF], 2007a, 2016a), [START_REF] Hr Clarke | Consumption/pollution tradeoffs in an environment vulnerable to pollutionrelated catastrophic collapse[END_REF]). Some recent literature examines also the ways to deal with the uncertainty. For this purpose, adaptation and mitigation policies and their implications under uncertainty is one of the major point of interest [START_REF] Zemel | Adaptation, mitigation and risk: An analytic approach[END_REF], [START_REF] Tsur | Policy tradeoffs under risk of abrupt climate change[END_REF], [START_REF] Mavi | Can a hazardous event be another source of poverty traps? FAERE Working Paper[END_REF].

Other than the studies on uncertainty and resource exploitation, another branch of the literature concentrates on the relationship between discounting and sustainability which is one of the lasting important debate in economics literature. Especially, the debate has been intensified in the context of climate change [START_REF] Stern | Report on the economics of climate change[END_REF], [START_REF] Weitzman | A Review of The Stern Review on the Economics of Climate Change[END_REF], [START_REF] Heal | Climate economics: a meta-review and some suggestions for future research[END_REF]). Some of the studies within this debate concern the role of individual time preferences (see [START_REF] Endress | Intergenerational equity with individual impatience in a model of optimal and sustainable growth[END_REF], [START_REF] Schneider | Trading off generations: Equity, discounting, and climate change[END_REF], Marini andScaramozzino (1995, 2008), [START_REF] Burton | Intertemporal Preferences and Intergenerational Equity Considerations in Optimal Resource Harvesting[END_REF]). The presence of individual time preferences in an economic model is appealing because infinitely lived agent model (ILA hereafter) is criticized for not respecting consumer sovereignty 1 . 1 The articles cited above which are incorporating individual discount rate are based on the famous framework proposed by [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF]. Authors are introducing individual time preferences (i.e individual discount rate) in an OLG model. Then, they find statically aggregate consumption level of all generations at a given time. Once aggregation is made, the model reduces to representative agent framework. This framework is well used in environmental economics to treat various problem as inter-generational equity2 by the above cited papers. In these models, there is any effect of individual discount rate on aggregate long term dynamics, neither on transitional dynamics. One may easily say that model introduces a clear dichotomy between an OLG and ILA model since individual discount rate has nothing to tell about the aggregate dynamics of the economy.

However, one one hand, the studies treating long term impacts of uncertainty on resource exploitation policy do not take into account sustainability and intergenerational equity. On the other hand, the strand of the literature on sustainability and intergenerational equity does not take into account uncertain events.

In this sense, to the best of our knowledge, we can argue that the link between sustainability and abrupt events is overlooked in environmental economics literature.

In this paper, our aim is to make a bridge between sustainability and uncertain abrupt events. To do this, we take a step further to answer the following central questions : What are the implications of abrupt events regarding the sustainability of an economy ? What can individual time preferences tell about sustainability in an economy subject to abrupt events ?

The contribution of this paper is threefold : First, by using the [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF], we show that an economy exposed to an abrupt event could have limit cycles (Hopf bifurcation) at the long run. On the one hand, the economy accumulates physical capital and creates wastes. In this sense, the environment is used as a "sink". This can be considered as an economic goal. On the other hand, since the inflicted damage after the abrupt event is proportional to remaining natural resource stock, the economy wishes to protect the natural resource. This is the environmental goal. When it becomes difficult to steer between two conflicting policies, it may be optimal to cycle around the optimal steady state3 .

One may argue that this trade-off between the economic and environmental goal is usual in growth models with environment. For this purpose, we show rigorously that with waste stemming from physical capital accumulation and the abrupt event probability, the preferences of the economy over time are not independent in the sense of Koopmans (see [START_REF] Tjalling | Stationary Ordinal Utility and Impatience Author ( s ): Tjalling C . Koopmans[END_REF]). Then, when the economy faces the above mentioned trade-off coupled with complementarity over time, it is likely that the economy is exposed to limit cycles (see [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF], [START_REF] Heal | Optimal Growth with Intertemporally Dependent Preferences[END_REF]).

Indeed, limit cycles are examined extensively in environmental economics. [START_REF] Wirl | Complex, dynamic environmental policies[END_REF][START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF] and Bosi andDesmarchelier (2016a, 2017) study the existence of limit cycles in models with representative agent framework. Nonetheless, none of these studies link limit cycles to sustainability.

At this point, the question to be addressed is : what are the implications of limit cycles regarding sustainability ? The Sustainable Development Criterion requires that the utility of consumption should have a non-decreasing path (i.e du(c(t)) dt ≥ 0). If the economy is exposed to limit cycles due to the trade-off between the environmental and the economic goal, consumption increases at one date and decreases at another date.

Then, the Sustainable Development Criterion is not respected since the utility has a cyclical behavior at the long run.

Secondly, contrary to the [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF] framework and to the articles using this framework, we show that individual time preferences could change the stability properties of the model. This result disproves the conventional result which says that aggregate dynamics are solely governed by the social planner's discount rate (see [START_REF] Endress | Intergenerational equity with individual impatience in a model of optimal and sustainable growth[END_REF], [START_REF] Schneider | Trading off generations: Equity, discounting, and climate change[END_REF], Marini andScaramozzino (1995, 2008), [START_REF] Burton | Intertemporal Preferences and Intergenerational Equity Considerations in Optimal Resource Harvesting[END_REF]). Indeed, we show that the individual discount rate has an important role concerning the sustainability of an economy.

Since the first part of the model is an OLG model, there is an intra-generational allocation of consumption which is stable over time. We show that intra-generational equity comes as a cost since it is more likely to have a stable equilibrium when social planner allocates unequally the consumption over different cohorts of age. It is shown that when the consumption is allocated equally between generations, the economy faces limit cycles at the long run which compromises the sustainability of the economy.

The reason why individual discount rate can raise limit cycles is the following : when the social planner aggregates the consumption over all agents, she finds an aggregate utility function and post value function which are dependent on social and individual discount rate. Different levels of individual discount rate changes the weight of utility function and post value function in the objective function of the social planner. This is equivalent to say that the trade-off between the economic and the environmental goal changes with respect to different individual discount rate.

It is argued that sustainability and intergenerational equity are generally perceived as normative inquiries [START_REF] Solow | Sustainability, an economist perspective[END_REF][START_REF] Solow | An almost practical step toward sustainability[END_REF]). Then, a social planner who pays attention to the sustainability and intergenerational equity should seek to avoid limit cycles. We show that the social planner could avoid limit cycles by protecting the environment. This is due to the fact that a higher natural resource stock implies a lower marginal utility of consumption. As a result, different levels of individual discount rate could not change so much the tradeoff between the economic and the environmental goal. Consequently, there are less chances that an economy traps to limit cycles at the long run.

The remainder of the paper is organized as follows. The section 2 presents the benchmark model and explains all economic mechanisms behind limit cycles in detail. Section 3 explains the model with environmental protection and focuses from a normative point of view in order to avoid limit cycle and last section concludes the paper.

Model

Similar to [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF], the model is constituted by two parts. In the first part, the social planner allocates consumption across different cohorts of age in an overlapping generations model (OLG).

At each time t, the economy is composed by individuals ranged between age 0 and N and also no two individuals have the same age. In the second part, the model is solved first by aggregating the consumption over all agents. At this point, the model reduces to an infinitely-lived agent model and one can maximize intertemporal utility and find optimal aggregate trajectories of the economy.

The lifetime utility of an individual born at date T (or equivalently, the aggregate utility of all living generations at a given date t) is written

U (c (t)) = N 0 u (c (t, τ )) e -βτ dτ (1)
where β is the individual discount. c (t, τ ) is the consumption level of an agent with age τ at time t. The utility of the society is a weighted sum of discounted utilities of individuals of the society augmented by the post-event value.

W C b = Eb b -N N 0 u (c (t, τ )) e -(β-ρ)τ dτ e -ρb db + e -ρ bϕ S b (2)
The social planner weights generations by her discount rate ρ. An abrupt event occurs at an uncertain date b at which a generation born. This generation can be considered as an "unlucky" generation. The variable S represents the natural capital stock. The expression ϕ (S) stands for the post-value function that describes the value of the problem after a catastrophic event. The constant parameter θ represents Poisson probability of catastrophic event. In our model, similar to [START_REF] Bretschger | Best policy response to environmental shocks : Applying a stochastic framework[END_REF], the use of constant hazardous event probability implies that the environmental shock is to a flow variable rather than a stock variable.

In the following model, note that natural capital accumulation is negatively affected by physical capital accumulation, which is similar to [START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF] and [START_REF] Ayong | Sustainable growth, renewable resources and pollution[END_REF]. In this sense, nature is considered as a "sink" for wastes coming from physical capital. Essentially, this feature creates a trade-off between capital accumulation and its negative effects to nature that we will analyze in depth within this study. The social planner seeks to choose the optimal level of consumption of agents c (t, τ ) for each time t subject to two physical constraints, i.e.,

             max c(t,τ ) W C b K (t) = f (K (t)) -δK (t) - N 0 c (t, τ ) dτ Ṡ (t) = G (S (t)) -γf (K (t)) (3)
The maximization program can be solved with optimal control techniques. Since the state variables K and S evolves in terms of pure time t, we can solve the problem of the social planner by taking the expectations of the objective function and reformulating it in terms of time t, instead of born date b of generations (see Appendix for calculations). For this purpose, we switch from (b, τ ) to (t, τ ) by using b = t -τ . We can reformulate the social planner's program in the light of these information.

max c(t) W = ∞ 0 {U (c (t)) + θϕ (S (t))} e -(ρ+θ)t dt (4) subject to    K (t) = f (K (t)) -δK (t) -c (t) Ṡ (t) = G (S (t)) -γf (K (t)) (5)
where the objective function consists of the allocation of consumption among different groups of age by the social planner.

U (c (t)) + θϕ (S (t)) = max N 0 u (c (t, τ )) + θϕ S b e -(β-(ρ+θ))τ dτ
The utility function of different cohorts is maximized statically with respect to constraint

N 0 c (t, τ ) dτ ≤ c (t) (6)
Once the social planner allocates the consumption between different cohorts of agent, she aggregates the consumption over all agents at a given date t. Otherwise, with the use of an endogenous hazardous event, since the utility becomes non-stationary over time (see [START_REF] Tjalling | Stationary Ordinal Utility and Impatience Author ( s ): Tjalling C . Koopmans[END_REF]), the social planner's decisions will be time inconsistent.

First stage optimization : A static problem

In the first stage, the social planner finds a relationship between c (t, 0) and c (t, τ ). Indeed, the social planner finds the optimal consumption allocation between different agents. This stage of the problem consists on a static optimization in an OLG model where c (t) is the aggregate consumption level at time t. We write the static optimization problem

L = N 0 u (c (t, τ )) e -(β-(ρ+θ))τ + θϕ (S (b)) dτ + λ (t) c (t) - N 0 c (t, τ ) dτ
For the sake of simplicity, after the abrupt event, the social planner does not maximize the utility since the consumption level is assumed to reduce to a subsistence level c min (see [START_REF] Bommier | Models-as-usual for unusual risks? On the value of catastrophic climate change[END_REF], Tsur and Zemel (2016a)). Then the post-event value does not depend on c (t, τ ). From this problem, we find easily

u c (c (t, τ )) = u c (c (t, 0)) e (β-(ρ+θ))τ (7) u (c (t, τ )) = c (t, τ ) 1-σ -c 1-σ min 1 -σ (8)
The maximization of the static program yields

c (t, τ ) = c (t, 0) e -(β-(ρ+σ))τ σ
The social planner gives more weight to younger generations when the individual discount rate is higher than her discount rate. This makes sense when one thinks about the sum of discounted utility of young and old agents. Since the utility of young agents are discounted over a longer time horizon, the social planner gives much more weight to a young agent than an old one. In addition, the social planner allocates equally consumption across agents when both individual and social discount rate coincide. Note that social planner tends to treat generations more equally when individual discount rate is closer to her own discount rate.

Graphically, it is easy to see how a social planner could allocate the consumption between different cohorts at some time t. The analysis of the effect of abrupt event probability on intragenerational allocation is also appealing. Regarding the explanation above, one can observe that a higher abrupt event probability shifts the consumption allocation from young to older cohorts.

Figure 3: Effect of risk on intragenerational allocation

This is due to the fact that a higher event probability, adjusted discount rate of social planner ρ + θ increases. Then, older agent's utility starts to take more weight with a higher abrupt event probability, which shifts the consumption from younger agents to older cohorts. We find aggregate utility function (see Appendix for calculations).

Once the social planner allocates consumption across agents, we need to find the utility of aggregate level of consumption and post-event value. With the CRRA form utility function, we have (see Appendix for details.)

U (c (t)) = c (t) 1-σ -c 1-σ min 1 -σ (Z 1 ) σ (9)
where Z 1 is the aggregation term which includes individual discount rate. The post value function ϕ (S) = u (c min ) -ψ (S) can be written in a similar way to Tsur andZemel (1998, 2016a). After the aggregation of post-event value, we have

ϕ (S) = Z 2 (u (c min ) -ψ (S)) (10)
where Z 2 is the aggregation term of the post-value function. We take the simplest possible form of post-value function in line with Tsur and Zemel (2016a). (see Appendix for the derivation concerning the aggregation of different cohorts of age.) The single occurrence event implies that consumption is reduced to a minimum level4 which does not provide utility. The penalty rate due to abrupt event is proportional to the level of natural capital stock after catastrophe.

ψ (S) = ψ (ω 1 -ω 2 logS) (11)
where ω 1 is the unrecoverable part of penalty. We put the condition that ω 1 > ω 2 logS to ensure that some part of the damage to the environment is permanent. The use of the penalty function implies that more the natural capital stock is protected, the less the economy suffers from penalty.

Second stage optimization : A dynamic problem

The social planner maximizes the program (4) subject to physical and natural capital accumulation constraints (5). The current-value Hamiltonian for maximizing W

H = U (c) + θϕ (S) + λ (f (K) -δK -c) + µ (G (S) -γf (K)) (12)
The first order conditions and dynamics of the economy are as follows

                     U c = Z -σ 1 c -σ = λ K = f (K) -δK -Z 1 λ -1 σ Ṡ = G (S) -γf (K) λ = (ρ + θ) λ -(f K -δ) λ + µγf K μ = (ρ + θ) µ -µG S -θϕ S (13)
By using the functional forms above, the steady-state of the economy can be written as a function of renewable resource stock S :

         K * (S) = (1-S * )S * γA λ * (S) = γAµ * (S) ((A-δ)-(ρ+θ)) µ * (S) = Z2 ψω2θ S * ((ρ+θ)-(1-2S * )) (14)

Proposition 1. In an economy with catastrophic event probability, a Hopf bifurcation is possible for different critical parameters of individual discount rate β.

It is important to figure out the economic reasons behind the occurrence of limit cycles. On one hand, more natural capital stock is high, less the penalty rate would be. (see [START_REF] Tsur | Pollution control in an uncertain environment[END_REF] This can be described as an environmental goal. On the other hand, the resource stock in this model is used as a sink for waste coming from physical capital accumulation5 . This represents the economic goal. Trying to decide between these two opposite strategies could lead to limit cycles around steady state when it becomes torn to steer between two goals. (see [START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF] One may argue that this trade-off is usual in growth models with environment and does not cause instability. Indeed, the limit cycles and unstable behavior in this economic model is due to conflicting economic and environmental goals coupled with complementarity of preferences over time (see [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF], [START_REF] Heal | Optimal Growth with Intertemporally Dependent Preferences[END_REF]). An incremental increase of consumption at time t 1 shifts the preferences from t 3 to t 2 (i.e. adjacent complementarity and distant complementarity vice versa.) if there is complementarity of preferences over time. If not, an incremental change of consumption at time t 1 does not have any effect on the preferences of other dates. This implies that preferences are intertemporally independent.

The model is shown to have complementarity of preferences over time (see Appendix). To understand the implications of complementarity, consider an increase of consumption near date t 1 which shifts a part of consumption of date t 3 to t 2 . At time t 1 , physical capital accumulation decreases. It follows that the economy accumulates less waste and protects the natural stock. On the contrary, since at date t 3 , consumption decreases and the economy accumulates more physical capital. Then, there is more waste accumulation which harms the environment and so on. From this mechanism, we can understand the intuition behind optimal limit cycles. Due to the complementarity over time, natural resource stock increases and decreases at consecutive dates, from which a cyclical behavior occurs.

Note that this mechanism does not hold if there is not waste stemming from physical capital accumulation. We also show that when waste rate γ approaches zero, the complementarity over time vanishes. (see Appendix.) We can also remark that the dynamic system (13) reduces to a block recursive system of (K, λ) and (S, µ) which admits a saddle path equilibrium when there is not waste accumulation in the economy.

The legitimate question to be addressed is why and how individual preferences have a role in explaining limit cycles and unstable behavior in this economy. It is evident that the individual discount rate β is crucial since it changes the difficulty of steering between economic and environmental goal. This can be remarked through the utility of aggregate consumption which depends on the individual discount factor. More importantly, not every level of individual discount rate causes limit cycles or unstable spirals but a special set of parameters causes. Despite the fact that we try simple functional forms, it is always a hard task to find analytically the critical value for bifurcation parameter. (see [START_REF] Wirl | Complex, dynamic environmental policies[END_REF][START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF])) As in [START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF], we refer to a numerical analysis to study in which cases bifurcation is a possible outcome. The steady state level is stable for low and high levels of individual discount rate β 7 . At β = 0.925644 and β = 0.974355, we have a pair of purely imaginary eigenvalues where Hopf bifurcation occurs. When we have β ∈ [0.925644.., 0.974355..], the economy is exposed to unstable spirals. As mentioned above, it is possible to understand the intuition behind why limit cycles and unstable spirals take place but we are unable to offer an economic explanation why Hopf bifurcation takes place at given two different critical values 8 .

It is clear that the economy becomes more aggressive to exploit natural resources when the individual discount rate is higher until the threshold after which individual discount rate exceeds the augmented discount rate of social planner. (see Figure ( 4) ) After this threshold, since the economy starts to be more precautionary since the social planner is more patient than individuals.

In the figure (4), we observe that instability and limit cycles occur in a region where individual discount rate decreases and increases the optimum steady state level of the natural resource stock. This is equivalent to say that a higher individual discount rate does not imply necessarily a lower stock of natural resources.

The main reason why the optimum steady state has a convex form with respect to individual discount rate relies also on the trade-off between environmental and economic goal. In order to understand why optimum steady state level of natural resource stock S with respect to different individual discount level has a convex form, we can reformulate the steady state level of consumption in the following way

c * (S) = Z [S * ((ρ + θ) -(1 -2S * )) ((A -δ) -(ρ + θ))] γA ψω 2 θ (A -δ) ( 15
)
where Z = Z1 (Z2) σ . With a higher individual discount rate β, the aggregation term for utility and post value function Z 1 and Z 2 are shown to decrease with respect to individual discount rate β (see Appendix).

The first effect is that in the case where Z decreases, consumption decreases, physical capital stock tends to increase and waste increases also. Then, the natural resource stock decreases. We have the reverse mechanism when the term Z increases.

According to the figure ((4)), the instability occurs in an area where it is hard to decide for the management of natural resources. When the individual discount rate is just below the social planner's adjusted discount rate, a higher individual discount rate makes decrease the natural resource stock. Contrary to this, when individual discount rate is above the social planner's discount rate, the natural resource stock starts to be conserved. One may understand from this explanation that it is difficult to steer between environmental and economic goal when the marginal utility of consumption is very sensible to small variations of individual discount rate β. We can also say that agents are in a sort of dilemma to consume more or less when their discount rate is closer to social planner's adjusted discount rate. Consequently, in this case, it becomes optimal to cycle around the steady state.

Definition 1. A path of utility respects the Sustainable Development criterion if dU

(t) /dt ≥ 0.
The Sustainable Development criterion states that utility of consumption should follow a non-decreasing or at least a constant path in order to ensure the sustainability of the economy. This means that limit cycles for consumption violate sustainability criterion since consumption decreases at some moments of time t.

An important point is that the model reduces to a representative agent model at the second stage.

Consequently, we can not talk about the intergenerational equity as in [START_REF] Schumacher | Pollution perception: A challenge for intergenerational equity[END_REF] for δ = 0.001, γ = 0.0001, ψ = 1,θ = 0.1 7 The Lyapunov number for both critical points are negative. This means that model shows a sub-critical Hopf bifurcation at two different critical individual discount rate β. For β = 0.925644 and β = 0.974355, the Lyanpunov numbers are -0.011645276 and -0.0116452956 respectively.

8 See Wirl (1991[START_REF] Wirl | CYCLICAL STRATEGIES IN TWO-DIMENSIONAL OPTIMAL CONTROL MODELS: NEC-ESSARY CONDITIONS AND EXISTENCE[END_REF] for a detailed discussion regarding the difficulty of explaining economically bifurcation points.

aggregate dynamics because we don't know the length of life cycle of a generation along time t. For example, each generation can "equally" face a decrease or an increase of the consumption. Then, the equity between generations is respected. For this reason, we are focusing on the link between limit cycles and sustainability notion in our analysis.

From a normative point of view, the social planner should think about the way of avoiding cycles. We will address this issue in the following section.

Proposition 2. When social planner treats all generations equally (i.e β = ρ+θ), sustainability criterion is not respected9 .

Another interesting result is that the way how a social planner allocates the consumption over different generations has an impact on long term dynamics of the economy (i.e sustainability). Treating generations equally comes at a cost since limit cycles and instability occur when individual discount rate is closer to social planner's adjusted discount rate (ρ + θ). When social discount rate is close to individual discount rate, the social planner tends to treat generations equally but she is also likely to trap the economy to limit cycles which compromises sustainability. . As a result, a new trade-off between intragenerational equity and sustainability arises in this framework when the economy faces a difficulty of steering between economic and environmental goal.

The economic intuition behind limit cycles can be understood from a mathematical point of view. (see Appendix for proof.) Note that limit cycles occur mostly when there exists complementarity over time. [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF]) This is equivalent to say that preferences of economic agents are interdependent, different than the growth models where the preferences of agents over time are independent in the sense of Koopmans (see [START_REF] Tjalling | Stationary Ordinal Utility and Impatience Author ( s ): Tjalling C . Koopmans[END_REF]). [START_REF] Wirl | CYCLICAL STRATEGIES IN TWO-DIMENSIONAL OPTIMAL CONTROL MODELS: NEC-ESSARY CONDITIONS AND EXISTENCE[END_REF] explains that limit cycles take place when synergistic between control and state and/or between states is strong. In the model, it is evident that waste accumulation creates a strong link between two state variables. We can also verify the claim of [START_REF] Wirl | CYCLICAL STRATEGIES IN TWO-DIMENSIONAL OPTIMAL CONTROL MODELS: NEC-ESSARY CONDITIONS AND EXISTENCE[END_REF] by looking at complementarity over time. It is shown that complementarity over time vanishes when there does not exist any waste coming from physical capital accumulation. (see Appendix.)

Proposition 3. (a) In an economy without abrupt event, individual discount rate has no effect on long term dynamics.

(b) The economy admits a saddle path equilibrium when there is any abrupt event probability.

Proof. See Appendix.

This result is quite plausible when one thinks about the absence of the trade-off between decreasing effect of natural stock on penalty rate (environmental goal) and the use of it as a sink for waste (economic goal).

Then, the abrupt event probability plays a crucial role since it changes the marginal utility which changes also the above mentioned trade-off between capital accumulation and the environment.

Intuitively, one may say that a very high catastrophic event may push the society to be more precautionary about environment since its importance is high in order to have a lower penalty rate. Equivalently, this means there would be less physical capital accumulation in this scenario.

On the contrary, when there is relatively very low catastrophe event probability, society becomes indifferent to environmental concerns and use the environment as a sink without caring much about it. In a nutshell, when social planner faces a catastrophic event probability between "too high" and "too low" probabilities, she can not decide firmly between these two strategies. Consequently, it becomes optimal to cycle arount the steady state. We also show in appendix that when the economy faces an abrupt event probability very high or closer to zero, the complementarity of preferences over time vanishes.

Another remark could be that individual discount rate β has no effect on stability of the equilibrium if there is not a catastrophic event probability. As t which extent this trade-off is important. mentioned above, the steady state is stable for relatively high and low levels of catastrophe probability. A closer intuitive explanation to the mechanism explained for limit cycles in this study is proposed by [START_REF] Heal | The use of common property resources[END_REF] and Bosi and Desmarchelier (2016a)which are based on compensation effect. Assume that economy is at steady state at a given date t and assume that natural capital stock S11 decreases exogenously.

The degradation of natural capital pushes agents to increase their consumption since they would like to compensate the disutility due to the decrease of natural capital stock. It follows that capital accumulation decreases as well, which increases the natural capital as there is less waste and so on : deterministic cycles arises. Since in our specification, the objective function is separable in consumption and natural resource stock, we are unable to justify the limit cycles by a compensation effect. However, we show that even with a non-additive objective function, we show that limit cycles exist for a set of parameters. 6) shows that at the first bifurcation point, the utility is exposed to cycles and does not converge to a stable point. From a normative point of view, a social planner that pays attention to sustainability should avoid any path that leads to bifurcations. The next section analyzes the necessary policy that avoids limit cycles.

Model with abatement activities

In this section, our focus will be on the implications of abatement activities on limit cycles. The social planner's program augmented by abatement activities is

W = ∞ 0 {U (c (t)) + θϕ (S (t))} e -(ρ+θ)t dt (16)
The dynamics of capital accumulation contains the cost of mitigation M .

   K (t) = f (K (t)) -δK (t) -c (t) -M (t) Ṡ (t) = G (S (t)) + Γ (M (t)) -γf (K (t)) ( 17
)
where

Γ (M ) = M α , α > 0
holds for abatement activities such as reforestation, desalination of water stock, enhancing carbon sinks etc. The specification of abatement is in the same line with [START_REF] Ariaster | Total factor productivity and the environmental Kuznets curve[END_REF]. Alternatively, function Γ (M ) can be thought as "environmental protection function". The expenditures for environmental protection may be directed not only toward pollution mitigation but also toward protection of forests, recovery of degraded areas. Equivalently, abatement activity in this model can be seen as a means of improving the environmental quality.

The first order conditions and dynamics of the economy augmented by abatement activity are as follows

                     U c = Z 1 λ -1 σ = λ K = f (K) -δK -Z 1 λ -1 σ -M Ṡ = G (S) + Γ (M ) -γf (K) λ = (ρ + θ) λ -(f K -δ) λ + µγf K μ = (ρ + θ) µ -µG S -θϕ S (18)
The steady-state of the economy as a function of natural resource stock is

                   K * (S) = (1-S * )S * + α µ * (S * ) λ * (S * ) α 1-α γA M * (S) = α µ(S * ) λ(S * ) 1 1-α λ * (S) = γAµ * (S * ) ((A-δ)-(ρ+θ)) µ * (S) = ψω2θ S * ((ρ+θ)-(1-2S * )) (19)
Proposition 5. Abatement activity makes limit cycles less likely to occur by increasing the determinant and decreasing the of sub-matrices of Jacobian, in which case the economy admits a saddle-path stable equilibrium.

Proof. See Appendix.

We also show numerically that model augmented by abatement activities admits always a saddle path equilibrium. whatever the level of individual discount rate is, the economy always admits a saddle path equilibrium. The economic explanation is the following : when the steady state level of natural resource stock is higher due to environmental protection activity, the marginal utility of consumption (see ( 19)) is lower. In this case, the variations of consumption due to different levels of individual discount rate β would have a lower impact on marginal utility. Indeed, when the economy is trapped to limit cycles due to a tight trade-off between economic and environmental goal, the abatement activity relaxes this tight trade-off by lowering the marginal utility. Then, the economy escapes the dilemma between environmental and economic goal by giving more weight to environmental goal.

Conclusion

In this paper, we showed the existence of limit cycles in an economy exposed to abrupt event probability.

The limit cycles are caused by conflicting economic and environmental goals coupled with complementarity of preferences over time. An interesting finding of the paper is that individual time preferences of agents other than the social planner's discount rate is crucial not only for intra-generational equity but also for the sustainability of an economy. It is shown that when individual discount rate is close to social planner's discount rate, intra-generational equity is ensured. However, the equality between living generations comes at a cost since it is more likely that the economy is trapped at limit cycles at the long run. This result also disproves a wide spread result in the literature which says that aggregate dynamics are solely governed by the social planner's discount rate. Since we adopt a normative point of view, the limit cycles are considered an undesirable result since it compromises the sustainability according to the prominent Sustainable Development Criterion. Then, we seek to avoid these cycles by protecting the environment. The result in the model augmented by environmental protection activities shows that it is less likely to have limit cycles at the long run with an environmental policy aiming at improving/protecting the environment.

Appendix

Aggregate Economy facing a catastrophic event

In order to understand clearly the problem, firstly we present the economy without a catastrophic event from the starting point as in [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF];

W N C = ∞ -N N 0 u (c (t, τ )) e -βτ dτ e -ρb db ( 20
)
where b is the born date of a cohort and script N C indicates the welfare in a no-catastrophe situation. A social planner, once she makes intra-generational allocation, decides for an optimal allocation of consumption between generations. In this form of the social welfare function W , time does not enter directly in the utility function of an agent. Since born date b is equal to t -τ where τ stands for age, we can reformulate the social planner's program in a way that reflects directly time t by switching the space of the integral from born date b to calendar time t.

W N C = ∞ 0 N 0 u (c (t, τ )) e -(β-ρ)τ dτ e -ρt dt (21)
Differently from [START_REF] Calvo | Optimal Time-Consistent Fiscal Policy with Finite Lifetimes[END_REF], what happens if there is an abrupt event probability in this economy ? The basic way to reformulate this problem is the following

W C = Eb b -N N 0 u (c (t, τ )) e -(β-ρ)τ dτ e -ρb db + e -ρ bϕ S b ( 22
)
where b is the uncertain date at which a catastrophic event occurs. This is equivalent to say that there is an "unlucky" generation which born with an abrupt event. In case where an abrupt event does not happen, we have b → ∞. Then, the maximization program reduces to W N C . In this problem, we seek to maximize expected future utility stream. With integration by parts of the equation ( 22), we find the deterministic equivalent of the expected value of the uncertainty problem in ( 22),

W C = ∞ 0 f b b -N N 0 u (c (t, τ )) e -βτ dτ e -ρt dt + e -ρ bϕ S b d b = ∞ 0 f b b -N N 0 u (c (t, τ )) e -βτ dτ e -ρt dt d b A + ∞ 0 f b e -ρ bϕ S b d b B
Note that probability distribution and density function are

f (b) = θe -θb and F (b) = 1 -e -θb
Integrating by parts A and keeping in mind the definition (1) and fundamental theorem of calculus yield

dX = f b =⇒ X = b -N f (s) ds Y = b -N N 0 u (c (t, t -b)) e -βτ dτ e -ρb db =⇒ dY = N 0 u c t, t -b e -βτ dτ e -ρ b
Using Y dX = XY -XdY yields and inflicted penalty

A = b -N f (s) ds b -N N 0 u (c (t, t -b)) e -βτ dτ e -ρb db ∞ -N - ∞ -N F b N 0 u c t, t -b e -βτ dτ e -ρ bd b Recall that ∞ -N f (s) ds = 1. The part A leads to ∞ -N (1 -F (b)) N 0 u (c (t, t -b)) e -βτ dτ e -ρb db
We can reformulate part B in the following manner

B = ∞ -N f b e -ρ b N 0 u (c min ) -ψ S t, t -b e -βτ dτ d b
Since the natural resource stock is independent of age τ12 , one can simplify the term in integral,

B = ∞ -N f b N 0 u (c min ) -ψ (S (t)) e -βτ dτ d b
Then, the social planner's program becomes ;

A+B = W = ∞ -N (1 -F (b)) N 0 u (c (t, t -b)) e -βτ dτ + f (b) N 0 u (c min ) -ψ (S (t)) e -βτ dτ e -ρb db
By changing the space of the integral by born date b to calendar time t , using probability distribution, density function, and taking the overall sum A + B, we have objective function to be maximized by the social planner,

∞ 0 [U (c (t)) + θϕ (S (t))] e -(ρ+θ)t dt
where the post-value function can be written after integration of the integral

ϕ (S (t)) = Z 2 (u (c min ) -ψ (S (t)))
where

Z 2 = 1-e -(β-(ρ+θ))N (β-(ρ+θ))

Aggregate Utility from OLG model

From equation (7) and using ( 8)

c (t, τ ) = c (t, 0) e -(β-(ρ+θ))τ σ (23)
Aggregating the consumption of all generations gives

c (t) = N 0 c (t, τ ) dτ = N 0 c (t, 0) e (β-(ρ+θ)) dτ (24)
This yields

c (t) = c (t, 0) σ 1 -e -(β-(ρ+θ))N σ (β -(ρ + θ)) (25)
In order to ease the presentation, we denote

Z 1 = σ 1-e -(β-(ρ+θ))N σ (β-(ρ+θ))
. Plugging equation ( 23) in aggregate utility 26gives

U (c (t)) = N 0 c (t, τ ) 1-σ -c 1-σ min 1 -σ e -(β-(ρ+θ)) dτ (26)
Plugging equation ( 23) in aggregate utility 26gives

U (c (t)) = N 0 c (t, 0) σ 1-e -(β-(ρ+θ))N σ (β-(ρ+θ)) 1-σ 1 -σ e -(β-(ρ+θ)) dτ (27)
Calculating the integral and using equation (25) gives

U (c (t)) = c (t) 1-σ -c 1-σ min 1 -σ (Z 1 ) σ

Proof of Proposition 1

We analyze an economy without catastrophic event to show how catastrophic event probability could have different implications on local stability of the economy. To see this, we equalize catastrophe probability to zero. The maximization program is indeed the reduced form of the social welfare function.

W N C = ∞ 0 U (c (t)) e -ρt dt (28) subject to    K (t) = f (K (t)) -δK (t) -c (t) Ṡ = G (S) -γf (K (t)) (29)
The differential system describing dynamics of the economy is

               K = f (K) -δK -Z 1 λ -1 σ Ṡ = G (S) -γf (K) λ = ρλ -(f K -δ) λ + µγf K μ = ρµ -µG S (30)
In this version of the model, an important aspect is that social planner does not face a trade-off mentioned in the text. She takes into account only the negative effect of capital accumulation to the environment. Then, environment is just served as a sink and does not represent any amenity value.

In an economy without catastrophic events, some straightforward calculations allow us to find analytically the steady state equilibrium.

               S * = 1-ρ 2 K * (S) = (1-S * )S * γA λ * (S) = 1 M1 N -(A-δ)(1-S * )S * γA µ * (S) = λ * (S)[(A-δ)-ρ] γA (31)
Then, just recall the det (J) becomes the following one when there is not an abrupt event

det (J) = 2 (Aγ) 2 (A -δ) (A -δ) -ρ 1 -ρ 2 1 - 1 -ρ 2 > 0 K = (A -δ) (ρ -(A -δ)) < 0
Then, it is easy to remark that the aggregation term Z 1 cancels out when there is not an abrupt event probability. It is also obvious that K does not depend on Z 1 neither.

Proof of Proposition 2

The Jacobian matrix of the differential system with using given functional forms is

J =       A -δ 0 Z1 σ λ -1 σ -1 0 -Aγ (1 -2S) 0 0 0 0 (ρ + θ) -(A -δ) Aγ 0 2µ + θ ψω2 S 2 0 (ρ + θ) -(1 -2S)       (32)
Following [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF], the characteristic polynomial associated with Jacobian is

υ 4 -trJυ 3 + b 2 υ 2 -b 3 υ + det (J) = 0 (33)
where b 2 and b 3 are the sum of second and third order minors of Jacobian respectively. We have

trJ = 2 (ρ + θ) and -b 3 + (ρ + θ) b 2 -(ρ + θ) 3 = 0 ( 34
)
υ i = (ρ + θ) 2 ± ρ + θ 2 2 - Ω 2 ± Ω 2 -det (J) (35)
The sum of the determinants of sub-matrices of Jacobian can be specified

Ω = A -δ Z1 σ λ -1 σ -1 0 (ρ + θ) -(A -δ) + (1 -2S) 0 2µ + θ ψω2 S 2 (ρ + θ) -(1 -2S) + 2 0 0 0 Aγ (36)
Then we have

Ω = (A -δ) [(ρ + θ) -(A -δ)] + (1 -2S) [(ρ + θ) -(1 -2S)] (37)
and

det (J) = [(ρ + θ) -(1 -2S)] [(A -δ) (1 -2S) [(ρ + θ) -(A -δ)]] + (Aγ) 2 Z 1 λ -1 σ -1 σ 2µ + θ ψω 2 S 2 (38)
For the possibility of a Hopf bifurcation, Ω should be positive. In this framework, this one is possible when following condition is ensured

ρ + θ > G S > 0 (39)
In this case, we can observe that an economy in which all levels of the natural resource stock exceeding maximum sustainable yield is stable. A Hopf bifurcation in a 4x4 dimension system occurs when two of eigenvalues have only imaginary parts. This means that real part of these two eigenvalues crosses zero for some parameters. More precisely, the derivative of the real part of eigenvalues with respect to chosen bifurcation parameter is non-zero. The necessary condition to have Hopf bifurcation

det (J) - Ω 2 2 > 0 (40) det (J) - Ω 2 2 -(ρ + θ) 2 Ω 2 = 0 (41)
are necessary and sufficient in order to have four complex eigenvalues and two having only imaginary parts. A saddle-path equilibrium is a possible outcome if and only if det (J) > 0 and Ω < 0 (42)

Besides, the inequality 0 < det (J) < 1 2 Ω 2 is sufficient in order to have all eigenvalues with real parts, which implies local monotonicity.

Proof of Proposition 4

The Jacobian matrix of the differential system with abatement activity at steady state

J =       A -δ 0 Z1 σ λ -1 σ -1 + 1 1-α α µ λ 1 1-α λ -1 -1 1-α α µ λ 1 1-α µ -1 -Aγ (1 -2S) 0 0 0 0 -α 1-α α µ λ α 1-α λ -1 α 1-α α µ λ α 1-α µ -1 0 2µ + θ ψω2 S 2 0 (ρ + θ) -(1 -2S)       (43)
The sum of the determinants of sub-matrices and determinant of Jacobian are

Ω = (A -δ) [(ρ + θ) -(A -δ)] + (1 -2S) [(ρ + θ) -(1 -2S)] -2µ + θ ψω 2 S 2 αM α (1 -α) µ (44)
and

det (J) = [(ρ + θ) -(1 -2S)] [(A -δ) (1 -2S) [(ρ + θ) -(A -δ)]] + (Aγ) 2 Z 1 λ -1 σ -1 σ 2µ + θ ψω 2 S 2 (45)
The proof is easy to follow.

= (A -δ) [(ρ + θ) -(A -δ)] + 1 -2S A (ρ + θ) -1 -2S A - αM α (1 -α) 1 -(ρ + θ) S A >0
Notice that M depends only on constant parameters in the long run. The presence of abatement activity decreases the sum of sub-matrices of Jacobian Ω, which makes less likely that unstable spirals and cycles occur. To see the effect of a higher natural stock level on the determinant of the Jacobian matrix, we look at the first derivative with respect to S ∂ (det (J))

∂S = ((ρ + θ) -(A -δ)) <0 2 (A -δ) 2 1 -2S A -(ρ + θ) - (Aγ) 2 σ Z 1 2 µ λ + θ ψω 2 S 2 λ ∂λ S A ∂S <0 + ((A -δ) -(ρ + θ)) (1 -(ρ + θ)) γAS 2 >0 > 0
The two last terms are unambiguously negative and positive respectively. The determinant increases unambiguously with respect to S if 2 1 -2S A -(ρ + θ) < 0.

Proof for complementarity effect between different time periods

This proof aims to show the existence of the complementarity of preferences over time within this framework. The objective of the proof is to point out the fact that complementarity over time vanishes when waste rate γ tends to zero.

We write down the objective function in form

J [c (.)] = ∞ 0 u (c (t) , S (t)) e -(ρ+θ)t dt (46)
In order to ease the calculations, we don't take the quadratic form for regeneration function G (S) but simply take a linear regeneration function G (S) = mS. 13 At this point, [START_REF] Wirl | CYCLICAL STRATEGIES IN TWO-DIMENSIONAL OPTIMAL CONTROL MODELS: NEC-ESSARY CONDITIONS AND EXISTENCE[END_REF][START_REF] Wirl | A new route to cyclical strategies optimal control models[END_REF] show that the occurrence of limit cycles does not have any link with the form of regeneration function. The author also shows that a linear regeneration function can also generate limit cycles. Our proof also supports this idea by showing the existence of complementarity of preferences over time. Then, the use of a linear regeneration function does not enter in a conflict with the aim of the proof. We can express physical and natural capital from constraints

S (t) = e mt γA ∞ t K (s) e -ms ds (47) K (t) = e (A-δ)t ∞ t c (s) e -(A-δ)s ds (48)
To understand the complementarity effects between state variables and control and state variables, we refer to Volterra derivatives (Heal andRyder (1973), Dockner and[START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF]). This requires to look at the marginal rate of substitution between different time periods t 1 ,t 2 , t 3 etc. For example, the marginal utility at time t 1 is J [c (.) , t 1 ] which is a Volterra derivative. To sum up, a small incremental increase of consumption at the neighborhood of time t 1 can be calculated by using Volterra derivatives. The concept of Volterra derivative is useful to show how a change in consumption at a given date shifts the allocation of consumption between other dates. The marginal rate of substitution between consumption at dates t 1 and

t 2 R [c (.) , t 1 , t 2 ] = J [c (.) , t 1 ] J [c (.) , t 2 ] (49)
In order to see the effect of an incremental change of consumption near date t 3 , we take the Volterra

derivative of R [c (.) , t 1 , t 2 ]. R [c (.) , t 1 , t 2 ; t 3 ] = J [c (.) , t 2 ] J [c (.) , t 1 , t 3 ] -J [c (.) , t 1 ] J [c (.) , t 2 , t 3 ] (J [c (.) , t 2 ]) 2 (50)
If R [c (.) , t 1 , t 2 ; t 3 ] > 0, an incremental increase of consumption at date t 3 shifts consumption from t 2 to t 1 , in which case there is complementarity between t 1 and t 3 . This represents a distant complementarity.

If R [c (.) , t 1 , t 2 ; t 3 ] < 0, the preferences shift from t 1 to t 2 where two neighboring dates hold for adjacent complementarity. Taking the derivatives of ( 46) , ( 47) and ( 48 where f c is the derivative of S with respect to c. Regarding the equation ( 47), we know that S is a function of c since K depends on c . For example, when consumption changes marginally, capital accumulation changes and consequently the trajectory of natural capital accumulation. In this sense, we can simply say S = f (c) but we can not know the form of this function analytically out of steady state (see [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF]). Note also that since we have an additive objective function, u ss does not depend directly on c in our model.

For the sake of simplicity regarding the analysis, we can restrict our attention to a constant investment path (i.e steady state) similar to [START_REF] Heal | Optimal Growth with Intertemporally Dependent Preferences[END_REF] and [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF]. Using equations ( 47) and ( 48), we can write mS * = γK * = γc A . It is obvious that with this simplification, we can find the form of f c which is equal to γ mA > 0. We write simplified form of equations ( 51) and ( 52), J [c (.) , t 1 ] = e -(ρ+θ)t1 u c (c (t 1 ) , S (t 1 )) + f c (c (t 1 )) u s (c (t 1 ) , S (t 1 )) (ρ + θ -m) J [c (.) , t 1 , t 2 ] = f c (c (t 1 )) f c (c (t 2 )) e -m(t1+t2)-(ρ+θ-2m)t2 u ss ρ + θ -2m With all these elements, it is easy to express the effect of marginal increase of consumption at date t 3 on the marginal rate of substitution between t 1 and t 2 , R [c (.) , t 1 , t 2 ; t 3 ] = [α (t 3 -t 1 ) -α (t 3 -t 2 )] where 0 < t 1 < t 2 . and u ss < 0. We know that the date t 2 and t 3 is placed after t 1 but we don't know the order of dates for t 2 and t 3 . In order to understand the effect of a small increase near date t 3 , we claim that date t 3 is situated before t 2 . Otherwise, it is evident on the expression above that there is no effect of a variation of consumption at date t 3 on the marginal rate of substitution of consumption between dates t 1 and t 2 . This makes sense because when we are near date t 3 , all decisions at t 1 and t 2 are already made. Similar to [START_REF] Heal | Optimal Growth with Intertemporally Dependent Preferences[END_REF], we write α (t) = e -(ρ+θ-m)t f or t > 0 α (t) = e mt f or t < 0

We observe that R < 0 which means that there exists an adjacent complementarity between dates t 2 and t 3 .

t 3 < (ρ + θ -m) t 1 -mt 2 (ρ + θ -2m) (53)
Note that right hand side of the inequality increases when θ increases. Limit cycles can appear in the model with both adjacent and distant complementarity, which is in line with the results of [START_REF] Ej Dockner | On the optimality of limit cycles in dynamic economic systems[END_REF].

It is important to point out the presence of waste in the model. It is obvious to see

lim γ→0 R = 0 ( 54
)
lim θ→0 R = 0 ( 55
)
lim θ→∞ R = 0 (56)

The limit of the condition (54) means that when the rate of waste coming from physical capital accumulation is lower, the complementarity over time becomes lower as well. Also, since stated in the text, when the abrupt event probability is too low or too high, the complementarity over time vanishes as well. This gives a mathematical support to see how waste from physical capital is a source of limit cycles in the economy.

The limit cycles occur at the presence of waste but for which parameter sets, does it occur ?

Variation of Aggregation Term w.r.t β

We differentiate the aggregation term with respect to x = β -(ρ + θ)

∂Z 1 ∂x = σ e -xN σ 1 + N σ -1 x 2 < 0
The variation is unambiguously negative when individual discount rate is higher than social planner's discount rate. e -xN σ 1 + N σ attains 1 when x = 0 after it decreases monotonically.

6 References

Figure 1 :

 1 Figure 1: Allocation and aggregation of consumption

Figure 2 :

 2 Figure 2: Allocation of consumption across generations

Figure 4 :

 4 Figure 4: Steady state level of natural resource stock for variations in individual discount rate β

  Figure (4) shows the steady state values of natural resource stock with respect to individual discount rate β. Proposition 1. In an economy with abrupt event probability, a Hopf bifurcation is possible for different critical parameters β. 6

FigureFigure 5 :

 5 Figure 5: Limit cycles on a phase plane (K, S) with bifurcation parameter β

Figure 7 :

 7 Figure 7: det (J) and Ω in benchmark (blue) and augmented model

  ), J [c (.) , t 1 ] = e -(ρ+θ)t1 u c (c (t 1 ) , S (t 1 )) + f c (c (t 1 )) e -mt1 ∞ t1 e -(ρ+θ-m)t u s (c (t) , s (t)) dt (51) J [c (.) , t 1 , t 2 ] = f c (c (t 1 )) f c (c (t 2 )) e -m(t1+t2) ∞ t2e -(ρ+θ-2m)t u ss (c (t) , s (t)) dt (52)

  γ 2 c(t1)c(t2) A 2 u SS (c(t1),S(t1)) (ρ+θ-2m) u c (c (t 1 ) , S (t 1 )) + γc(t 1 ) A us(c(t1),S(t1)) (ρ+θ-m)

  Figure 6: Limit cycles across time for utility function What model stipulates is that not only the level of natural resource stock or an abrupt event as a fact would violate Sustainable Development criterion but even a probability of an extreme event and individual preferences could compromise sustainability. Figure (

	c t
	2471.15
	2471.10
	2471.05
	2471.00
	2470.95
	2470.90
	50	100	150	200	250	300	350	t

Proposition 4. Sustainability Development criterion is violated for a combination of parameters β.

  Let S A > S B where S A is the steady state level of natural stock with abatement activity and S B stands for the steady state level of natural stock for the benchmark model without abatement activity. By replacing steady state value of µ in (44), we can reformulate the sum of sub-matrices of Jacobian Ω

Contributions made in intergenerational equity discussions argue that ILA framework as a utilitarian social welfare function corresponds to a different generation at each point of time (see[START_REF] Schneider | Trading off generations: Equity, discounting, and climate change[END_REF]).

Note that cycles around steady state are optimal.

The minimum level of consumption can be thought as a subsistence consumption. The reduction of consumption can be imposed by the social planner as a political response to appease the public outrage. (SeeTsur and Zemel (2016a) for a detailed discussion.)

See Ayong Le Kama (2001) and[START_REF] Wirl | Sustainable growth, renewable resources and pollution: Thresholds and cycles[END_REF] for a similar model.

When we choose β as a bifurcation parameter, we use following set of parameters ; ρ = 0.85, N = 60, A = 1, ω 2 = 0.1,

Note that we can not put exactly β = ρ + θ since the aggregation term Z 1 and Z 2 turn out to be indeterminate. In the numerical analysis, we also skip the point where individual discount is exactly equal to the adjusted discount rate.

The phase diagram is plotted for the first critical bifurcation point where β = 0.925644.

Note thatBosi and Desmarchelier (2016a) use pollution stock instead of natural capital stock in their model.

Note that S t, t -b = S (t)

Otherwise, in case of the use of logistic growth function for natural regeneration, one should deal with Riccati differential equation which yields very tedious calculations.


Etienne De Saint Germain
email: etienne.de-saint-germain@argon-consult.com

Vincent Leclère
email: vincent.leclere@enpc.fr

Frédéric Meunier
email: frederic.meunier@enpc.fr

A stochastic multi-item lot-sizing problem with bounded number of setups

Keywords: Lot-sizing, Stochastic Optimization, Sample Average Approximation, Simulation

Within a partnership with a consulting company, we address a production problem modeled as a stochastic multi-item lot-sizing problem with bounded numbers of setups per period and without setup cost. While this formulation seems to be rather non-standard in the lot-sizing landscape, it is motivated by concrete missions of the company. Since the deterministic version of the problem is NP-hard and its full stochastic version clearly intractable, we turn to approximate methods and propose a repeated two-stage stochastic programming approach to solve it. Using simulations on real-world instances, we show that our method gives better results than current heuristics used in industry. Moreover, our method provides lower bounds proving the quality of the approach. Since the computational times are small and the method easy to use, our contribution constitutes a promising response to the original industrial problem.

INTRODUCTION

Fixing the production level for the forthcoming week is a basic decision to be taken when managing a production line. Usually, a demand has to be satisfied at due dates but the limited capacity of the line prevents last minute production. On the other hand, too early productions may lead to unnecessary high inventory costs. The challenge of this kind of problems, known as lot-sizing problems in the operational research community, consists in finding a trade-off between demand satisfaction and holding costs. This is a well-studied topic, with many variations (deterministic/stochastic, single/multi item, etc.). Recent surveys have been proposed: see (Gicquel et al., 2008,Quadt and[START_REF] Quadt | Capacitated lot-sizing with extensions: a review. 4OR[END_REF] for the deterministic version and [START_REF] Mula | Models for production planning under uncertainty: A review[END_REF][START_REF] Aloulou | A bibliography of non-deterministic lot-sizing models[END_REF][START_REF] Díaz-Madroñero | A review of discrete-time optimization models for tactical production planning[END_REF] for the stochastic version. When several references can be produced on a same line -the so-called multi-item lot-sizing problem -, the capacity is often all the more reduced as the number of distinct references produced over the current week is high. Indeed, changing a reference in production stops the line for a moment. This additional capacity reduction is usually modeled by setup costs contributing to the total cost.

The present work introduces a stochastic multiitem lot-sizing problem met by the authors within a partnership with a consulting company. A nonstandard feature of the problem is that the capacity reduction due to reference changes is not modeled by setup costs but instead by an explicit upper bound on the total number of references that can be produced over a week. According to the consulting company, many clients aim at minimizing mainly their inventory costs while keeping the number of distinct references produced over each week below some threshold. This is essentially because, contrary to inventory costs, setup costs are hard to quantify and a maximal number of possible changes per week is easy to estimate. To the best of authors' knowledge, the problem addressed in the present work is original and such a bound on the number of distinct references produced over a week has not been considered by academics yet, with the notable exception of [START_REF] Rubaszewski | A capacitated lot-sizing problem with limited number of changes[END_REF] but, contrary to our problem, their bound is an overall bound for the whole horizon and they still consider setup costs.

We consider two versions of the problem. In a first one, backorder costs are present and the objective is to minimize the overall inventory costs (sum of holding and backorder costs). In a second version, there are only holding costs and there is a service level constraint to be satisfied over the whole horizon. We propose for these two versions a method that can be easily used and maintained in practice. The efficiency of the method is proved via extensive numerical experiments on real industrial data. In particular, we compare our results with current heuristics used in industry and with lower bounds.

PROBLEM FORMULATION AND MODEL

The assembly line produces a set R of references over T weeks. The number of distinct references produced over a week cannot exceed N. There is also an upper bound on the total week production (summed over all references). We normalize all quantities so that this upper bound is equal to 1.

The production must satisfy a random demand. The demand of reference r over week t is a random parameter d r t , whose realization is known at the end of week t. When production of a reference r is not used to satisfy the demand, it can be stored but incurs a unit holding cost h r > 0 per week. When a demand is not satisfied by the production of the current week or by inventory, it can be satisfied later but incurs a unit backorder cost γ per week for some coefficient γ > 0. Note that there is no setup cost, as discussed in the introduction. For each reference r, there is an initial inventory s r 0 ∈ R + . At the beginning of each week, before the demand of each reference is revealed, the production of the week has to be fixed. The objective is to minimize the total expected inventory cost (holding cost plus backorder cost) over the whole horizon of T weeks.

Regarding randomness, we assume that for any r and t, realizations of (d r t , d r t+1 , . . . , d r T ) have finite expectation and can be efficiently sampled, knowing a realization of d r t-1 . Before turning to the modeling of the problem, we make two remarks. The first one is about its position in the literaure. Our problem (in its deterministic version) is a variation of the Capacitated Lot-Sizing Problem (see [START_REF] Freeman | The capacitated lot sizing problem: a review of models and algorithms[END_REF] for a review of models and algorithms). The single differences with the usual basic version of that problem is that there is no setup cost and that there is this upper bound on the total number of distinct references that can be produced over a week. The second remark is about the somehow non-conventional assumption that the demand is revealed at the end of the week. There are indeed industrial cases where this occurs, like in computer manufacturing. However, the methods developed in this paper can be easily adapted to the case where the demand is revealed at the beginning of the week.

The production problem at the beginning of week t can be modeled as the following stochastic program (S):

min E T ∑ t =t ∑ r∈R (h r sr t + γb r t ) s.t. s r t = sr t -b r t t ∈ t, T , r ∈ R s r t = s r t -1 + q r t -d r t t ∈ t, T , r ∈ R ∑ r∈R q r t ≤ 1 t ∈ t, T q r t ≤ x r t t ∈ t, T , r ∈ R ∑ r∈R x r t ≤ N t ∈ t, T x r t ∈ {0, 1} t ∈ t, T , r ∈ R q r t , sr t , b r t ≥ 0 t ∈ t, T , r ∈ R σ(q r t ) ⊂ σ (d r t , . . . , d r t -1 ) r ∈R t ∈ t, T , r ∈ R ,
where the variable sr t (resp. b r t ) models the inventory (resp. the backorder) of reference r at the end of week t , the variable q r t models the quantity of reference r produced over week t , and the variable x r t takes the value 1 if the reference r is produced over week t and 0 otherwise. All these variables are random. The last constraint of the program, written as a measurability constraint, means that the values of the variables q r t can only depend on the values taken by the demand before time t (the planner does not know the future). The reader would have probably expected similar constraints for the other variables, but it is easy to see that adding them changes neither the optimal value, nor the production decisions. Note that in this program the inventory s r t-1 is not a variable and has a known deterministic value at the beginning of week t.

A feasible solution (q r t ) t ≥t, r∈R of (S) provides a deterministic production (q r t ) r∈R for the current week t.

FOCUS ON THE DETERMINISTIC VERSION

The deterministic version of (S) is obtained by removing the measurability constraint and by considering all d r t 's deterministic, denoted then simply by d r t . The variables are then also denoted with non-bold characters. The main message of this section is that this version is already difficult.

The deterministic version is NP-hard in the strong sense for any fixed N ≥ 3, since there is a straightforward reduction from 3-PARTITION. Given a set of 3n positive integers {a 1 , . . . , a 3n }, this latter problem consists in deciding whether this set admits a partition into n triples of same sum. It is a notorious strongly NP-complete problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. The problem 3-PARTITION can be reduced to the case N = 3 where each index in 1, 3n is a distinct reference r whose demand d r t is equal to 0 for all period t, except for the last period T = n where it is equal to It is also worth noting that the optimal value of the continuous relaxation is independent of N. This can be seen by considering the mixed integer program obtained by removing the constraints ∑ r∈R x r t ≤ N: any feasible solution of that program provides a feasible solution of the original program of same value, simply by redefining x r t = q r t for all t and all r if necessary. In an attempt of improving the quality of the continuous relaxation, one may consider the extended formulation with the binary variables y p t for each p ∈ R N and each t ∈ 1, T in place of the x r t 's (indicating the references produced on period t) but it can be shown with a nontrivial proof that it does not improve the quality of the continuous relaxation. Similarly, the extended formulation with the binary variables z r τ for each τ ⊆ 1, T (indicating the periods of production of the reference r) does not improve the quality of the relaxation (this time, with an easy proof).

METHOD

As explained in Section 3, the deterministic version of (S) is difficult. Therefore, we cannot expect a quick algorithm solving exactly the problem, and this holds especially for the full stochastic version. Moreover, one of the requests of the partner was to have an easy to understand method, which can be used and maintained in practice, with short computation times. We propose a two-stage approximation consisting in replacing the measurability constraint by

σ(q r t ) ⊂ σ(∅) r ∈ R σ q r t ⊂ σ d r t , . . . , d r T r ∈R t ≥ t + 1, r ∈ R ,
which provides a relaxation of the initial program: the production decisions for the current week t can still not depend on the future, but now the subsequent production decisions depend on the future demand. We denote this relaxation by (2SA). This approximation is a two-stage approximation as we distinguish between two levels of information over the uncertainty: production decisions for the first week are the first stage variables, while all other decisions are second stage variables. Three-stages or more generally multistage approximation would give better approximations of (S) but increases exponentially the number of variables. We chose for practicability reasons to stick to the two-stage approximation.

The (2SA) relaxation is then solved by a classical sample average approximation, see [START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF] for a presentation of the method. We build a set Ω of m scenarios sampled uniformly at random. Each of these scenarios is a possible realization of (d r t , d r t+1 , . . . , d r T ) for each r. The parameter m is fixed prior to the resolution.

We get the following mixed integer program (2SA-m), solved by any standard MIP solver.

min 1 m ∑ ω∈Ω T ∑ t =t ∑ r∈R h r sr t ,ω + γb r t ,ω s.t. s r t ,ω = sr t ,ω -b r t ,ω t ∈ t, T , r ∈ R , ω ∈ Ω s r t ,ω = s r t -1,ω + q r t ,ω -d r t ,ω t ∈ t, T , r ∈ R , ω ∈ Ω ∑ r∈R q r t ,ω ≤ 1 t ∈ t, T , ω ∈ Ω q r t ,ω ≤ x r t ,ω t ∈ t, T , r ∈ R , ω ∈ Ω ∑ r∈R x r t ,ω ≤ N t ∈ t, T , ω ∈ Ω x r t,ω = x r t r ∈ R , ω ∈ Ω q r t,ω = q r t r ∈ R , ω ∈ Ω x r t , x r t ,ω ∈ {0, 1} t ∈ t, T , r ∈ R , ω ∈ Ω q r t , q r t ,ω , sr t ,ω , b r t ,ω ≥ 0 t ∈ t, T , r ∈ R , ω ∈ Ω
At week t, the production is then set to be the solution (q r t ) r∈R found by the solver. The validity of this method for solving (2SA) is supported by the following proposition.

Proposition 1. The following three properties hold when m goes to infinity:

(i) The value of (2SA-m) converges almost surely to the optimal value of (2SA). (ii) For every m, we consider the values ( qr t,m , xr t,m ) r∈R of the decision variables for week t of an optimal solution of (2SA-m). Any limit point of these values is an optimal solution of (2SA). (iii) Let ε > δ > 0. Assume that the random demand

(d r t ) t ≥t,r∈R is such that ∃C, K, ∀u ∈ R, E[e u d ] ≤ Ce u 2 K . (1)
Denote by Q δ m (resp. Q ε ) the set of all possible values of ( qr t,m ) r∈R in a δ-optimal solution of (2SA-m) (resp. in an ε-optimal solution of (2SA)). Then for every α ∈ (0, 1), we have

P(Q δ m ⊆ Q ε ) > 1 -α for m large enough.
If the random demand d is bounded or Gaussian then it satisfies (1). (Some companies work with possibly negative demands. Assuming a normal distribution is hence not irrelevant.)

The proof of part (iii) relies on the following technical lemma.

Lemma 2. Consider g(d) = inf y∈Y G(y, d), where Y is non-empty and where the function G is non-negative and κ-Lipschitz with respect to d. If the random variable d satisfies (1), then g(d) also satisfies (1) with C = max{1, e G(y 0 ,0) C} and K = Kκ + G(y 0 , 0), for y 0 ∈ Y .

Proof. For u ≥ 0 we take y 0 ∈ Y yielding g(d) ≤ G(y 0 , 0) + κ d .

We then have E[e ug (d) ] ≤

E[e uG(y 0 ,0) e κu d ] ≤ C e K u 2 . For u ≤ 0, by nonnegativity of G we have g(d) ≥ 0, hence E[e ug (d) ] ≤ 1.

Proof of Proposition 1. Let Q t be the (bounded) set of feasible values for the first-stage variables q t = (q r t ) r∈R in (2SA). Denote by F(q t , d) the minimal cost of (2SA) that can be reached when the first stagevariables are fixed to q t ∈ Q t and the realization of the demand is d = (d r t ) t ≥t,r∈R . We introduce the map f : q t → E[F(q t , d)], which associates to a given choice of q t ∈ Q t the expected minimal cost, and similarly the map fm , which associates to a q Q t the minimal cost of (2SA-m) when the first-stage variables are set to this q t . The map f is continuous and Q t is compact. Thus f is bounded, and, by (Shapiro et al., 2009, Theorem 7.48), we have that ( fm (q t )) m∈Z + converges to f (q t ) uniformly on Q t . Then (i) and (ii) are direct consequences of (Shapiro et al., 2009, Theorem 5.3).

By Lemma 2, there exist K and C such that for any q t ∈ Q t , F(q t , d) satisfy (1). Consequently there exists σ > 0 such that for all q t , q t ∈ Q t , the random variable [F(q t , d)f (q t )] -[F(q t , d)f (q t )] is σsubgaussian, (see e.g., [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF]). Furthermore, for any demand d, the map F(•, d) is Lipschitzcontinuous on Q t . Then, according to (Shapiro et al., 2009, Theorem 5.18), for every α ∈ (0, 1), there exists M ∈ Z + such that for m ≥ M, we have

P(Q δ m ⊆ Q ε ) > 1 -α.

NUMERICAL EXPERIMENTS

Instances

C++11 has been chosen for the implementations and Gurobi 6.5.1 was used to solve the model on a PC with Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz and 8Go RAM.

The instances used are realistic and have been provided by a client of the partner. The client gave actually the figures of seven assembly lines but we give the results for only two of them: L2 and L6. The line L2 experiences overcapacity: the ratio expected demand over total capacity is smaller than 1 (actually equal to 0.7). The line L6 experiences undercapacity: the ratio expected demand over total capacity is larger than 1 (actually equal to 1.1). The horizon T is the typical one used in practice by this client, namely T = 13 weeks (a quarter). The demand is obtained via d r t = ( dr t + e r t ) + , where the dr t 's are historical data and where e r t is a generalized autoregressive process so that e r t = 0.25e r t-1 + 0.75ε t , where ε t is a Gaussian white noise process with zero mean and standard deviation equal to v dr t where v is the "volatility" and is chosen in {0.2, 0.5}. The initial inventory is set to s r 0 = 1 3 ( dr 1 + dr 2 + dr 3 ). The other parameters are provided in Table 1. In particular, for each value of v, we have considered three possible values of the unit backorder cost γ, which have been determined following a procedure described in Section 6. At the moment, we do not discuss further these values and take them as part of the input, as required by the problem formulation. The parameter C is the capacity of the line before normalization. (Recall that problem and the model have been formulated in Section 2 after normalization.) In the column hr , we indicate the range of the holding costs before normalization. We obtain the h r 's by dividing these costs by C.

The number m of scenarios used to solve (2SA-m) is fixed to 20, determined by preliminary experiments showing that it is a good trade-off between accuracy and tractability. The time limit of the solver has been set to 90 seconds.

Other heuristics

Our method is compared with three other heuristics.

The first heuristic is the deterministic version of (S), where the random demand is replaced by its expectation.

The second one, the lot-size heuristic, consists in determining before the first week once and for all a value * r for each reference r ∈ R . At time t, if the inventory of reference r is below a precomputed safety level, the quantity q r t is chosen so that the inventory of reference r exceeds the safety level of exactly * r . In case of capacity issues, the production is postponed and thus backorder costs appear. In addition, if some capacity issues are easily anticipated, the production of a reference r can be activated even if the inventory is not below the safety level.

The third one, the cover-size heuristic is almost the same, but instead of precomputing a fixed quantity for each reference, a duration T * r is fixed before the first week. When the inventory of reference r is below the safety level, the quantity q r t is computed so that the inventory of reference r exceeds the safety level of the expected demand for the next T * r weeks. The values * r and T * r are determined as follows. (T * r ) r∈R is actually chosen to be the optimal solution of the following convex program, which somehow considers the problem at a "macroscopic" level. (Similar convex programs in the same context have been considered in the literature; see [START_REF] Ziegler | Solving certain singly constrained convex optimization problems in production planning[END_REF] for example.)

min ∑ r∈R h r dr T r s.t. ∑ r∈R 1 T r ≤ N T r > 0 r ∈ R ,
where dr = E ∑ T t=1 d r t . The parameter * r of the lot-size heuristic is then set to dr T * r . The safety levels have been provided by the partner and are those used in practice.

The cover-size heuristic adapts the production to the realization of the demand, contrary to the lotsize heuristic. According to our partner, it makes the cover-size heuristic more suitable for situations with low short term volatility of demand or for overcapacitated lines, while the lot-size heuristic is expected to behave better with high short term volatility of demand or for undercapacitated lines.

Notice that the backorder costs are not taken into account at all for determining the values of the parameters * r and T * r . But playing with safety levels allows to prevent too large backorder costs. However, in real life it is usually the other way round: the company does not associate costs to backorder and aims at keeping the total amount of unsatisfied demand below some predetermined level. We come back to this point later in Section 6.

Results

The results are provided in Table 2. All quantities are in Me and are given with a confidence interval at 95%.

The column LB provides the lower bound obtained by the optimal value at time t = 1 of program (2SA-m) (with m = 1000 and a time limit of 24 hours for the solver).

The column 2SA-m is the cost of the method proposed in Section 4. (We remind the reader that we propose m = 20 in this case.) The three next columns provide the results for the three heuristics described in Section 5.2. For these four columns, between 27 and 30 runs have been used for each instance.

Complementary numerical results of these experiments are provided in Tables 3 and4, discussed in Section 6 where other criteria are considered. Table 3 can already by of interest since it deals with holding costs obtained for the L2 instances.

Comments

Our method clearly outperforms lot-size and coversize heuristics and is better than the deterministic approximation for all but one instance. By running our method instead of a usual heuristic at the beginning of each week, the inventory costs can be reduced often by more than 50%. For the instance L2_v20_13, the inventory costs have been divided by more than 6 (which corresponds to several Me). We may note that the cover-size heuristic behaves better than the lot-size heuristic on the L2 instances, while the converse holds on the L6 instances. It is in line with our partner's opinion regarding their behaviors with respect to the capacity of the line (see Section 5.2). Regarding their behaviors with respect to the volatility, it does not seem to be possible to draw any concrete conclusion.

The 2SA-m algorithm requires 90 seconds to output a solution, while lot-size and cover-size heuristics take less than a second and the deterministic approximation less than 10 seconds. The 2SA-m algorithm is thus slower, but note that 90 seconds to be run only once at the beginning of each week remains very short. Moreover, even with improved computers, the lot-size and cover-size heuristics and the deterministic approximation will not change their output (in all our experiments, Gurobi always found the optimal solution of the deterministic approximation). This is not the case for 2SA-m, which means that it would benefit from improved computational capacities.

Otherwise, complementary figures not provided here show that the deterministic approximation and 2SA-m use almost every available setups (between 97% and 100%) whereas lot-size and cover-size heuristics cannot take full advantage of this flexibility (only 84% to 98% of setups are used).

Our method dramatically reduces the inventory costs, while remaining a quite simple approach. It provides without any doubt a positive answer to the industrial request. Moreover, as explained in the next section, our method can be used to address the problem of minimizing the holding costs, while keeping the backorder at a reasonable level.

CONTROLLING THE FILL RATE SERVICE LEVEL

The problem

In practice, except when they are enshrined through contracts with the client, backorder costs can be hard to estimate. Thus, we consider here an alternative version of the problem, which often meets the objectives in industry. Everything remains the same, except that there is no backorder cost anymore, which means that we remove the corresponding part in the objective function of (S) and that the following constraint is added:

E ∑ r∈R ∑ T t=1 min(d r t , q r t + sr t-1 ) ∑ r∈R d r ≥ β, (2)
where β ∈ [0, 1] is the desired fill rate service level defined by the company, and where d r = ∑ T t=1 d r t . This constraint implies that, in expectation, the proportion of the production delivered on-time is at least β. Indeed, at week t, we deliver on-time the minimum between the demand and the sum of the real inventory at the end of week t -1 and of the production of week t.

First, note that this constraint does not make a difference between a one-week delay and a two-week delay. Second, this constraint is a global constraint to satisfy over the whole horizon: it does not depend on the current week. This is a matter of modeling point of view since we were not able to get a precise formulation from our partner and from its clients. A drawback of this formulation with respect to the one given in Section 2 is that this additional constraint cannot always be satisfied. In the next section, we propose a way to address this new constraint in an approximate way, while being able to provide solutions in any case.

Surrogate backorder costs

We address the alternative formulation by defining a surrogate backorder coefficient γ before the first week, with the idea to use the algorithm of Section 4 and to heuristically entice it to choose solutions satisfying constraint (2).

γ := max r∈R γ r γ r := P[d r ≤ q r (β)] P[d r > q r (β)] q r (β) := min q ∈ R + E min(d r , q) d r ≥ β (3a) (3b) (3c)
where d r = ∑ T t=1 d r t is the demand of reference r aggregated over time. Since d r is integrable, q r (β) is well-defined (we set 0 0 = β so that references with no demand would not impact the constraint). Computing an approximate value of q r (β) at an arbitrary precision can easily be performed by binary search.

To justify this choice, consider the original problem of Section 2 with only one reference and for a horizon of one week. Assuming no initial inventory, it takes then the form of the famous newsvendor problem (see e.g., [START_REF] Shapiro | Lectures on stochastic programming: modeling and theory[END_REF], Chapter 1))

min q≥0 E h r (q -d r ) + + γ r (d r -q) + , (4)
where γ r is a unit backorder cost specific to reference r. The next proposition means that with the right choice for γ r , the alternative version with a desired fill rate service level β (which takes the form of (3b) since h r > 0) is equivalent to the original one (4). Of course, it holds only for the case with one reference and a horizon of one week. Many references lead to several possible values γ r for the surrogate backorder cost. The choice in (3a) of taking the maximum of them for our actual γ aims at facing the additional uncertainty induced by the limited number of setups per week and the length of the horizon.

Proposition 3. Define γ r as in (3b). Then q r (β) is the smallest optimal solution to (4).

Proof. The aggregated production problem ( 4) is known to have the optimal solution q r * = F -1 d r (γ r /(γ r +h r )), where F -1 d r is the left-inverse of the cumulative distribution function of d r , i.e., F -1 d r (κ) = inf{q | P(d r ≤ q) ≥ κ}. Since we have set γ r = P[d r ≤q r (β)] P[d r >q r (β)] , we have q r * = inf{q | P(d r ≤ q) ≥ P(d r ≤ q r (β)}, which implies that q r * ≤ q r (β). Now if this inequality were strict, then it would mean that P(d r ∈ (q r * , q r (β))) = 0, which contradicts the minimality assumption in the definition of q r (β) (Equation (3c)).

Note that this formulation does not take into account the capacity constraint. It is therefore probably better suited to overcapacited production lines. Remark. If instead of controlling the fill rate service level, we want to control the cycle service level, defined as the probability of satisfying the whole demand, then we can choose

γ = β 1 -β max r∈R h r . (5)
Indeed, in this case, the optimal solution q r * of (4) satisfies P(q r * ≥ d r ) = β. Interestingly, Equation ( 5) does not depend on the distribution of the demand, which contrasts with the fill rate service level.

Numerical results

We discuss briefly the results of the experiments described in Section 5 in the new context of this section. Tables 3 and4 provide the numerical results: the first table shows the holding costs obtained when using the method proposed in this work; the second table shows the average (over the 30 runs) of the fill rate service level measured ex post, which we define by the realization of

∑ r∈R ∑ T t=1 min(d r t , q r t + sr t ) ∑ r∈R d r .
We plot on Figure 1 the values of the holding costs obtained for each of the 30 runs for the instance L2 with β = 95% and v = 0.5. On Figure 2, we plot the average of these values for three desired fill rate service levels (85%, 95%, and 98%).

Again, the inventory costs (which are in this version also the holding costs) are dramatically better than the ones obtained by the heuristics used by the clients of our partner. These two heuristic are however able to provide almost always very good fill rate service level but it comes at the price of very high holding costs

When β ∈ {95%, 98%}, the deterministic approximation reaches better holding costs than 2SA-m but worse fill rate service levels. This is especially true when the volatility is high (v = 0.5). Further experiments would probably be necessary to understand better how the deterministic approximation and 2SA-m compare but since they both use the surrogate backorder cost, their behavior regarding the fill rate service level shows that the way we compute it is relevant.

For high desired fill rate service levels (≥ 95%), the method 2SA-m behaves well: the fill rate service level measured in practice is close to the desired level. Since constraint ( 2) is an approximation of what is sought in practice, these results can be considered as positive and the method of Section 4 combined with our rule of thumb for determining γ forms an effective way to compute the production levels in real-world settings.

CONCLUSION

In this paper, we focus on a multi-item lot-sizing problem with a constraint not yet addressed by academic works, while being often met in practice: an explicit upper bound on the weekly number of setups. We model this problem as a stochastic program that includes this constraint and we propose a repeated two-stage approximation to solve it. Our method proves its efficiency on real-word instances and outperforms the heuristics currently used.

The costs in the previous model originate from holding and backorder. According to the consulting partner we work with, companies approach their production decisions by considering only holding costs, while trying to keep the backorder below some threshold. Building on the method proposed for the previous model, we propose a rule of thumb for determining surrogate backorder costs that bias the solutions towards the desired service levels. This rule relies on the famous newsvendor problem, which is the special case of our problem when there is one reference and one week. There exist generalizations of the newsvendor problem with more references or more weeks; it would be worth checking whether they could improve the way this surrogate cost is computed.

  cases N ≥ 4 are dealt with via a similar reduction.) The complexity status of the deterministic version when N = 1 or N = 2 seems to be a challenging open question.

Figure 1 :

 1 Figure 1: Results for the 30 runs of L2 with v = 0.5 and β = 95%

Figure 2 :

 2 Figure 2: Average results for L2 for various values of desired service level β and v = 0.5

Table 1 :

 1 Instance characteristics

	Instances			Instance characteristics
		|R |	max( dr t )	C	N	hr	v	γ
	L2_v20_13	21	4992	10562	7	35-61	0.2	13
	L2_v20_81							81
	L2_v20_203							203
	L2_v50_48						0.5	48
	L2_v50_154							154
	L2_v50_341							341
	L6_v20_3	22	8640	13299	8	16-23	0.2	3
	L6_v20_19							19
	L6_v20_55							55
	L6_v50_11						0.5	11
	L6_v50_42							42
	L6_v50_98							98

Table 2 :

 2 Results -Inventory costs (in Me)

	Instances	LB	2SA-m	Det.	Cover-size	Lot-size
	L2_v20_13	0.53	0.89 ± 0.03	1.17 ± 0.10	6.95 ± 0.17	7.79 ± 0.14
	L2_v20_81	0.94	2.29 ± 0.06	2.36 ± 0.07	8.12 ± 0.19	9.65 ± 0.14
	L2_v20_203	1.00	3.05 ± 0.07	3.25 ± 0.08	9.35 ± 0.29	10.99 ± 0.19
	L2_v50_48	0.97	2.73 ± 0.11	3.06 ± 0.21	8.03 ± 0.26	8.37 ± 0.21
	L2_v50_154	1.36	4.54 ± 0.20	5.06 ± 0.33	10.83 ± 0.53	11.20 ± 0.38
	L2_v50_341	1.51	5.91 ± 0.25	7.90 ± 0.66	15.17 ± 1.21	14.65 ± 0.77
	L6_v20_3	0.54	0.61 ± 0.01	0.70 ± 0.02	1.71 ± 0.08	1.74 ± 0.08
	L6_v20_19	1.41	1.81 ± 0.06	1.86 ± 0.06	3.51 ± 0.12	3.20 ± 0.08
	L6_v20_55	2.67	3.57 ± 0.24	3.71 ± 0.30	7.49 ± 0.39	6.24 ± 0.34
	L6_v50_11	1.33	2.00 ± 0.11	2.14 ± 0.12	3.42 ± 0.15	3.03 ± 0.13
	L6_v50_42	2.99	4.45 ± 0.53	4.48 ± 0.51	7.99 ± 0.62	6.57 ± 0.60
	L6_v50_98	6.13	8.29 ± 1.23	7.94 ± 1.04	16.34 ± 1.61	12.96 ± 1.45

Table 3 :

 3 Results -Holding costs (in Me)

	Instances

Table 4 :

 4 Results -Fill rate service level (in %)

	Instances