From lakes and glades to viability algorithms: Automatic classification of system states according to the Topology of Sustainable Management

TL;DR

This paper integrates viability theory with the Topology of Sustainable Management framework to visualize and analyze system states, introducing algorithms for better understanding complex environmental and socioeconomic models.

Contribution

It connects the Topology of Sustainable Management to viability theory, extending algorithms to compute capture basins and applying them to a coupled environmental-socioeconomic model.

Findings

Viability algorithms effectively visualize system state partitions.

Coordinate transformations address unbounded state spaces.

The approach enhances understanding of manageable dynamical systems.

Abstract

The framework Topology of Sustainable Management by Heitzig et al. (2016) distinguishes qualitatively different regions in state space of dynamical models representing manageable systems with default dynamics. In this paper, we connect the framework to viability theory by defining its main components based on viability kernels and capture basins. This enables us to use the Saint-Pierre algorithm to visualize the shape and calculate the volume of the main partition of the Topology of Sustainable Management. We present an extension of the algorithm to compute implicitly defined capture basins. To demonstrate the applicability of our approach, we introduce a low-complexity model coupling environmental and socioeconomic dynamics. With this example, we also address two common estimation problems: an unbounded state space and highly varying time scales. We show that appropriate coordinate transformations can solve these problems. It is thus demonstrated how algorithmic approaches from viability theory can be used to get a better understanding of the state space of manageable dynamical systems.

Equations

$$\dot{\text{mathcalboondox}x}$$

$$\forall \text{mathcalboondox}x\in \text{mathcalboondox}V\exists q_{\text{mathcalboondox}x}\forall t\in \mathbb{T}:q_{\text{mathcalboondox}x}(t)\in \text{mathcalboondox}Y.$$

$${\text{VN}}_{u_f}(\text{mathcalboondox}Y):={\text{Viab}}_{\{u_f\}}(\text{mathcalboondox}Y).$$

$$\begin{array}{rl}{\text{Capt}}_{\mathcal{U}}^{\text{mathcalboondox}Y}\left(\text{mathcalboondox}Z\right)=& \{{\text{mathcalboondox}x}_0\in \text{mathcalboondox}X\text{leavevmode}|\text{leavevmode}\exists \text{ solution}{q}_{{\text{mathcalboondox}x}_0}\exists T\in \mathbb{T}:\\ & q_{{\text{mathcalboondox}x}_0}(T)\in \text{mathcalboondox}Z\wedge \forall t<T:q_{{\text{mathcalboondox}x}_0}(t)\in \text{mathcalboondox}Y\}.\end{array}$$

$$\text{mathcalboondox}S:={\text{VN}}_{u_0}\left({\text{mathcalboondox}X}^{+}\right).$$

$$\text{mathcalboondox}M:={\text{Viab}}_{\mathcal{U}}\left({\text{mathcalboondox}X}^{+}\right).$$

$$\text{mathcalboondox}U:={\text{Capt}}_{\mathcal{U}}\left(\text{mathcalboondox}S\right)$$

$$\text{mathcalboondox}D:={\text{Capt}}_{\mathcal{U}}\left(\text{mathcalboondox}M\right)-\text{mathcalboondox}U.$$

$$\text{mathcalboondox}G={\text{Capt}}_{\mathcal{U}}^{{\text{mathcalboondox}X}^{+}}\left(\text{mathcalboondox}S\right)-\text{mathcalboondox}S.$$

$$\text{mathcalboondox}L=\text{mathcalboondox}U\cap \text{mathcalboondox}M-\text{mathcalboondox}S-\text{mathcalboondox}G.$$

$$\text{mathcalboondox}L{}_u$$

$$\text{mathcalboondox}L{}_l$$

$${\text{mathcalboondox}U}^{(+)/-}=(\text{mathcalboondox}U-\text{mathcalboondox}M)\cap {\text{mathcalboondox}X}^{+/-}.$$

$$\text{mathcalboondox}W=\text{mathcalboondox}M-\text{mathcalboondox}U.$$

$${\text{mathcalboondox}D}^{(+)/-}=\left(\text{mathcalboondox}D-\text{mathcalboondox}W\right)\cap {\text{mathcalboondox}X}^{+/-}.$$

$${\text{mathcalboondox}E}^{+/-}\subseteq {\text{mathcalboondox}X}^{+/-}-\text{mathcalboondox}U-\text{mathcalboondox}D,$$

$${\text{mathcalboondox}E}^{+/-}\subseteq {\text{Capt}}_{\mathcal{U}}({\text{mathcalboondox}E}^{-/+}),$$

$$\text{mathcalboondox}E={\text{mathcalboondox}E}^{+}+{\text{mathcalboondox}E}^{-}.$$

$$\Theta =\text{mathcalboondox}X-{\text{Capt}}_{\mathcal{U}}\left({\text{mathcalboondox}X}^{+}\right).$$

$${\rm Y}=\text{mathcalboondox}X-\text{mathcalboondox}U-\text{mathcalboondox}D-\text{mathcalboondox}E-\Theta .$$

$$S(\text{mathcalboondox}x,u)=\left\{\text{mathcalboondox}y\in {\text{mathcalboondox}Y}_h\text{leavevmode}|\text{leavevmode}\Vert \text{mathcalboondox}y-(\text{mathcalboondox}x+f(\text{mathcalboondox}x,u)\cdot \Delta t)\Vert \le h+\frac{Ml}{2}\cdot {\left(\Delta t\right)}^2\right\}.$$

$$K_{i+1}=\left\{x\in K_i\text{leavevmode}|\text{leavevmode}\exists u\in \mathcal{U}:S(x,u)\cap K_i\ne \varnothing \right\}$$

$$K_{i+1}^{\prime }=\left\{x\in {\text{mathcalboondox}X}_h\text{leavevmode}|\text{leavevmode}\exists u\in \mathcal{U}:S(x,u)\cap K_i^{\prime }\ne \varnothing \right\}.$$

$$\dot{\text{mathcalboondox}x}=f(\text{mathcalboondox}x)$$

$$\displaystyle\begin{array}[]{@{}r@{\ }c@{}c@{\enspace}l@{}}\Phi:&\mathbb{R}_{\geq 0}^{n}&\longrightarrow&\left[0,1\right)^{n}\\ &({\mathcalboondox{x}}_{i})&\longmapsto&(\frac{{\mathcalboondox{x}}_{i}}{{\mathcalboondox{x}}_{i,mid}+{\mathcalboondox{x}}_{i}}),\end{array}$$

$$\dot{\text{mathcalboondox}y}=F(\text{mathcalboondox}y):=$$

$$=$$

$$\dot{\text{mathcalboondox}y}=F(\text{mathcalboondox}y).$$

$$\dot{\text{mathcalboondox}y}=\tilde{F}(\text{mathcalboondox}y):=\frac{F(\text{mathcalboondox}y)}{\Vert F(\text{mathcalboondox}y)\Vert +ϵ}.$$

$$\Vert \tilde{F}(\text{mathcalboondox}y)\Vert =\frac{\Vert F(\text{mathcalboondox}y)\Vert }{\Vert F(\text{mathcalboondox}y)\Vert +ϵ}\approx 1.$$

$$\tilde{F}(\text{mathcalboondox}y)=\frac{F(\text{mathcalboondox}y)}{\Vert F(\text{mathcalboondox}y)\Vert +ϵ}\approx \frac{1}{ϵ}F(\text{mathcalboondox}y).$$

$${{\text{mathcalboondox}E}^{+}}_0$$

Full text

1 Introduction

Charting pathways to a sustainable future for humanity needs to account for different biophysical and social constraints. First, the concept of planetary boundaries developed by Rockstrom2009 (1) and extended by Steffen2015 (2) defines a set of biophysical boundaries that ensure human development under Holocene-like conditions on Earth. They set limits to a safe operating space in terms of the planet’s capacity to accommodate humanity’s metabolism. Raworth2012 (3) extended the framework to include social foundations, i.e. a decent livelihood for every member of society enabling just and equitable human development. Combined with the planetary boundaries, they define the safe and just operating space, within which a sustainable development is possible.

There is much research on refining the current definitions of the planetary boundaries, e.g. for freshwater Gerten2013 (4) and phosphorus Carpenter2011 (5), on downscaling Hayha2016 (6) and on extending them, e.g. to the terrestrial net primary plant production Running2012 (7). Only a limited number of studies focus on their interaction due to the system’s intrinsic dynamics Anderies2013 (8, 9).

Heitzig2016 (10) introduced a mathematical framework called the Topology of Sustainable Management (tsm) to study the interactions of boundaries and their consequences for a dynamical system. The framework can be used to analyze models that distinguish between a default dynamic and one or several management options. This reflects the observation that social-ecological systems are often governed by relatively stable rules. To change the resulting dynamics, an extra effort is necessary. Also, the dynamics can only be influenced to a certain degree by active management. Complete control is usually not possible. The boundaries and the dynamics induce a possibly complex partition of the model’s state space into different regions, corresponding to a hierarchy of safety levels. The regions differ qualitatively in how secure they are and how much management they need to either stay within a desirable region or reach it. The tsm framework allows for the qualitative analysis of a model and thus precedes other analyses like quantitative optimization.

In this paper, we introduce a variant definition of tsm based on viability theory, a subfield of control theory. This allows us to operationalize the tsm-framework for computations. Aubin and his collaborators Aubin1990 (11, 12, 13, 14) developed viability theory to address the problem of maintaining a dynamical system within a set of desirable states. Viability theory has been applied in many domains, e.g. economics Aubin1997 (15), wireless sensor node Kone2017 (16), language competition Chapel2010 (17), social networks Mathias2017b (18) and resilience modeling Martin2004a (19, 20, 21).

There is a growing literature using viability concepts to analyze sustainability in social-ecological models Delara2008 (22). For example, Baumgartner2009 (23) apply viability concepts to define conditions of strong sustainability under uncertainty. The review by Oubraham2018 (24) provides about 80 examples using viability in sustainability studies. In this context, viability theory has been applied for example to fisheries Cury2005 (25, 26, 27, 28), water Richter2003 (29), forests Mathias2015 (30), biodiversity Doyen2013 (31) and climate change Mathias2017a (32). Mathias2018 (33, 34) study the effect of knowledge infrastructure on safe operating spaces in stochastic social-ecological systems. Some research also focuses on the relationship between viability theory and other mathematical concepts of sustainability, for example the Maximin criterion Doyen2012 (35).

In contrast to tsm, viability theory does not distinguish between default dynamics and management options. Even though Heitzig2016 (10) already indicated relations between tsm and viability theory, they did not spell out the details. The analysis of the models in Heitzig2016 (10) was possible to do without computational analysis as they were all two-dimensional. However, to apply the framework to more complex models, we need to automatize the identification of tsm regions in the state space. Our definition of tsm differs from the original definition because it is based on the notion of reachability as defined in control theory instead of “safe reachability” as defined in Heitzig2016 (10). Also, we require a target set be reached in finite but arbitrary large time, instead of infinite time because it simplifies the computation and is more realistic. For computationally supported estimations of specific models these differences are usually not relevant.

Our operationalization of the tsm-framework is based on two main concepts from viability theory. The first one that we use is the viability kernel. It is defined as the largest set of states from where it is possible to keep the system’s trajectory within a desirable set and thus avoid the transgression of the set’s boundaries. The second main concept is the capture basin. It is defined with respect to a target set, i.e. a chosen set of states that one wants to reach. The capture basin is the largest set of states from which the target set can be reached. This concept is closely related to reachability in control theory.

Our variant definition of the tsm framework allows us to apply computational algorithms from viability theory to estimate different regions of the tsm partition. In particular, we use the Saint-Pierre algorithm Saint1994 (36) for the estimation of viability kernels. This algorithm is well-established and has been extended with classification methods based on machine learning deffuant2007 (37), to viability problems based on reachability Maidens2013 (38) and optimal control problems Bokanowski2006 (39, 22, 21). We also develop an extension of the Saint-Pierre Pierre algorithm to find specific tsm regions that can only be defined as implicit target sets. Furthermore, we propose a nonlinear, local time homogenization solving the problem of vastly differing time scales when estimating viability kernels and capture basins using the Saint-Pierre algorithm. It fully homogenizes the time scale of the dynamics, while keeping the major properties of the system invariant. With the help of these tools, we can visualize the shape and compute the volume of the tsm partition.

We demonstrate the applicability of the proposed methods by analyzing a three dimensional example model, which focuses on the interaction of stylized environmental and socioeconomic processes similar to other social-ecological models Wiedermann2015 (40, 41, 42, 43, 8, 44, 45). As part of the environment, we model the atmospheric carbon stock and its dynamics. The global economic output and a knowledge stock on renewable energy production, that can induce technological change, describe parts of the socioeconomic system. In our model, there are two boundaries. The first one is the planetary boundary on climate change Rockstrom2009 (1, 2) limiting the amount of atmospheric carbon. The second one is a social foundation for prosperity prescribing a lower threshold on yearly economic output. We analyze the model with two management options: (i) low growth and (ii) climate mitigation by inducing an energy transition.

We limit the model to key variables to keep its dimensionality low. Thereby, we ensure that it serves to illustrate the power of our automatic tsm estimation. Furthermore, we design the model to demonstrate that our methods can deal with common model properties that impede a direct application of computational algorithms such as vastly differing time scales and diverging dynamics. A recent study by strnad2020 (46) used the model to demonstrate an application of deep reinforcement learning for identifying sustainable management strategies.

The remainder of the article is structured as follows. First, we recall the central notions from viability theory in Section 2. These are used to derive the variant definition of tsm in Section 3. In Section 4, we show our means to deal with an unbounded state space and then introduce the nonlinear, local time-homogenization. Next, we shortly recall the Saint-Pierre algorithm and present its extension for implicitly defined target sets. In Section 5 we introduce the example model with its different managements and analyze it using the tools developed before. We close with a discussion and an outlook for future work.

2 Viability theory

In this section, we shortly introduce the notions of viability theory needed in this article, following Aubin1991 (12, 47). We start with a time-continuous ($t\in\mathbb{T}:=\mathbb{R}_{+}$) dynamical control system

[TABLE]

with ${\mathcalboondox{x}}\in{\mathcalboondox{X}}=\mathbb{R}^{n}$, the state space, and $u\in{\mathcal{U}}$, the set of all possible values for the control parameter $u$. We call $f$ the right-hand side function (RHS). Note that no dependency of ${\mathcal{U}}$ on ${\mathcalboondox{x}}$ is assumed for ease of notation, but the following could be extended for such a case in a straightforward manner.

A function $q_{{\mathcalboondox{x}}_{0}}\colon\mathbb{T}\rightarrow{\mathcalboondox{X}}$ is called a solution for an (arbitrary) initial condition ${\mathcalboondox{x}}_{0}\in{\mathcalboondox{X}}$ if there exists a measurable function of time $\pi\colon\mathbb{T}\rightarrow{\mathcal{U}}$, called policy, such that for any time $t\in\mathbb{T}$ the condition $\frac{dq}{dt}(t)=f(q(t),\pi(t))$ is fulfilled and $q_{{\mathcalboondox{x}}_{0}}(0)={\mathcalboondox{x}}_{0}$.

A viable set ${\mathcalboondox{V}}$ of a constraint set ${\mathcalboondox{Y}}\subseteq{\mathcalboondox{X}}$ is then defined as a set of initial conditions ${\mathcalboondox{x}}$ for which there exists a viable solution $q_{\mathcalboondox{x}}{}$ that stays within ${\mathcalboondox{Y}}$ forever:

[TABLE]

The largest viable set of ${\mathcalboondox{Y}}$ is called the viability kernel $\text{Viab}_{\mathcal{U}}({\mathcalboondox{Y}})$. The set of possible controls ${\mathcal{U}}$ is given as a subscript as we will distinguish different controls later. If one assumes only some constant control $u_{f}\in{\mathcal{U}}$, i.e. any solution has a policy $\pi(t)=u_{f}$, the corresponding viability kernel is called a viability niche

[TABLE]

The capture basin of a target set ${\mathcalboondox{Z}}\subseteq{\mathcalboondox{X}}$ is the subset of the state space in ${\mathcalboondox{X}}$ for which there exists a solution of Eq. 1 that reaches ${\mathcalboondox{Z}}$ in finite time

[TABLE]

In case no constraint set is given, the constrained set is assumed to equal the entire state space, i.e. $\text{Capt}_{\mathcal{U}}\left({\mathcalboondox{Z}}\right):=\text{Capt}^{\mathcalboondox{X}}_{\mathcal{U}}\left({\mathcalboondox{Z}}\right)$.

3 A variant of the topology of sustainable management based on viability theory

In this section, we present a variant definition of the tsm-partition based on viability theory. As viability theory, tsm builds on the distinction between a desirable and an undesirable region (set), called ”sunny“ and ”dark“ regions in the tsm framework. However, it additionally distinguishes between default dynamics and management options in the system. This distinction induces a fine partition of the state space into qualitatively different regions from which the desirable region can or cannot be reached. This also gives rise to multiple dilemmas. The decision tree in Fig. 1 shows the different regions of the tsm-partition. A state belongs to a specific region depending on whether it is feasible to stay within it or reach the desirable region and whether management options are needed to do so. From the perspective of viability theory, tsm can thus be seen as an extension building on the distinction of controls between default and management. For the full framework, we refer the reader to Ref. Heitzig2016 (10).

We recapture some central notions of tsm using an illustrative example of ducklings in Fig. 2. The water region represents the state space and the streamlines represent the dynamics. The ducklings can either swim with the flow (default dynamics) or struggle and swim against it (management). However, it is not possible to swim up a waterfall once they have dropped down. The desirable region, a safe environment providing enough food and nesting places for the ducklings, is on the left and the undesirable region full of predators on the right. In the following, we introduce the general concept for each region and then illustrate it with the ducklings example.

Our definitions of tsm-regions are based on a general control system as Eq. 1, for which we additionally require that a default control $u_{0}\in{\mathcal{U}}$ is separated out from the set of all possible controls ${\mathcal{U}}$. ${\mathcal{U}}_{m}={\mathcal{U}}-\{u_{0}\}$ is the set of manageable controls.111Here and in the following we use the lax difference and union notation with “$-$” and “$+$” for sets. Hence we call $f(x,u_{0})$ the default flow giving rise to the default dynamics and the dynamics corresponding to the manageable controls is called management options. Furthermore, we require a division of the state space into desirable ${{\mathcalboondox{X}}^{+}}\subseteq{\mathcalboondox{X}}$ and undesirable ${{\mathcalboondox{X}}^{-}}:={\mathcalboondox{X}}-{{\mathcalboondox{X}}^{+}}$ regions.

As already sketched in Heitzig2016 (10), the viability perspective allows us to base the partition of state space upon two basic regions, the shelters and the larger manageable region. In the shelters ${\mathcalboondox{S}}$, the system will stay in the desirable region forever without management, i.e. only by following the default dynamics described by the default control $u_{0}$. Thus, they are the safest regions in state space. In terms of viability theory, they are the viability niche of ${{\mathcalboondox{X}}^{+}}$ under $u_{0}$,

[TABLE]

For the ducklings in Fig. 2 the shelter is the nest, a place in which they can stay in a safe environment forever without swimming against the stream.

The manageable region ${\mathcalboondox{M}}$ is the part of the desirable set where one can stay by making use of all management options. With that, it is the viability kernel of ${{\mathcalboondox{X}}^{+}}$ when taking into consideration the full set of controls ${\mathcal{U}}$.

[TABLE]

In Fig. 2, this corresponds to all parts of the river in the desirable region on the left. However, there can be other parts of the desirable region that do not allow to stay in the desirable region forever, even with management (“remaining sunny downstream”, not depicted in Fig. 2).

In contrast to the definition here, Heitzig2016 (10) define the shelters as the “invariant open kernel” of ${{\mathcalboondox{X}}^{+}}$, which is slightly smaller than $\text{VN}_{u_{0}}\left({{\mathcalboondox{X}}^{+}}\right)$ because it does not contain its boundary. Likewise, they define the manageable region in terms of the “sustainability kernel” instead of the viability kernel, which again is a slight subset that does not contain the boundary. Analogous remarks apply to all the other regions defined below.

In tsm, the upstream is defined as the region from where it is possible to reach the shelter(s) with or without management. This can include leaving the desirable region for a limited time period. In viability theory, the capture basin formalizes the idea of being able to reach a certain set. Hence, the upstream is the capture basin of the shelter,

[TABLE]

The downstream are those regions from where the manageable region is reachable but it is not possible to reach the shelter(s). From the downstream it is possible to reach a region where one can stay within the desirable region forever using management. The downstream region is therefore the capture basin of the manageable region without the upstream,

[TABLE]

tsm partitions the state space into non-intersecting regions. This is why our definitions often need to remove certain sets so they are not classified twice. In the above definition, the removal of the upstream ensures that the shelter cannot be reached from the downstream.

As both the upstream and the downstream may have qualitatively different dynamics inside, tsm introduces a finer subdivision of these regions. The upstream is further divided into glades, lakes and the remaining sunny and dark upstream. Glades are the region from where it is possible to reach the shelter without leaving the desirable set. This means there is only a temporary need for management to stay in the desirable region forever. After reaching the shelter, the system will remain in the desirable set without management. In terms of viability theory, the glades are the capture basin of the shelter (excluding the shelter itself) when only considering the sunny region,

[TABLE]

In our ducklings example, the water around the nest is a glade.

From the lakes, the shelters are only reachable trough the undesirable region. Thus, they are the intersection between the upstream (shelter is reachable) and the manageable region with the shelters and the glades removed,

[TABLE]

Being located in a lake imposes a qualitative choice, the lake dilemma: Inside a lake one can choose between eventual safety and uninterrupted desirability. Deciding for the first requires crossing the undesired region to finally reach the shelter. But deciding for the second implies that there will always be a need for management options.

As a further refinement, let us distinguish two kinds of lakes, time-limited lakes ${\mathcalboondox{L}}{}_{l}$ and time-unlimited lakes ${\mathcalboondox{L}}{}_{u}$, which can be defined as

[TABLE]

Because time-unlimited lakes are the viability kernel of the lakes, there is no time pressure for a decision on the dilemma. In contrast, there is a fixed deadline in time-limited lakes. This leads to a stronger form of the dilemma called the pressing lake dilemma. In Fig. 2, the time-unlimited lake is the waterfall starting from the glade and the subsequent calmer waters. The ducklings can swim in the calmer waters as long as they want to. The time-limited lake is the following waterfall that splits into two streams. Here, the ducklings drop down the waterfall and they have only a moment to decide for the left or the right.

The rest of the upstream is split into the remaining sunny ${{\mathcalboondox{U}}^{(+)}}$ and the dark upstream ${{\mathcalboondox{U}}^{(+)}}$, depending on whether a state is in the desirable region or not. They can be defined as the intersection of the upstream without the manageable region with the desirable or undesirable region,

[TABLE]

The downstream is subdivided into backwaters, remaining sunny and dark downstream. The backwaters are those regions where it is possible to remain in the desirable region using management but the shelters are not reachable. The first condition means that they are a part of the manageable region and the second condition excludes that it is in the upstream, thus

[TABLE]

The backwaters can also be defined as the intersection of the manageable region and the downstream, ${\mathcalboondox{W}}={\mathcalboondox{M}}\cap{\mathcalboondox{D}}$. The calm waters in the left lower corner of Fig. 2 is a backwater because the ducklings can swim against the stream and stay inside the backwater but due to the waterfalls, they cannot reach the nest.

Within the capture basin of the backwaters, there are no further regions such as the lakes in the upstream. Analogously to the upstream, the rest of the downstream is divided into the remaining sunny ${{\mathcalboondox{D}}^{(+)}}$ and the dark downstream ${{\mathcalboondox{D}}^{-}}$, depending on whether they are in the desirable region or not. These regions are therefore the intersection of the downstream without the backwater with the desirable or undesirable region, respectively,

[TABLE]

If the desirable region can be reached over and over again one is inside the eddies ${\mathcalboondox{E}}$ that are divided into sunny eddies ${{\mathcalboondox{E}}^{+}}$ and dark eddies ${{\mathcalboondox{E}}^{-}}$ (these and the following regions are not depicted in Fig. 2). The metaphorical image behind the naming is that of a circular flow of which one segment is in the desirable and the other one is in the undesirable region. They are the maximal pair of sets fulfilling

[TABLE]

Section 4.4 and Appendix A provide details why eddies can be defined like this. Trenches are regions from which the desirable region cannot be reached ever again. In terms of viability theory, they are the regions that do not belong to the capture basin of the desirable region,

[TABLE]

To complete the partition of the state space, the abysses are defined as the rest of the state space not belonging to any region mentioned before,

[TABLE]

Inside the abysses, on can reach the desirable region a finite number of times only. There are sunny and dark abysses ${{\Upsilon}^{+/-}}={\Upsilon}\cap{{\mathcalboondox{X}}^{+/-}}$. This completes the main partition of the tsm framework.

4 Estimation

To estimate the tsm-partition, we use the definitions derived in Section 3 and the established Saint-Pierre algorithm Saint-Pierre2001 (48) to find the viability kernels and capture basins that our definitions build upon. We apply the Saint-Pierre algorithm because it is the standard algorithm in the field and operationalizes the most fundamental idea how to estimate viability kernels. Furthermore, we chose it as a first step for automatic TSM region detection because it is comparatively simple to implement. The Saint-Pierre algorithm always converges to the viability kernel under well-defined conditions and when decreasing adequately the time step and the space resolution. Other algorithms may be more efficient on specific examples but the conditions of their convergence to the viability kernel are more restrictive Krawczyck2011 (49, 39, 37, 50).

The Saint-Pierre algorithm is based on finitely discretizing the state space and then using local linear approximations of the dynamics. Hence, it is applicable to bounded state spaces only. In case the relevant domain of the state space is unbounded, we need to map it to a bounded space first. Also, vastly differing time scales might be problematic for the linear approximations, so there is a need to homogenize the time scales. In this section, we explain the basic ideas of the Saint-Pierre algorithm and introduce the two coordinate transformations to make it applicable. Finally, we extend the algorithm to estimate implicitly defined capture basins for identifying the eddies of the tsm-partition.

4.1 The Saint-Pierre algorithm

The Saint-Pierre algorithm Saint1994 (36) was developed in order to estimate the viability kernel of a control system Eq. 1.

It starts with a discretization ${\mathcalboondox{Y}}_{h}$ of the constraint set ${\mathcalboondox{Y}}$ where a point $x\in{\mathcalboondox{Y}}$ is at most at a distance $h$ of a point ${\mathcalboondox{y}}\in{\mathcalboondox{Y}}_{h}$222Often, a regular grid with resolution $h$ is chosen for this discretization, but this is by no means a necessity. We use a regular grid because we already apply a coordinate transformation (see Section 4.2) that best resolves the scales of interest in our model.. Furthermore, a small time step $\Delta t>0$ is chosen and it supposes that the set of controls ${\mathcal{U}}$ is discrete (if not, it is also discretized). It supposes that $f$ is $l$-Lipschitz and there exists an upper bound $M$ of ${\mathcalboondox{Y}}$.

For each point ${\mathcalboondox{x}}\in{\mathcalboondox{Y}}_{h}$ and for each control $u\in{\mathcal{U}}$, the algorithm computes the successors $S({\mathcalboondox{x}},u)$ of $x$ when applying control $u$ for a linearized, extended dynamics defined from $f$. Note that the linearization is done in the ${\mathcalboondox{x}}$ variable only. It is extended in the sense that the successors include all the points located in a ball around ${\mathcalboondox{x}}+f({\mathcalboondox{x}},u)\cdot\Delta t$. The successors $S({\mathcalboondox{x}},u)$ of $x$ when applying control $u$ are given by

[TABLE]

This extension of the dynamics guarantees that the algorithm described below converges to the actual viability kernel when $\Delta t$ and the resolution of the grid decrease to 0. Computing and storing all the successors for each point of the grid rapidly becomes computationally heavy when the dimensionality of the the state space is large and the grid resolution is small (this is an example of the famous curse of dimensionality).

Then, the algorithm builds a series of discrete sets (subsets of the grid ${\mathcalboondox{Y}}_{h}$) $K_{0}={\mathcalboondox{Y}}_{h},K_{1},...,K_{n}$ such that $K_{i+1}\subset K_{i}$, defined as follows:

[TABLE]

After a finite number of steps, the algorithm reaches a fixed point, i.e. $K_{n+1}=K_{n}$. The set $K_{n}$ is the viability kernel of ${\mathcalboondox{Y}}_{h}$ for the linearized, extended discrete dynamics. Saint1994 (36) shows that this set converges to the viability kernel of the continuous time dynamics when $\Delta t$ and $h$ tend to 0 appropriately. Note that the approximations are done from the exterior of the viability kernel: generally, the approximation includes points that do not belong to the actual viability kernel (their proportion decreases when the resolution of the grid decreases).

This algorithm has been extended by deffuant2007 (37) for using continuous sets $K_{i}$, using a machine learning algorithm that takes as input the points of the grid that belong to $K_{i}$ and the ones that do not, and derives an approximation of its boundary. This opens the possibility to represent continuous viability kernels that are defined more conveniently than a huge set of points.

A slight modification of the algorithm described above enables us to approximate the capture basin (see Saint-Pierre2001 (48) for details). We start with a discretization of the state space ${\mathcalboondox{X}}_{h}$, analogously to ${\mathcalboondox{Y}}_{h}$ above, and define the discretized target set ${\mathcalboondox{Z}}_{h}={\mathcalboondox{Z}}\cap{\mathcalboondox{X}}_{h}$ for a target set ${\mathcalboondox{Z}}$. Again, we create a series of discrete sets $K^{\prime}_{i}$ with $K^{\prime}_{0}={\mathcalboondox{Z}}_{h}$ and where the successors of all elements in $K^{\prime}_{i+1}$ are in $K^{\prime}_{i}$

[TABLE]

Again, after a finite number of steps, the algorithm reaches a fixed point, i.e. $K^{\prime}_{n+1}=K^{\prime}_{n}$, and $K^{\prime}_{n}$ is the capture basin of ${\mathcalboondox{Z}}_{h}$ in ${\mathcalboondox{X}}_{h}$ for the linearized extended discrete dynamics. In contrast to the viability estimation, this is an approximation from the interior.

Improvements and extensions to this algorithm are currently under intensive research. Relations to dynamical programming Frankowska1989 (51) and other extensions Bayen2002 (52, 48) can provide the minimal time to reach the target set. Also, one can even find controllers that drive the system to the target set Cardaliaguet1998 (53, 54, 55).

4.2 Dealing with an unbounded state space

There are multiple ways to map an unbounded state space to a bounded one, depending on the specific need for the system. In case of the example system analyzed later, each coordinate is bounded from below and unbounded from above. This is rather common in socioeconomic models, in particular due to unbounded economic growth. Hence, we propose a solution that maps each coordinate separately.

We assume a general dynamical system given by a set of ordinary differential equations

[TABLE]

with ${\mathcalboondox{x}}\in\mathbb{R}_{\geq 0}^{n}$. In contrast to Eq. 1, there is no control parameter, because the dependence on the control parameter is irrelevant here and the mapping can be done for general ordinary differential equations. We imply an extension to control systems by considering fixed controls. Then we propose the coordinate transformation

[TABLE]

where ${\mathcalboondox{x}}_{i,mid}\in\mathbb{R}_{\geq 0}$ are parameters. Applying this transformation on the dynamics leads to a new set of ordinary differential equations

[TABLE]

for ${\mathcalboondox{y}}\in[0,1)^{n}$, i.e. inside a bounded space, where $D\Phi$ is the Jacobian of $\Phi$ and $\circ$ is the symbol for function composition. The parameters $x_{i,mid}$, summarized to the vector $x_{mid}$, define the scales for each coordinate that is resolved best when discretizing the state space for numerical estimations because $\Phi(x_{mid})=(\frac{1}{2},\frac{1}{2},\dots)^{T}$. So they should be taken to be around the main region of interest.

For the dynamical system including the management options, we assume the control parameters to be fixed and let the system switch between the default dynamics and the management options instantaneously. Therefore, the transformation given in this section can also be applied to the controlled system.

4.3 Nonlinear, local time-homogenization

A problem during the estimation of viability kernels from Section 2 is a possibly inhomogeneous time scale, i.e. that the (norm of the) RHS function of the control system (1) can have values through several orders of magnitude. For instance, models like Eq. 27 often lead to divergences at the upper boundary of a coordinate.

This problem can be addressed by rescaling the time of the system in a nonlinear way. The used definitions of viability theory depend only qualitatively but not quantitatively on time: It does not matter for the definition of viability kernels and capture basins how long the system takes to reach the target set or desirable region, it only matters that it is possible to reach in a finite time period. Hence, viability kernels and capture basins in the rescaled system are equivalent to the original ones. Because there is no functional difference between these quantities in the original and the homogenized system, it is technically possible to use the latter in order to estimate them in the former.

As the control parameter is not necessary for the rescaling of the system, we use a differential equation

[TABLE]

We propose to use the new system

[TABLE]

Assuming $\epsilon$ is small enough, this new system generally fulfills three criteria:

The systems (28) and (29) are orbitally equivalent (cf. (Kuznetsov1998, 56, p. 42), Definition 2.4), i.e. the trajectories of solutions with the same initial conditions follow the same path. In other words, only the time has been rescaled. 2. 2.

Everywhere away from fixed points $\|F({\mathcalboondox{y}})\|\gg\epsilon$ holds and hence the time scale is properly homogenized

[TABLE] 3. 3.

Within a small enough environment of a fixed point, $\|F({\mathcalboondox{y}})\|\ll\epsilon$ holds. We can therefore approximate

[TABLE]

This implies that the function $\tilde{F}$ goes to zero at fixed points of the original system with the same properties as $f$ at these points (e.g. $\mu$-Lipschitz or $C^{\infty}$, same Lyapunov-Exponents etc.).

Because the units of the coordinates of ${\mathcalboondox{y}}$ might be different from each other, there is no real physical interpretation of $\tilde{F}$. But that is not necessary either as it is only an auxiliary system for the estimation with the Saint-Pierre Algorithm.

4.4 Estimation of implicitly defined capture basins: eddies

For capture basins we only care about entering a target set at least once. However, eddies are defined by being able to reach the desirable region again and again, so the most natural definition is an implicit one as in Eqs. 16a and 16b. In order to estimate them, we find an alternate definition in terms of a limit process.

We start by defining the largest sets that could contain eddies

[TABLE]

and then use the iteration step

[TABLE]

for $i=1,2,\dots$. Note that ${{\mathcalboondox{E}}^{-}}_{i}$ used in Eq. 33b is already computed in Eq. 33a. So this can really be seen as a step-by-step prescription. Thus, the eddies can be recovered as

[TABLE]

The limit exists because both sequences are monotone and nonincreasing. The Saint-Pierre algorithm works on a discretized state space with finitely many elements. Hence, there exists an $k\in\mathbb{N}$ such that ${{\mathcalboondox{E}}^{+}}={{\mathcalboondox{E}}^{+}}_{k}={{\mathcalboondox{E}}^{+}}_{k-1}$ and ${{\mathcalboondox{E}}^{-}}={{\mathcalboondox{E}}^{-}}_{k}={{\mathcalboondox{E}}^{-}}_{k-1}$. Thus, the existence of the limit results in the convergence of the algorithm after a finite number of steps.

This iteration process follows the idea of being able to visit the desirable region again and again and is an algorithmic description for the estimation of eddies. Other similarly implicitly given sets can be estimated by adjusting this basic idea.

5 Example: the AYS low-complexity model of climate change, wealth and energy transition

Here, we demonstrate the applicability of our operationalization of the tsm-framework using a three dimensional example model. We start with describing the model’s default dynamics, the management options and the desirable region. Then, we analyze the model’s attractors and estimate the current state within the state space. Finally, we apply the coordinate transformations and discuss the resulting tsm regions found with the Saint-Pierre algorithm.

5.1 Model default dynamics

To develop our low-complexity example model, we took inspiration from Kellie-Smith2011a (43) and added a renewable energy sector with a learning-by-doing dynamics. The model describes the interactions between the atmosphere and a stylized economy whose energy supply comes from either fossil fuels and renewable resources. While the atmosphere reacts to CO2 emissions by increasing temperatures, the higher temperatures have negative impact on economic output via climate damages. The share between the two energy sources is determined by their relative cost, which depends on previous use of renewables and optional policies. The model structure is depicted in Fig. 3. It comprises three state variables.

1. The first variable is the excess atmospheric carbon stock $A$ [GtC = giga tons of carbon], measured w.r.t. a pre-industrial level $A_{0}\approx 600\,$GtC. It increases with current CO2 emissions $E$ [GtC/a = GtC per year]. Taking $A_{0}$ as an estimate for the long-term no-emissions equilibrium value, we assume $A$ approaches zero if $E=0$, due to carbon uptake by oceans, plants and soil. To keep the complexity of the model as low as possible, we do not explicitly model a carbon cycle as in Anderies2013 (8) but simply assume the carbon uptake leads to an exponential relaxation towards equilibrium on a characteristic time scale of $\tau_{A}\approx 50\,$a [a = years]. Hence our first model equation is

[TABLE]

where $E$ will be derived below from economic assumptions.

2. The second variable is economic output or production $Y$ [US$/a] of the global economy. It thus represents the gross world product as a common indicator. We assume the economy to have a positive basic growth rate$\beta\approx 3,%[1/a] and additional climate impacts as in Kellie-Smith2011a ([43](#bib.bib43)). As a proxy for temperature we simply useA$, effectively assuming an infinitely fast greenhouse effect. Hence this terms is represented by$-\theta AY$where$\theta\approx 8.57\cdot 10^{-5},/$(GtC a) is a temperature sensitivity parameter chosen such that the total growth rate$\beta-\theta A$becomes negative when$A$exceeds the level corresponding to a global warming of$+2,^{\circ}\text{C}$. This gives

[TABLE]

3. The third state variable is the renewable energy knowledge stock $S$ that indicates how much knowledge is available for the production of renewable energy $R$ [GJ/a = giga joule per year]. In accordance with Wright’s law (e.g., Nagy2013 (57)) of learning-by-doing, we basically identify $S$ with the past cumulative production of renewables and thus measure it in units of [GJ]. To account for the human capital component, we additionally assume that knowledge depreciates on a characteristic time scale of $\tau_{S}\approx 50\,$a. We choose this value to be in the upper range of a typical professional life span. Cumulation and depreciation then give

[TABLE]

where $R$ will be derived below.

Finally, to determine $E$ and $R$, we use the following simplistic economic assumptions. The energy demand $U$ [GJ/a] is proportional to the economic output

[TABLE]

where $\epsilon\approx 147\,$US$/GJ is an energy efficiency parameter. This demand is satisfied by a mix of fossil and renewable energy which are assumed to be perfect substitutes (and ignoring other energy sources such as bioenergy). Their respective shares are determined by a price equilibrium. We assume convex monomial cost functions and unit costs of renewable energy that show a power-law decay with growing$S$ Nagy2013 (57). This implies that the fossil sector has a share given by the sigmoidal function

[TABLE]

where $\sigma\approx 4\cdot 10^{12}\,GJ$ is the break-even knowledge level at which renewable and fossil extraction costs become equal, and $\rho\approx 2$ is a dimensionless parameter determined from the cost convexity and learning rate. $\Gamma$ approaches unity (no renewables) as $S\to 0$ and zero (no fossils) as $S\to\infty$. Fossil and renewable energy flows and emissions are then

[TABLE]

where our final parameter $\phi\approx 4.7\cdot 10^{10}\,$GJ/GtC is the fossil fuel combustion efficiency. This completes the equations of the AYS model.

The three dynamical variables $A$, $Y$ and $S$ are interrelated due to the various connecting equations and nonlinearities arise particularly due to Eq. 39 and Eq. 36. The model’s default dynamics is also specified by the eight parameters $\tau_{A}$, $\beta$, $\theta$, $\tau_{S}$, $\epsilon$, $\rho$, $\sigma$ and $\phi$. Appendix C contains details on how we estimated these parameters. The resulting flow is depicted in Fig. 4, where the basins of the system’s two attractors (discussed in Section 5.5) are already colored differently.

5.2 Management options (controls)

The above parameter values define what we consider the default dynamics because they represent a business-as-usual case. This means humanity applies no specific management that would alter the way things usually go.

In addition to the default dynamics, we consider two distinct management options. They represent possible policy choices that may be combined in any way, leading to more or less shifted trajectories. Because the management options are defined as changes in the parameters $\beta$ and $\sigma$, the control parameter $u$ of the model is two-dimensional and the particular control $u_{0}$ for the default flow is

[TABLE]

1. The first distinct management option we consider is low growth (lg). This corresponds to policies that voluntarily reduce the basic growth rate $\beta$ to half its value $\beta_{\textsc{lg}}=1.5\%/$a, leading to the control

[TABLE]

2. The second management option we introduce aims at mitigating climate change by inducing an energy transition (et), e.g. via taxing fossils and/or subsidizing renewable energy use. These policy instruments shift the relative costs of fossil and renewable energy, which according to Eq. 39 can be effected in our model by a reduction of $\sigma$. Hence, we represent this option by reducing $\sigma$ to approx. $(1/2)^{1/\rho}=1/\sqrt{2}$ of its default value, i.e. to $\sigma_{\textsc{et}}=2.83\cdot 10^{12}$ GJ, corresponding to dividing the renewable to fossil cost ratio by half and leading to the control

[TABLE]

With these two management options, the set of possible managements is given by ${\mathcal{U}}_{m}=\left\{u_{\textsc{lg}},u_{\textsc{et}}\right\}$.

5.3 Desirable states

Steffen et al. Steffen2015 (2) set the planetary boundary for climate change (pb-cc) to $350\,$ppm [ppm = parts per million] with an uncertainty zone until $450\,$ppm. We chose the desirable region to be where $A>A_{PB}=345\,$GtC (above the pre-industrial level $A_{0}$), corresponding to the looser 450ppm boundary (see Appendix C for the conversion).

The safe and just operating space by Raworth2012 (3) complements the planetary boundaries with social foundations. Here we use a lower bound on the world gross product for illustration of a social foundation, even though we are well aware of the serious shortcomings of this indicator to measure decent livelihoods and basic human need satisfaction (e.g. Costanza2014 (58, 59, 60)). Most importantly, the aggregate product does not capture any distributional aspects between countries and individuals. The exact value for such a lower boundary is open for discussion. For pragmatic reasons, we use here $Y_{SF}=4\cdot 10^{13}\,$US$ (sf-y), the economic production of the year 2000.

Fig. 5 shows the location of the boundaries in the three dimensional state space as two shaded planes. The desirable region is defined by the two boundaries and located on the left in all panels of Fig. 5.

5.4 Current state

The current state ${\mathcalboondox{x}}_{c}=\left(A_{c},Y_{c},S_{c}\right)$ can be estimated as follows. $A_{c}$ is currently around 240 GtC, corresponding to a concentration of 400 ppm Betts2016 (61), and the world gross product of 2015 is around 70 Trillion USWorldbankGDP ([62](#bib.bib62)).S_{c}$is estimated on the basis of the total past renewable energy consumption of roughly$2\cdot 10^{12} GJ Mongabay2015 ([63](#bib.bib63)). This number has accumulated over roughly the same time as the characteristic knowledge depreciation time, \tau_{S}=50$a. We assume roughly half of it has already depreciated, leaving$10^{12}$GJ. Further accounting for the large error margins involved in estimating this number and because it contains hydroelectricity, whose growth potential is strongly limited, we take a conservative estimate by setting$S_{c}$to a quarter of this values,$5\cdot 10^{11}$ GJ.

[TABLE]

This estimation is quite rough and only provides limited information about the state of the target system because it reduces it to the three dimensions of a simplified model. However, it shows that the current state is still well inside the desirable region.

5.5 Attractors

With the above parameter values, the default dynamics of the AYS model has two fixed points. The first fixed point has the coordinates at

[TABLE]

It represents a fossil-based economy without renewable energy use, reduced economic output and constant climate damages. Therefore, we call it the “black fixed point” ${\mathcalboondox{x}}_{b}$.

The second fixed point is located at the boundaries of the state space,

[TABLE]

with “$+\infty$” describing the boundary of the state space, as becomes clear in Section 5.6. This attractor corresponds to an economy with unbounded exponential growth of economic output and renewable knowledge. Also, fossil fuel usage and emissions decline towards zero. Therefore, we call this attractor the “green fixed point”.

Both the black and green fixed points represent extreme asymptotic states. Even if they might correspond to unrealistic special cases of the modeled system, we find this acceptable because our focus lies on the model transients and the model only serves as an example. With regard to the desirable regions, our analysis of the default dynamics shows that the green fixed point ${\mathcalboondox{x}}_{g}$ is not violating either of the two boundaries, while the black one ${\mathcalboondox{x}}_{b}$ violates the pb-cc as can be seen in Fig. 5a.

Applying the low growth management option moves the black fixed point to ${\mathcalboondox{x}}_{b,\textsc{lg}}=(A,Y,S)=(175\,$GtC$,2.42\cdot 10^{13}\,$US$/a$,0), no longer violating the pb-cc (see [Fig. 5b](#S5.F5.sf2)) but now violating the sf-y. Furthermore, for our choice of parameters and boundaries, no choice of the basic growth rate \beta$ can make the black fixed point lie within the desirable region, thus satisfying both the pb-cc- and sf-y-boundaries at the same time.

Applying the energy transition management option does not affect the location of the two attractors. But, more importantly, it changes the shape of the basins of attractions. When carefully inspecting Fig. 5c, one can see that the volume the basin of attraction of the green fixed point is enlarged in comparison to the default flow in Fig. 5a. Within the concept of Basin Stability Menck2013 (64, 65) the volume of the basin of attraction has been found to be an important indicator for an attractor’s stability hence we will use a similar approach for the bifurcation analysis in Section 5.8.

5.6 Dealing with the unbounded state space

The model described in Section 5 and given by Eqs. (35)-(40) is defined on an unbounded state space. In this section, we map it to a bounded state space in order to deal with the diverging dynamics and apply the Saint-Pierre algorithm. We map the original state space ${\mathcalboondox{X}}=\mathbb{R}_{\geq 0}^{3}$ with the dynamical variables ${\mathcalboondox{x}}=(A,Y,S)$ to a bounded space ${\mathcalboondox{Y}}=[0,1)^{3}$ with transformed variables ${\mathcalboondox{y}}=(a,y,s)$ and then add the point ${\mathcalboondox{y}}_{g}=(0,1,1)$ which is the equivalent of ${\mathcalboondox{x}}_{g}$ in the new coordinates. So we perform a change of coordinates

[TABLE]

where we switch from the old right-hand side (RHS) $f$ to the new RHS $F$. Following the explanations in Section 4.2, we use the transformation

[TABLE]

where the parameters ${\mathcalboondox{x}}_{mid}=(A_{mid},W_{mid},S_{mid})^{T}$ are such that $\Phi({\mathcalboondox{x}}_{mid})=(\frac{1}{2},\frac{1}{2},\frac{1}{2})^{T}$. They can be understood as the scale where the “resolution is the best”. Hence, changing this value does not qualitatively influence the result, but a good choice can make them clearer.

This is exactly the transformation that has been used to create Figs. 4, 5a, 5b and 5c and will be used for all the following figures, too. As we care most about the current state of the world, we choose ${\mathcalboondox{x}}_{mid}={\mathcalboondox{x}}_{c}$.

Using Eq. 27, we get a new set of ODEs with ${\mathcalboondox{y}}=(a,y,s)$ as coordinates

[TABLE]

where $\gamma$ is the equivalent of $\Gamma$ in Eq. 39 but in the ${\mathcalboondox{y}}$-coordinates. The fixed points in the new $y$-coordinates are

[TABLE]

Now, we formally extend the dynamics such that $F({\mathcalboondox{y}}_{g})=0$.

5.7 Results

We analyze the model in its compactified form Eqs. 59a, 59b, 59c and 59d using the nonlinear, local time-homogenization (Section 4.3) and the Saint-Pierre algorithm (Section 4.1). Here, we do not provide the equations for the time-homogenized version because they are lengthy, their calculation straight-forward and they do not give much insight. The most important identified regions are depicted in Fig. 6. We discuss their properties in the following.

The first regions to note are the shelter and the backwater shown in Fig. 6a. In our model, the shelter, in which the system can stay in the desirable region without management, is the invariant kernel of the green fixed point’s basin of attraction when restricting it to desirable states ${{\mathcalboondox{X}}^{+}}$ only. States in this region correspond to a system that has accumulated enough knowledge for the energy production with renewable resources such that they become so cheap that there is no need for fossil fuels anymore. So the remaining (excess) CO2 (above long-term equilibrium) is removed over time due to natural carbon uptake, leading the system to the green fixed point. The model also contains a glade, i.e., a region from which one can reach the shelter through the desirable region. It is a thin layer under the shelter and has not been visualized in Fig. 6.

The backwater, where one can stay in the desirable region forever but needs to apply management over and over again, is the partition of the desirable region where the growth of economic output and hence of emissions can be restricted. That way, the atmospheric carbon concentration can be kept within the planetary boundary. Within the backwater, decarbonization of the economy is impossible because the given maximal carbon tax or renewable subsidy policies are too weak to make renewables competitive with fossil fuel. Instead, one can use the low growth management option to stay in a state that corresponds to a carbon-based economy, in which climate impacts are compensated by low economic growth. The atmospheric carbon level is relatively high but still within the boundary and in equilibrium with the emissions. Using the low growth option, one can therefore stay within the desirable region but cannot reach the green fixed point.

For simplicity, we did not include the option to choose a value of the base growth rate lying between the two options $\beta$ and $\beta_{\textsc{lg}}$. So formally, the management strategy required to stay in the desirable region involves switching between $\beta$ and $\beta_{\textsc{lg}}$ because either of these two values alone leads to a black fixed point in the undesirable region. Still, it is easy to see that this management strategy is equivalent to using a constant intermediate value of $\beta_{m}$ (e.g. $\beta_{m}=2.7\,\%/$a) instead. The dynamics with such an intermediate growth rate has a black fixed point that lies in the desirable region. So this option is implicitly included in the model. In other words, one only needs to model the maximal and minimal values of the option space explicitly because the tsm-framework allows for arbitrarily fast switches between the management options. Intermediate options, more precisely all convex combinations of options, are implicitly taken into account. Still, replacing the value of $\beta_{\textsc{lg}}$ with $\beta_{m}$ would introduce further changes. The maximal management in the transient would be restricted and thus the size of the backwaters reduced. We show this in Section 5.8. Here, it becomes obvious that the tsm-partition depends both on the asymptotics and the transient dynamics of the model.

According to our analysis, the current state ${\mathcalboondox{x}}_{c}$ (as estimated in Section 5.4) lies between the shelter and the backwater discussed above, in a time-limited lake as the zoom into the state space in Fig. 6b shows. Therefore, the results suggest that humanity currently faces a so-called lake dilemma in the AYS model. In such a dilemma, we have to make the qualitative decision between staying within the desirable region uninterruptedly but depending on management forever or going through the undesirable region to finally end up in the shelter. In the model, the choice is between using the energy transition option in order to speed up the knowledge accumulation on renewable resource use and finally reach the green fixed point or the low growth option that restricts the total energy demand and helps reaching the black fixed point. Even with a combined usage of the different options, the model does not allow to reach the shelter without transgressing the boundaries when going for the first option. This classification of ${\mathcalboondox{x}}_{c}$ depends on the finite resolution $h$ of the applied algorithm. An alternative analysis of the AYS model with deep reinforcement learning suggest that ${\mathcalboondox{x}}_{c}$ is in a glade region instead of a lake region strnad2020 (46).

5.8 Bifurcation analysis

Varying the parameters corresponding to the two management options may lead to what can be understood as tsm* bifurcations* with respect to the topology of tsm regions in state space. A full theory of tsm bifurcations would be based on a suitable notion of topological equivalence between control systems with default trajectories and specified desirable regions, just like the standard theory of bifurcations is based on such an equivalence between (uncontrolled) dynamical systems Kuznetsov1998 (56). Because this is beyond the scope of this paper, we use a working definition of tsm bifurcation here that is based on sudden, qualitative changes of the tsm partition at certain critical parameters (see Appendix B for details).

This allows us to use the relative volume of each region an indicator for tsm bifurcations, motivated by the concept of Basin Stability Menck2013 (64, 65) and its extensions Hellmann2015 (66, 67, 68, 69, 70, 71, 72). We use uniformly distributed points in state space for the Saint-Pierre algorithm. Therefore, we can estimate the relative volume of one region with the number of points associated to this region over the total number.

When varying $\beta_{\textsc{lg}}$ corresponding to the low growth management option (lg) from $1.5\,\%/$a to $3.5\,\%/$a, a downstream-eddies bifurcation occurs. Until the fixed point of the lg flow crosses the planetary boundary at the critical value of $\beta_{\textsc{lg}}=$ $\beta_{PB}\approx$ $2.95\%/$a there exists a backwater333Because of the discretization of state space and time in the estimation, Fig. 7a suggests that the bifurcation occurs already at $\approx 2.8\%$/a. However, this is a numerical inaccuracy.. Beyond this point, there are only eddies left because the focus of the default flow and the one of the low-growth flow are both in the undesirable region. Thus, for $\beta_{\textsc{lg}}>\beta_{PB}$ it is only possible to switch between the two flows such that the system circles around both foci and visits the desirable region again and again, while having to pass through the undesirable region in between.

As discussed in Section 5.7, the backwater can be explained by the fact that there is an in-between value $\beta_{m}=2.7\,\%/$a, for which the black fixed point lies in the desirable region. Fig. 7a shows that the volume of the backwater for this intermediate value $\beta_{m}=2.7\,\%/$a is smaller than for $\beta_{\textsc{lg}}=1.5\,\%/$a. This is due to the restriction of the maximally available management such that it becomes impossible to stay within the desirable region forever for a large share of trajectories starting of in the downstream.

Fig. 7b depicts the changes in relative volume of the tsm-partition over the energy transition control $\sigma_{ET}$. Especially for a strong management corresponding to small values of $\sigma_{ET}$, we observe strong changes in the relative volume of the glade, the lake and the backwater. The glade increases in size, because it is possible to reach the green fixed point from more initial conditions without transgressing the boundaries.

6 Summary & Outlook

In this paper, we introduce a variant definition of the topology of sustainable management (tsm) in terms of viability theory that operationalizes the partition of state space into qualitatively different regions. We apply the Saint-Pierre algorithm to estimate the volume and visualize the position and shape of different tsm regions. Because the algorithm works on bounded spaces only, we use a coordinate transformation to a bounded space that best resolves a fixed scale of interest. Furthermore, we address the problem that time scales vary by orders of magnitude by introducing a nonlinear, local time homogenization.

We apply these novel concepts to a three dimensional example model of interactions between climate change, economic growth and energy transition. Even though we design the model in a stylized way, the analysis reveals a rich tsm topology. This results from the two fixed points in the model corresponding to a fossil and a green economy. We estimate the current state of the world to lie inside a finite-time lake. In the model, humanity is therefore facing a pressing lake-dilemma: it has limited time to make a choice between two qualitatively different options. Moreover, we analyze how the relative volume of different tsm regions change with the management parameters. We identify a downstream-eddies bifurcation, in which a downstream region becomes an eddie after one of the fixed points crosses the boundary of the desirable region.

The example model presented here has a comparatively low-complexity that allows us to focus on the application of the novel methods for estimating the tsm regions. Relevant processes, like the carbon cycle, economic dynamics and different energy productions, are represented only in a stylized way. Future work can apply our framework to extended versions of our model. For example, it could incorporate other management options, such as carbon capture and storage, or introduce higher dimensional representations of the carbon cycle or economic system. The models developed in Refs. Nitzbon2017 (45, 73) could to be good candidates for such extensions because they have only a limited number of dynamical variables.

Extending the framework also to high-dimensional models would require improved algorithms for the estimation of viability kernels and capture basins. While there is work on efficient algorithms for such problems Bokanowski2006 (39, 38, 74, 75), they need to be adjusted and extended to be applicable for the computation of the tsm partition. Furthermore, future theoretical research could further formalize our working definition for tsm bifurcations and extend the tsm framework to stochastic dynamical systems Anishchenko2007 (76). Such further developments of the presented methods would be highly valuable to advance our understanding of complex social-ecological models of the Earth system Mueller-Hansen2017 (77, 78) that can guide sustainable development within a safe and just operating space.

Software availability

For this project, an open source Python-library for viability computations called pyviability has been developed and is available online at https://timkittel.github.io/PyViability/ (published July 7, 2017). The source code for the model-specific computations including a description to reproduce the results presented in this article is available at https://github.com/timkittel/ays-model/ (published July 7, 2017).

Acknowledgements.

Acknowledgement

This article was developed within the scope of the IRTG 1740/TRP 2011/50151-0, funded by the DFG/FAPESP. This work was conducted in the framework of PIK’s flagship project on coevolutionary pathways (copan). The authors thank CoNDyNet (FKZ 03SF0472A) for their cooperation. The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research. The authors thank the developers of the used software: Python python (79), Numerical Python numpy (80) and Scientific Python scipy (81).

The authors thank Sabine Auer, Wolfram Barfuss, Karsten Bölts, Catrin Ciemer, Eduardo Costa, Jonathan Donges, Reik Donner, Jaap Eldering, Jasper Franke, Roberto Gueleri, Frank Hellmann, Jakob Kolb, Julien Korinman, Till Koster, Chiranjit Mitra, Jan Nitzbon, Ilona Otto, Tiago Pereira, Camille Poignard, Francisco Rodrigues, Edmilson Roque dos Santos, Stefan Ruschel, Paul Schultz, Fabiano Berardo de Sousa, Lyubov Tupikina, and Kilian Zimmerer for helpful discussions and comments.

Appendix

Appendix A Existence of eddies

Definition of eddie-like sets:

We call a pair of sets $\mathcalboondox{A}^{+/-}\subseteq{\mathcalboondox{X}}$ eddie-like if and only if they fulfill the following two conditions: (i) $\mathcalboondox{A}^{+/-}\subseteq{{\mathcalboondox{X}}^{+/-}}-{\mathcalboondox{U}}-{\mathcalboondox{D}}$ and (ii) $\mathcalboondox{A}^{+/-}\subseteq\text{Capt}\left(A^{-/+}\right)$. Note the inverted order of the signs in the last term.

Union of two eddie-like pair sets are also eddie-like:

For two eddie-like pairs of sets ${{\mathcalboondox{E}}^{+/-}}_{1}$ and ${{\mathcalboondox{E}}^{+/-}}_{2}$, the union pair ${{\mathcalboondox{E}}^{+/-}}_{3}={{\mathcalboondox{E}}^{+/-}}_{1}\cup{{\mathcalboondox{E}}^{+/-}}_{2}$ is eddie-like, too.

Proof:

The first condition is trivially fulfilled and the second one follows straight away from $\text{Capt}\left(\mathcalboondox{A}\right)\cup\text{Capt}\left(\mathcalboondox{B}\right)=\text{Capt}\left(\mathcalboondox{A}\cup\mathcalboondox{B}\right)$ for two state sets $\mathcalboondox{A},\mathcalboondox{B}\subseteq{\mathcalboondox{X}}$. Hence, the union of all eddie-like pairs of sets is maximal and eddies exist. ∎

Appendix B Working definition of tsm bifurcations

For the analysis in Section 5.8, we use the following working definition of tsm bifurcations. Let $\psi=(\psi_{1},\dots,\psi_{k})$ be a vector of parameters of the control system, let ${\mathcalboondox{S}}(\psi),{\mathcalboondox{G}}(\psi),\dots,{\Theta}(\psi)\subseteq{\mathcalboondox{X}}$ be the various tsm regions in dependence on parameter values, and let $f(\psi)=({\mathcalboondox{S}}(\psi),{\mathcalboondox{G}}(\psi),\dots,{\Theta}(\psi))$ be the function that maps parameter values to the tuple of these regions. Now, we say that at a certain point $\psi^{\ast}$ in parameter space a tsm* bifurcation* happens if and only if $f(\psi)$ is discontinuous at $\psi^{\ast}$ with respect to the Hausdorff metric between subsets of ${\mathcalboondox{X}}$. In particular, we have a tsm bifurcation if the volume of some of the regions changes discontinuously at $\psi^{\ast}$.

Appendix C Parameter estimation

To get a roughly realistic setting, we estimated the parameters of the model using several publicly available data sources.

$A_{0}$ was taken from Ciais2013 (82) and slightly rounded. $\tau_{A}$ and $\beta$ were taken from Kellie-Smith2011a (43). $\phi$ was based on the ton oil equivalent of various fossil fuels and a typical mass share of 90% carbon in fossil fuels, as described in Nitzbon2016thesis (73).

Assuming that two degrees warming correspond to a carbon concentration of around 450 ppm IPCC2013 (83) and thus to a carbon stock of 950 GtC (both being 1.6 times their pre-industrial value), we require that the total growth rate $\beta-\theta A_{1}$ becomes zero for $A_{1}\approx 950\,$GtC${}-A_{0}=350\,$GtC, hence $\theta$ was taken to be $\beta/A_{1}$ $\approx$ $8.57\cdot 10^{-5}/{}$(GtC a).

$\epsilon$ was estimated from the World Bank’s primary energy intensity data WorldbankEnergy (84).

For $\tau_{S}$, the characteristic depreciation time of renewable energy knowledge, no reliable source was found, so we made a very coarse guess by setting it roughly to the length of an average working life of 50 a.

The break-even knowledge level $\sigma$ was also estimated very coarsely. According to Mongabay2015 (63), past cumulative world consumption of renewable energy is $\approx 2\cdot 10^{18}\,$Btu $\approx 2\cdot 10^{12}\,$GJ or roughly 20 years of world energy consumption. To be on the conservative side and avoid overestimating the potential of renewables, we took $\sigma$ to be two times that value.

$\rho$ was set as follows. We assume fossil and renewable energy production costs of $C_{F}\propto F^{1+\gamma}$ and $C_{R}\propto R^{1+\gamma}/S^{\lambda}$, where $\gamma>0$ is a convexity parameter and $\lambda>0$ is a learning exponent. Then energy prices are $\pi_{F}\propto\partial C_{F}/\partial F$ $\propto F^{\gamma}$ and $\pi_{R}\propto\partial C_{R}/\partial R$ $\propto R^{\gamma}/S^{\lambda}$. In the price equilibrium, $\pi_{F}=\pi_{R}$, hence $R/F\propto S^{\lambda/\gamma}$, and thus $\rho=\lambda/\gamma$. According to Rubin2015 (85), the learning rate $LR=1-2^{-\lambda}$ of several renewables is around 1/8, hence $\lambda\approx\log_{2}(8/7)\approx 0.2$. Assuming a mild convexity of $\gamma\approx 0.1$, we get $\rho\approx 2$.