**Contents**


## **About the Special Issue Editor**

**Ioannis Dassios** is currently a UCD Research Fellow/Assistant Professor at AMPSAS, University College Dublin, Ireland. His research interests include dynamical systems, mathematics of networks, differential and difference equations, singular systems, fractional calculus, optimization methods, linear algebra, and mathematical modeling (materials, electrical power systems, economic models, etc). He studied Mathematics, completed a two-year MSc in Applied Mathematics and Numerical Analysis, and obtained his Ph.D. degree in Applied Mathematics at University of Athens, Greece with the grade "Excellent" (highest mark in the Greek system). He had positions at the University of Edinburgh, U.K; University of Manchester, U.K.; and University of Limerick, Ireland. He has published more than 55 articles in internationally leading academic journals and has participated in several international collaborations. He has served as a reviewer more than 500 times in more than 75 different journals, he has been member in scientific and organizing committees of international conferences, and he is also member of editorial boards of peer-reviewed journals. Finally, he has received several awards such as travel grants, for reviewing and for his contributions to his institute.

## *Editorial* **Special Issue on Mathematical Modeling Using Differential Equations and Network Theory**

#### **Ioannis Dassios**

AMPSAS, University College Dublin, Dublin 4, Ireland; ioannis.dassios@ucd.ie

Received: 1 March 2020; Accepted: 3 March 2020; Published: 10 March 2020

#### **1. Introduction**

This special issue collects the latest results on differential/difference equations, the mathematics of networks, and their applications to engineering, and physical phenomena. The Special Issue has 42 submissions and eight high-quality papers which got published with original research results. The Special Issue brought together mathematicians with physicists, engineers, as well as other scientists. Topics covered in this issue:


#### **2. Acoustic Wave Equations Using Fractional-Order Differential Equations**

In [1], the authors present a newly developed technique, defined as a variational homotopy perturbation transform method in order to solve fractional-order acoustic wave equations. The basic idea behind this article is to extend the variational homotopy perturbation method to the variational homotopy perturbation transform method.

The proposed method is an accurate and straightforward technique to solve fractional-order partial differential equations, and can be considered as a practical analytical technique to solve non-linear fractional partial differential equations compared to other analytical techniques existing in the literature. Several illustrative examples verify the method.

#### **3. Analytical Solutions of Dimensional Physical Models Using Modified Decomposition Method**

In [2], the authors present a new analytical technique based on an innovative transformation in order to solve (2+time fractional-order) dimensional physical models. The proposed method is based on the hybrid methodology of Shehu transformation along with the Adomian decomposition method.

The solutions of the targeted problems are represented by graphs and are obtained in a series form that has the desired rate of convergence. The method is, in general, a practical analytical technique to solve linear and non-linear fractional partial differential equations. Numerical examples are given using the proposed method.

#### **4. Multi-Switching Combination Synchronization of Fractional-Order Delayed Systems**

In [3] the authors investigate multi-switching combination synchronization of three fractional-order delayed systems. This is actually a generalization of previous multi-switching combination synchronization of fractional-order systems by introducing time-delays.

Based on the stability theory of linear fractional-order systems with multiple time-delays, the article provides appropriate controllers to obtain multi-switching combination synchronization of

three non-identical fractional-order delayed systems. In addition, numerical simulations show that they are in accordance with the theoretical analysis given.

#### **5. An Overview of Early Developments of the Hardy–Cross-Type Methods**

In [4], the authors provide an overview of early developments of the Hardy–Cross-type methods for computation of flow distribution in pipe networks.

Cross originally proposed a method for analysis of flow in networks of conduits or conductors in 1936. His method was the first really useful engineering method in the field of pipe network calculation. Only electrical analogs of hydraulic networks were used before the Hardy–Cross method. A problem with flow resistance versus electrical resistance makes these electrical analog methods obsolete. The method by Hardy–Cross is taught extensively at faculties, and it remains an important tool for the analysis of looped pipe systems. Engineers today mostly use a modified Hardy–Cross method that considers the whole looped network of pipes simultaneously (use of these methods without computers is practically impossible).

In addition, in [4] a method from a Russian practice published during the 1930s, which is similar to the Hardy–Cross method, is also described. Some notes from the work of Hardy–Cross are also presented. Furthermore, an improved version of the Hardy–Cross method, which significantly reduces the number of iterations, is presented and discussed.

Finally, the authors present results on tested multi-point iterative methods, which can be used as a substitution for the Newton–Raphson approach used by Hardy–Cross.

#### **6. Parametrical Non-Complex Tests to Evaluate Partial Decentralized Linear-Output Feedback Control Stabilization Conditions**

In [5], the authors formulate sufficiency-type linear-output feedback decentralized closed-loop stabilization conditions if the continuous-time linear dynamic system can be stabilized under linear output-feedback centralized stabilization.

The provided tests are simple to evaluate, while they are based on the quantification of the sufficient smallness of the parametrical error norms between the control, output, interconnection and open-loop system dynamics matrices and the corresponding control gains in the decentralized case related to the centralized counterpart.

The tolerance amounts of the various parametrical matrix errors are described by the maximum allowed tolerance upper-bound of a small positive real parameter that upper-bounds the various parametrical error norms. Such a tolerance is quantified by considering the first or second powers of such a small parameter.

The results are seen to be directly extendable to quantify the allowed parametrical errors that guarantee the closed-loop linear output-feedback stabilization of a current system related to its nominal counterpart. Several numerical examples are included and discussed in the article.

#### **7. Transient-Flow Modeling of Vertical Fractured Wells with Multiple Hydraulic Fractures**

Massive hydraulic fracturing of vertical wells has been extensively employed in the development of low-permeability gas reservoirs. The existence of multiple hydraulic fractures along a vertical well makes the pressure profile around the vertical well complex.

In [6], the authors study the pressure dependence of permeability in order to develop a seepage model of vertically fractured wells with multiple hydraulic fractures. Both transformed pseudo-pressure and perturbation techniques have been employed to linearize the proposed model.

The proposed work further enriches the understanding of the influence of the stress sensitivity on the performance of a vertical fractured well with multiple hydraulic fractures and can be used to more accurately interpret and forecast the transient pressure.

Some key points in the article are the superposition principle and a hybrid analytical-numerical method that are used to obtain the bottom-hole pseudo-pressure solution, the type curves for pseudo-pressure that are presented and identified, and finally, the discussion that is included on the effects of the relevant parameters on the type curve and the error caused by neglecting the stress sensitivity.

#### **8. Policy-Compliant Maximum Network Flows**

Computer network administrators are often interested in the maximal bandwidth that can be achieved between two nodes in the network, or how many edges can fail before the network gets disconnected. Classic maximum flow algorithms that solve these problems are well-known. However, in practice, network policies are in effect, severely restricting the flow that can actually be set up. These policies are put into place to conform to service level agreements and optimize network throughput, and can have a large impact on the actual routing of the flows.

In [7], the authors model the problem and define a series of progressively more complex conditions and algorithms that calculate increasingly tighter bounds on the policy-compliant maximum flow using regular expressions and finite-state automata. This is the first time that specific conditions are deduced, which characterize how to calculate policy-compliant maximum flows using classic algorithms on an unmodified network.

#### **9. The Fractional Form of the Tinkerbell Map Is Chaotic**

In [8], the authors are concerned with a fractional Caputo-difference form of the well-known Tinkerbell chaotic map. The dynamics of the proposed map are investigated numerically through phase plots, bifurcation diagrams, and Lyapunov exponents considered from different perspectives.

In addition, a stabilization controller is proposed, and the asymptotic convergence of the states is established by means of the stability theory of linear fractional discrete systems. Numerical results are employed to confirm the analytical findings.

**Funding:** This Editorial was supported by: Science Foundation Ireland, by funding Ioannis Dassios under Grant No. SFI/15/IA/3074.

**Acknowledgments:** This issue would not have been possible without the help of a variety of talented authors, professional reviewers, and the dedicated editorial team of Applied Sciences. Thank you to all the authors and reviewers for this opportunity. Finally, thanks to the Applied Sciences editorial team.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Fractional Form of the Tinkerbell Map Is Chaotic**

**Adel Ouannas 1, Amina-Aicha Khennaoui 2, Samir Bendoukha 3, Thoai Phu Vo 4, Viet-Thanh Pham 5,\* and Van Van Huynh <sup>6</sup>**


Received: 3 December 2018; Accepted: 12 December 2018; Published: 16 December 2018

**Abstract:** This paper is concerned with a fractional Caputo-difference form of the well-known Tinkerbell chaotic map. The dynamics of the proposed map are investigated numerically through phase plots, bifurcation diagrams, and Lyapunov exponents considered from different perspectives. In addition, a stabilization controller is proposed, and the asymptotic convergence of the states is established by means of the stability theory of linear fractional discrete systems. Numerical results are employed to confirm the analytical findings.

**Keywords:** fractional discrete calculus; discrete chaos; Tinkerbell map; bifurcation; stabilization

#### **1. Introduction**

Throughout the last 50 years, chaotic dynamical systems have attracted increasing attention due to their applicability in a range of diverse and multidisciplinary fields. A dynamical system is said to be chaotic if its states are extremely sensitive to small variations in the initial conditions. Another important property of chaotic systems is that they have attractors characterized by a complicated set of points with a fractal structure commonly referred to as a strange attractor. This chaotic behavior was first observed in continuous dynamical systems and was thought to be an undesirable property. The first chaotic system encountered in the modeling of a real-life phenomena is that of Lorenz [1], which describes atmospheric convection. Soon after, researchers found that chaotic systems can also be discrete. A number of chaotic maps were proposed throughout the years including the Hénon map [2], the logistic map [3], the Lozi map [4], the 3D Stefanski map [5], the Rössler map [6], and many more. Recently, nonlinear oscillations on Riemannian manifolds that can exhibit a chaotic behavior were introduced in [7,8]. Other related works include an investigation of the chaotic dynamics in a fractional love model with an external environment, as in [9], and an extension using a fuzzy function [10].

In recent years, with the growing interest in fractional discrete calculus [11], people have started looking into fractional chaotic maps. Although fractional maps come with considerable added complexity, they provide better flexibility in the modeling of natural phenomena and lead to richer dynamics with more degrees of freedom. Among the fractional chaotic maps that have been proposed, studied, and applied over the last five years are the fractional logistic map [12], the fractional Hénon

map [13], the generalized hyperchaotic Hénon map [14], and the fractional unified map [15]. Perhaps the main concern of the research community has been the possibility of controlling and synchronizing these types of maps [15–20]. An application of a generalized fractional logistic map to data encryption and its FPGA implementation was achieved in [21].

In this paper, we are interested in the Tinkerbell discrete-time chaotic system, which is of the form:

$$\begin{cases} \mathbf{x}\left(n+1\right) = \mathbf{x}^2\left(n\right) - y^2\left(n\right) + a\mathbf{x}\left(n\right) + \beta y\left(n\right),\\ \mathbf{y}\left(n+1\right) = 2\mathbf{x}\left(n\right)\mathbf{y}\left(n\right) + \gamma\mathbf{x}\left(n\right) + \delta y\left(n\right), \end{cases} \tag{1}$$

where *α*, *β*, *γ*, and *δ* are system parameters and *n* represents the discrete iteration step. It is rumored that the map (1) derives its name from the famous Cinderella story, as the trajectory followed by the map resembles that of Tinkerbell appearing in the movie adaptation of the fairy tale. The Tinkerbell map has been studied by many as it exhibits very rich dynamics including a chaotic behavior and a range of periodic states. For instance, its bifurcation subject to different scenarios and initial settings has been studied in [22–25]. A more comprehensive study was performed in [26]. The authors identified conditions for the existence of fold bifurcation, flip bifurcation, and Hopf bifurcation in the Tinkerbell map.

In order to visualize the dynamics of the map (1), we resort to phase plots, bifurcation diagrams, and Lyapunov exponent estimation. We assume parameter values *α* = 0.9, *β* = −0.6013, *γ* = 2, and *δ* = 0.5 and initial states (*x* (0), *y* (0)) = (−0.72, −0.64). The results are depicted in Figure 1. The Tinkerbell map's phase plot is depicted in Figure 1a. Based on Figure 1b, we can see that the estimated Lyapunov exponents of (1) are given by *λ*<sup>1</sup> ≈ 0.2085 and *λ*<sup>2</sup> ≈ −0.4925. It is well known that a positive Lyapunov exponent indicates a chaotic behavior. The remaining parts of Figure 1 depict the bifurcation diagrams of the map (1) with respect to different parameters. These diagrams confirm that the map exhibits a range of different behaviors.

It should be clear to the reader that the Tinkerbell map has rich dynamics and is heavily dependent on its parameters, as well as the initial setting. The main objective of this paper is to investigate the fractional Caputo-difference form of the Tinkerbell map in order to benefit from the added degrees of freedom due to the fractional nature. It is expected that the fractional form will have even richer dynamics and may consequently be more suitable for applications that require a higher entropy level such as data/image encryption.

**Figure 1.** (**a**) Attractor of the Tinkerbell map (2) with (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5) and initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64). (**b**) Estimated Lyapunov exponents by means of the Jacobian matrix method. (**c**) Bifurcation plot with *α* ∈ [−0.5, 1] as the critical parameter and Δ*α* = 0.0075. (**d**) Bifurcation plot with *β* ∈ [−0.6, −0.1] as the critical parameter and Δ*β* = 0.0025. (**e**) Bifurcation plot with *γ* ∈ [0, 2.1] as the critical parameter and Δ*γ* = 0.01. (**f**) Bifurcation plot with *δ* ∈ [−1, 0.6] as the critical parameter and Δ*δ* = 0.008.

#### **2. Fractional Tinkerbell Map**

In this section, we use recent developments in fractional discrete calculus to define the Caputo-difference fractional map corresponding to (1). First, let us define the *υ*th fractional sum of anarbitrary function *X* (*t*) [27] as:

$$
\Delta\_a^{-\upsilon} X\left(t\right) = \frac{1}{\Gamma\left(\upsilon\right)} \sum\_{s=a}^{t-\upsilon} \left(t - s - 1\right)^{\left(\upsilon - 1\right)} X\left(s\right), \tag{2}
$$

for *<sup>t</sup>* <sup>∈</sup> <sup>N</sup>*a*+*n*−*<sup>υ</sup>* and *<sup>υ</sup>* <sup>&</sup>gt; 0, where <sup>N</sup>*<sup>a</sup>* :<sup>=</sup> {*a*, *<sup>a</sup>* <sup>+</sup> 1, *<sup>a</sup>* <sup>+</sup> 2, ...}. Note that the term *<sup>t</sup>* (*υ*) is known as the falling function and may be defined by means of the Gamma function Γ as:

$$t^{(v)} = \frac{\Gamma\left(t+1\right)}{\Gamma\left(t+1-v\right)}.\tag{3}$$

Based on this definition of the fractional sum, we may define the Caputo-like fractional difference operator.

In this section, we would like to produce a fractional difference form of the Tinkerbell map (1). First, we take the difference form, which for function *<sup>x</sup>* (*t*) : <sup>N</sup>*<sup>a</sup>* <sup>→</sup> <sup>R</sup> with fractional order *<sup>υ</sup>* ∈ <sup>N</sup> is given by:

$$\, \, ^C \Delta\_a^v \ge (t) = \Delta\_a^{-(n-v)} \Delta^n \ge (t) \,. \tag{4}$$

Substituting yields the final form proposed in [28], which is defined as:

$$\, \, ^C \Delta\_{\mathfrak{a}}^v \ge \left( t \right) = \frac{1}{\Gamma \left( n - v \right)} \sum\_{s=a}^{t - \left( n - v \right)} \left( t - s - 1 \right)^{\left( n - v - 1 \right)} \, \, \Delta^v \ge \left( s \right), \tag{5}$$

where *<sup>t</sup>* <sup>∈</sup> <sup>N</sup>*a*+*n*−*<sup>υ</sup>* and *<sup>n</sup>* <sup>=</sup> *υ* <sup>+</sup> 1.

We are now ready to examine the fractional map. First, we take the difference form of (1) to obtain:

$$\begin{cases} \Delta \mathbf{x}\left(n\right) = \mathbf{x}^2\left(n\right) - y^2\left(n\right) + \left(a - 1\right)\mathbf{x}\left(n\right) + \beta y\left(n\right),\\ \Delta y\left(n\right) = 2\mathbf{x}\left(n\right)y\left(n\right) + \gamma \mathbf{x}\left(n\right) + \left(\delta - 1\right)y\left(n\right). \end{cases} \tag{6}$$

We may replace the standard difference in (6) with the Caputo-difference, which yields:

$$\begin{cases} \, \, ^C \mathrm{A}\_d^v \mathbf{x} \left( t \right) = \mathbf{x}^2 \left( t - 1 + \upsilon \right) - y^2 \left( t - 1 + \upsilon \right) \\ \qquad + \left( \mathfrak{a} - 1 \right) \mathbf{x} \left( t - 1 + \upsilon \right) + \beta y \left( t - 1 + \upsilon \right), \\ \, \, ^C \mathrm{A}\_d^v y \left( t \right) = 2 \mathbf{x} \left( t - 1 + \upsilon \right) y \left( t - 1 + \upsilon \right) + \gamma \mathbf{x} \left( t - 1 + \upsilon \right) \\ \qquad + \left( \delta - 1 \right) y \left( t - 1 + \upsilon \right), \end{cases} \tag{7}$$

for *<sup>t</sup>* <sup>∈</sup> <sup>N</sup>*a*+1−*υ*, 0 <sup>&</sup>lt; *<sup>υ</sup>* <sup>≤</sup> 1, *<sup>a</sup>* is the starting point, and *<sup>C</sup>*Δ*<sup>υ</sup> <sup>a</sup>* is a Caputo-like difference operator. The case *υ* = 1 corresponds to the non-fractional scenario (1).

#### **3. Dynamics of the Fractional Tinkerbell Map**

In this section, we will employ numerical tools to assess the dynamics of the proposed fractional Tinkerbell map (7). For that, we will need a discrete numerical formula that allows us to evaluate the states of the map in fractional discrete time. According to [29] and other similar studies, we can evaluate (7) numerically as:

$$\begin{cases} \mathbf{x}\left(n\right) = \mathbf{x}\left(0\right) + \frac{1}{\Gamma(v)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v)}{\Gamma(n-j+1)} \left[\mathbf{x}^2\left(j-1\right) - y^2\left(j-1\right)\right] \\ \qquad + \left(n-1\right)\mathbf{x}\left(j-1\right) + \beta y\left(j-1\right)\right], \\\ y\left(n\right) = y\left(0\right) + \frac{1}{\Gamma(v)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v)}{\Gamma(n-j+1)} \left[2\mathbf{x}\left(j-1\right)y\left(j-1\right)\right] \\ \qquad + \gamma \mathbf{x}\left(j-1\right) + \left(\delta - 1\right)y\left(j-1\right)\right], \end{cases} \tag{8}$$

where we assumed *a* = 0 for simplicity. This yields an initial-value problem similar to that of [30], which allows us to use a similar discrete integral equation.

Using Formula (8), we may obtain the states of the fractional Tinkerbell map and consequently produce time series plots of the states, phase-space plots, and bifurcation diagrams. We start with a simple case where the parameters and initial conditions are identical to those adopted in the standard case, i.e., (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5) and (*x* (0), *y* (0)) = (−0.72, −0.64). Given the fractional order *υ* = 0.98, Figure 2 depicts the discrete time evolution of the states. Since the time series in Figure 2 do not indicate the existence or absence of chaos definitively, it is more convenient to show the trajectories followed by the map in state space. Figure 3 shows the phase plots for different values of the fractional order *υ* ∈ {0.995, 0.99, 0.97, 0.952}. We see that the overall Tinkerbell shape remains

valid for a short range of fractional orders. As the order gets close to 0.95, the trajectory almost completely disappears.

**Figure 2.** Time evolution of the fractional Tinkerbell map's states with parameters (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and fractional order *υ* = 0.98.

**Figure 3.** Phase plots of the fractional Tinkerbell map (7) for parameters (*α*, *β*, *γ*, *δ*) = (0.9,−0.6013, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and different fractional orders.

Although the phase plots give an indication of the behavior of the map, it is not until we visualize the bifurcation of the map subject to different parameters that a more complete picture forms. We choose the parameter *β* as the critical parameter and varied it over the range *β* ∈ [−0.6, −0.1] in steps of Δ*β* = 0.0025. The process may be easily repeated for other parameters. The bifurcation

diagrams obtained using the same parameter and initial condition values from earlier are depicted in Figure 4. We observe that although the general dynamics remain similar, the intervals seem to become shorter as the fractional order is decreased.

Even though these bifurcation diagrams suggest the existence of chaos in the fractional Tinkerbell map, they are not definitive. Generally, in order to prove the existence of chaos, we must use multiple tools including time series, phase portraits, Poincaré maps, power spectra, bifurcation diagrams, Lyapunov exponents, etc. The next tool at our disposal is Lyapunov exponents. We calculate these exponents by means of the Jacobian method. It is well known that when *λmax* is positive and the points in the corresponding bifurcation diagram are dense, the map is highly likely to be chaotic. Figure 5 shows the largest Lyapunov exponents corresponding to the same bifurcation diagrams depicted in Figure 4 in the *x*-*β* plane. We can observe clearly that for certain ranges of the parameter *β*, chaos exists.

**Figure 4.** Bifurcation diagrams of the fractional Tinkerbell map (7) with *β* ∈ [−0.6, −0.1] being changed in steps of Δ*β* = 0.0025, parameters (*α*, *γ*, *δ*) = (0.9, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and different fractional orders.

**Figure 5.** The largest Lyapunov exponent as a function of parameter *β* for different values of the fractional order.

Another interesting aspect is the effect of the fractional order on the dynamics of the map for a specific set of parameter values. We fix the parameters and initial conditions at (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5) and (*x* (0), *y* (0)) = (−0.72, −0.64), respectively. Figure 6 shows the bifurcation plot with the critical parameter *υ* ∈ [0, 1] being changed in steps of Δ*υ* = 0.005. This is interesting in that it shows that although the chaotic behavior disappears when the fractional order drops close to 0.95, it is observed again over intermittent intervals. Chaos does not disappear totally until the fractional order is very low.

**Figure 6.** Bifurcation diagram of the fractional Tinkerbell map (7) with *υ* ∈ [0, 1], Δ*υ* = 0.005, (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5), and (*x* (0), *y* (0)) = (−0.72, −0.64).

The largest Lyapunov exponent corresponding to the this bifurcation diagram in the *x*-*υ* plane is depicted in Figure 7. From the figure, we observe that for a fractional order larger than 0.952, *λmax* is positive, which implies that the fractional Tinkerbell map is chaotic. During the interval (0.6609, 0.952),

*λmax* is observed to change intermittently between positive and negative signs, which means that chaos starts to appear and disappear. Finally, for values lower than 0.6609, chaos disappears completely. These results agree with the bifurcation diagram in Figure 6.

**Figure 7.** The largest Lyapunov exponent as a function of the fractional order *υ* for the same parameters and initial conditions in Figure 6.

#### **4. Control of the Fractional Tinkerbell Map**

In this section, we show that the proposed fractional Tinkerbell can be stabilized by means of a simple adaptive feedback controller. In order to be able to establish the asymptotic convergence of the controlled states towards zero, we first need to recall some important results from the literature concerning the asymptotic stability of fractional discrete systems. Since fractional discrete calculus is still relatively new, the existing literature related to stability is very limited. There are two main ways of establishing asymptotic stability. The first relies on the linearity of the system and places conditions on the eigenvalues of the Jacobian [31]. The second scheme is a generalization of the well-known Lyapunov direct method [32]. Although, the Lyapunov method is powerful and can support different types of systems, its has yet to be established for delayed fractional discrete systems, which renders it unusable for the system at hand. Hence, our objective here is to design the control laws to linearize the system, which will allow us to use the stability theory of linear systems. The following theorem summarizes the result of [31].

**Theorem 1.** *The zero equilibrium of the linear fractional discrete system:*

$${}^{\mathbb{C}}\Delta\_{a}^{\upsilon}F\left(t\right) = MF\left(t+\upsilon-1\right),\tag{9}$$

*where F*(*t*) = (*f*1(*t*), ..., *fn*(*t*)) *T,* <sup>0</sup> <sup>&</sup>lt; *<sup>υ</sup>* <sup>≤</sup> <sup>1</sup>*, and <sup>M</sup>* <sup>∈</sup> <sup>R</sup>*n*×*n, is asymptotically stable if the eigenvalues <sup>λ</sup> of M satisfy:*

$$\lambda \in \left\{ z \in \mathbb{C} : \ |z| < \left( 2 \cos \frac{|\arg z| - \pi}{2 - \upsilon} \right)^{\upsilon} \text{ and } |\arg z| > \frac{\upsilon \pi}{2} \right\} \tag{10}$$

*for all t* <sup>∈</sup> <sup>N</sup>*a*+1−*υ.*

Consider the controlled version of (7) given by:

$$\begin{cases} \, \, ^\mathbb{C} \Delta\_\mathbf{u}^\mathbf{v} \left( t \right) = \mathbf{x}^2 \left( t - 1 + \upsilon \right) - y^2 \left( t - 1 + \upsilon \right) + \left( \mathfrak{a} - 1 \right) \mathbf{x} \left( t - 1 + \upsilon \right) \\ \qquad + \beta y \left( t - 1 + \upsilon \right) + u\_\mathcal{x} \left( t - 1 + \upsilon \right), \\ \, ^\mathbb{C} \Delta\_\mathbf{u}^\mathbf{v} y \left( t \right) = 2 \mathbf{x} \left( t - 1 + \upsilon \right) y \left( t - 1 + \upsilon \right) + \gamma \mathbf{x} \left( t - 1 + \upsilon \right) \\ \qquad + \left( \delta - 1 \right) y \left( t - 1 + \upsilon \right) + u\_\mathcal{y} \left( t - 1 + \upsilon \right), \end{cases} \tag{11}$$

where *ux* (*t*) and *uy* (*t*) are adaptive control terms. The following theorem presents the proposed control laws.

**Theorem 2.** *The states of the controlled 2D fractional Tinkerbell map (11) are guaranteed to converge towards zero asymptotically subject to:*

$$\begin{cases} \boldsymbol{u}\_{x}\left(t\right) = \boldsymbol{y}^{2}\left(t\right) - \boldsymbol{x}^{2}\left(t\right),\\ \boldsymbol{u}\_{y}\left(t\right) = -2\boldsymbol{x}\left(t\right)\boldsymbol{y}\left(t\right) - \gamma\boldsymbol{x}\left(t\right). \end{cases} \tag{12}$$

**Proof.** Substituting (12) into (11) yields the dynamics:

$$\begin{cases} \, \, ^{\mathbb{C}}\Delta\_{a}^{v}\mathbf{x}\left(t\right) = \left(\mathfrak{a} - 1\right)\mathbf{x}\left(t - 1 + \upsilon\right) + \beta\mathfrak{y}\left(t - 1 + \upsilon\right),\\ \, ^{\mathbb{C}}\Delta\_{a}^{v}\mathbf{y}\left(t\right) = \left(\delta - 1\right)\mathbf{y}\left(t - 1 + \upsilon\right), \end{cases} \tag{13}$$

or more compactly:

$$\prescript{C}{}{\Delta\_a^v} \left( \mathbf{x} \left( t \right), \mathbf{y} \left( t \right) \right)^T = A \left( \mathbf{x} \left( t \right), \mathbf{y} \left( t \right) \right)^T,\tag{14}$$

with:

$$A = \begin{pmatrix} \alpha - 1 & \beta \\ 0 & \delta - 1 \end{pmatrix}. \tag{15}$$

The eigenvalues of Aare simply *λ*<sup>1</sup> = *α* − 1 and *λ*<sup>2</sup> = *δ* − 1. It is straight forward to see that these eigenvalues satisfy the conditions of Theorem 1. Consequently, the zero solution of (13) is asymptotically stable, and the states of the controlled map (11) are asymptotically stabilized.

The result of Theorem 2 can be easily put to the test. Consider, for instance, parameters (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and fractional order *υ* = 0.98. Using a modified version of the numerical formula (8), we obtain the states depicted in Figure 8. Clearly, the states do converge towards the all-zero solution. Although the convergence was only established for the commensurate case, experiments have shown that the proposed control laws are also valid for the incommensurate case. Figure 9 shows the controlled states with the same parameters and initial conditions from above, but with different fractional orders (*υ*1, *υ*2) = (0.99, 0.95). Again, we see that the states do in fact converge towards zero, indicating successful stabilization. However, it is apparent that the convergence happens faster in the commensurate case where *υ*<sup>1</sup> = *υ*2.

**Figure 8.** Stabilized states of the controlled fractional Tinkerbell map (11) with parameters (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and fractional order *υ* = 0.98.

**Figure 9.** Stabilized states of the controlled fractional Tinkerbell map (11) with parameters (*α*, *β*, *γ*, *δ*) = (0.9, −0.6013, 2, 0.5), initial conditions (*x* (0), *y* (0)) = (−0.72, −0.64), and fractional orders (*υ*1, *υ*2) = (0.99, 0.95).

#### **5. Conclusions**

In this paper, we have considered a fractional Caputo-difference form of the standard Tinkerbell chaotic map, which is well known for its rich dynamics and interesting characteristics. The dynamics of the fractional Tinkerbell map were investigated numerically using phase plots, bifurcation diagrams, and Lyapunov exponents. Through this investigation, we observed that the fractional order has a significant effect on the fractional map's dynamics. This confirms what has been reported previously in the literature and suggests that the fractional map is superior to the standard one, as it includes more degrees of freedom.

We have also introduced a feedback linearization stabilizing controller for the proposed map and established the asymptotic convergence of the states towards the all-zero solution by means of the stability theory of linear fractional discrete systems. The success of the proposed scheme was demonstrated through numerical simulations both in the commensurate and incommensurate cases.

Although feedback linearization is simple to design and implement, its practicality has been challenged by many in the control and cybernetics research communities. For future work, we plan to investigate other control schemes that can perform better in terms of the power consumption and other essential criteria. The main challenge that we anticipate is the limited literature concerning the stability of fractional discrete systems, especially on the Lyapunov method. This would have to be addressed in order to be able to establish the convergence of any new control scheme.

**Author Contributions:** Conceptualization, A.O.; investigation, A.-A.K.; methodology, A.O. and A.-A.K.; project administration, V.-T.P.; resources, S.B.; software, S.B. and T.P.V.; supervision, V.V.H.; validation, T.P.V.; writing, original draft, V.-T.P.; writing, review and editing, V.V.H.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Policy-Compliant Maximum Network Flows**

#### **Pieter Audenaert \*, Didier Colle and Mario Pickavet**

Faculty of Engineering and Architecture, Department of Information Technology, Ghent-University-imec, 9052 Ghent, Belgium; didier.colle@ugent.be (D.C.); mario.pickavet@ugent.be (M.P.)

**\*** Correspondence: Pieter.Audenaert@UGent.be; Tel.: +32-(0)9-33-14900

Received: 19 January 2019; Accepted: 23 February 2019; Published: 28 February 2019

**Abstract:** Computer network administrators are often interested in the maximal bandwidth that can be achieved between two nodes in the network, or how many edges can fail before the network gets disconnected. Classic maximum flow algorithms that solve these problems are well-known. However, in practice, network policies are in effect, severely restricting the flow that can actually be set up. These policies are put into place to conform to service level agreements and optimize network throughput, and can have a large impact on the actual routing of the flows. In this work, we model the problem and define a series of progressively more complex conditions and algorithms that calculate increasingly tighter bounds on the policy-compliant maximum flow using regular expressions and finite state automata. To the best of our knowledge, this is the first time that specific conditions are deduced, which characterize how to calculate policy-compliant maximum flows using classic algorithms on an unmodified network.

**Keywords:** communication networks; maximum flow; network policies; algorithms

#### **1. Introduction**

Connecting two nodes in a computer communication network involves setting up paths, possibly more than one. Often, this is done with resiliency in mind: clearly, having more than one path between nodes, it might be possible to route around a failed link or node, depending on the paths themselves and the failed node or link. Also, network operators provide multiple paths between nodes in order to increase the maximum achievable throughput or flow in between these nodes.

In practice, however, network policies are in effect, effectively restricting the flow that can be set up between two nodes in the network. Specific restrictions can be implemented to optimize routing, due to security constraints or because of policies and agreements between different network operators. As one example, inter-domain paths in the internet often fulfill a valley-free routing constraint [1,2] which is a simple condition that models real-life business agreements between different operators. Compared to a policy-free routing model, these conditions severely restrict allowed paths and might have a large impact on the overall throughput between two nodes when the conditions are strictly enforced. Conditions also influence path diversity, which in turn has a large impact on the overall resiliency of the connection.

Calculating the maximum throughput that can be achieved between two nodes in a network can easily be done via classic techniques solving the maximum flow problem, such as the algorithm of Ford–Fulkerson or the push-relabel method [3,4]. However, these algorithms do not take into account policies and do not enforce any condition or constraint at all. In this work, we will adapt generic flow algorithms, such that the flow that is achieved actually is a policy-compliant maximum flow. Policies are defined using finite state automata (FSA) [5,6], as their expressive power is sufficient for many purposes in real life, while still being conceptually simple.

We will define an algorithm that exactly solves the policy-compliant maximum flow problem, however, at the expense of a large computational footprint. Then, we will continue to adapt the classic

algorithms and datastructures such that they take into account policies that can be expressed using FSA. These algorithms, however, do not guarantee exact solutions to the problem anymore; they merely provide lower bounds to the exact answer. At each step, we are able to tighten the lower bound, such that it comes closer to the exact solution. However, each step also implies additional computational work, so we are trading off the quality of the results with the time needed to run the algorithms.

The rest of the paper is organized as follows. First, in Section 2 we provide some background about policies and how they can model certain conditions and constraints that are applied in real life. We also point out some difficulties that arise when trying to adapt the classic algorithms. Then, in Section 3 we discuss the formal model and the classic algorithms from literature, and how we can formally model policy-compliant connections. Next, in Section 4 we adapt the algorithms that were discussed in Section 3 to obtain a series of techniques that allow us to calculate lower bounds to the problem. Finally, Section 5 discusses the results obtained, and concludes this work, after which we provide an outlook to future work.

#### **2. Background and Context**

Policies are very important in today's internet. Specific agreements between operators are often hidden from prying eyes, but they can have a huge impact on how specific paths are set up between them. It turns out that many network policies, used in practice, can be translated in very simple formal languages, all belonging to a class of formal languages called "regular languages" [5,6]. This class is one of the most basic classes of languages used in computer science, as it is non-trivial but contains a very broad range of languages with widely differing properties that are very useful in practice.

The complexity of properties which we can express using regular languages is larger than would be expected at first sight. To provide a non-trivial, albeit contrived example, consider a government that wants to contact its embassy located in another country. Some links are trusted, as they are operated by friendly states, but some links in between may be eavesdropped upon by a rival state. Traffic crossing such links needs to be encrypted; once encrypted, traffic may be routed through any link in the network. The receiving party will need to decrypt the message, before it can be distributed to internal data-processing departments. The government thus places into force the following conditions: traffic from the government to the embassy should either be sent through links owned exclusively by friendly states, or it should pass the encryption box, after which it can be routed through any link. However, if encrypted, the message has to be decrypted before delivery at the internal departments at the embassy, i.e., pass through a specific decryption box.

It is clear that enforcing such policies is important, moreover such policies should be formalized instead of specified in vague terms as above. Regular languages fit that job quite well, as they allow policies like the one above, next to other much more involved policies, to be specified without any possibility for misunderstanding. Moreover, once specified, network-routes can be automatically validated against a set of policies, and compliance can be trivially checked.

Taking into account policies in the network, a natural question is to ask for their impact on the routing between nodes in the network. More specifically, in this work we are mainly interested in their impact on the maximum flow achievable between two nodes. Indeed, previously acceptable paths suddenly do not obey the specific policies, and as such, some routes, maybe all, will become invalid. Thus, the policy-compliant maximum flow between two nodes needs to be redefined and recalculated. Gaining insight into the flow-structure between two nodes also tells us something about resiliency in the network. Operators might be convinced their network has sufficient spare capacity, however, depending on the policy at hand and the specific structure of the network, a failing link might have tremendous consequences as re-routing traffic around the failure might invalidate the policy which is unacceptable.

We thus come to the following main problem-statement that will be treated in this work: how to calculate or approximate the policy-compliant maximum flow between two nodes?

Clearly, policy-compliance and both the theoretical and practical consequences are important in the design and operation of networks. From a practical point of view, much information can be found about policies and being compliant, and multiple frameworks for the management and monitoring of policies exist and are in current use. In contrast, few theoretical models have been developed, and literature is sparse when taking a more fundamental approach to the problem.

Caesar and Rexford [7] discuss how routing policies came into existence, and how the protocols have evolved over time and became increasingly complex. They discuss how common policies are implemented and address the problems that arise in applying and supporting policies. Feamster et al. [8] discuss fundamental objectives for interdomain routing and traffic engineering and provide practical guidelines. They show how greater flexibility can be gained in several situations and demonstrate the manipulation of traffic via small changes in specific routes. The paper by Hu et al. [9] proposes an approach to overcome the inherent constraints of compliant recovery schemes. Adapting protocols, they succeed in improving route diversity, in turn increasing resilience.

Klöti et al. [10] proposed a graph-transformation technique that constructs a tensor product of a graph and a finite state automaton. They show how to model policies using FSA and apply standard flow maximization algorithms on the transformed graph in order to obtain bounds for the policy-compliant maximum flow problem. Sobrinho et al. [11] provide a deep mathematical analysis to gain insight in the inner workings of route-vector protocols. They relate this to a class of routing policies and quantify how much intrinsic connectivity is lost due to typical routing policies. Erlebach et al. [12] approximate the maximum number of edge- and node-disjoint valid paths between two nodes, via an involved mathematical theory and specific approximation-algorithms.

Soulé et al. [13] provide a declarative formal language based on logical formulas and regular expressions to express network policies in the Merlin framework. After compilation, a constraint solver allocates paths. Tools are provided to verify whether conditions are violated. Hinrichs et al. [14] introduce a declarative policy management language. They focus on expressive power so that policies can be expressed naturally while still being able to enforce policies efficiently. They focus on enterprise networks, and their design using formal mathematical languages is attractive.

Raghavan et al. [15] provide a practical authenticated source routing system that allows for fine-grained path selection and enforcement of cryptographic policies. Capabilities can be defined and composed resulting in complex statements defining specific policies.

Godfrey et al. [16] introduce so-called pathlet routing, pathlets being building blocks which can be concatenated to obtain end-to-end routes fulfilling policy constraints. Pathlet routing can handle typical policies, but also other recent multipath routing proposals or source-routing approaches. Batista et al. [17] propose a policy-based OpenFlow network management framework using the Ponder policy definition language. They focus on simplicity of the system, trying to infer both static and dynamic conflicts. In [18] they update their work and propose a more theoretical approach to conflict detection using first-order mathematical logic to model flows. A Prolog rule-engine applies condition-action rules to infer problems. Xu and Rexford [19] discuss problems concerning multipath interdomain routing and introduce MIRO (Multi-path Interdomain ROuting), which gives transit domains control over the flow across their infrastructure. MIRO remains backwards compatible with other technologies and offers large flexibility with reasonable overhead.

#### **3. Graph Model and Basic Algorithms**

We now introduce basic concepts that are needed in this work. Intuitively, we can think about a flow network like an oil-producing plant, where oil is drilled at one site and needs to be transported to an oil depot via a set of pipes. Following the flow from the source to the sink, the oil can be split or joined at junctions. Most important here is the fact that the path that was followed by a certain molecule of oil is of no importance: it carries no history and passes through junctions anonymously. In contrast, in this work, we need to pay attention to the specific path that a certain amount of flow has followed: it carries its history with itself, see Figure 1. The amount of history that accompanies

the flow is, however, limited: it is represented via its state, and the number of possible different states is limited by the size of the automaton (cf. infra). To understand how to model policies using FSA, we refer to [10].

**Figure 1.** It doesn't matter for oil whether it followed path abc, path def or path dgc, but in this paper, we do care about the specific paths followed. Policies about which path is acceptable and which path is not, are expressed via constraints on the labels across the paths, i.e., abc, def and dgc.

#### *3.1. Flow Networks and Flows*

We now introduce flow networks and flows, following the approach taken by Cormen et al. [4]. These definitions might differ from other approaches in several aspects, e.g., the fact that the capacity function *c* and the flow *f* are total functions, allowing more concise and elegant definitions of constraints. The rationale behind this way of defining the context is included in Cormen et al. [4] and we refer interested readers to study the mathematical discussions in that book. We thus define a flow network *<sup>G</sup>* = (*V*, *<sup>E</sup>*) as a directed graph together with a capacity function *<sup>c</sup>* : *<sup>V</sup>* <sup>×</sup> *<sup>V</sup>* <sup>→</sup> <sup>R</sup> such that each edge (*u*, *v*) ∈ *E* has a nonnegative capacity *c*(*u*, *v*) ≥ 0. Note that *c* is a total function such that, if (*u*, *v*) ∈/ *E*, we have *c*(*u*, *v*) = 0. Moreover, we require that if *E* contains a certain edge (*u*, *v*), then there is no reverse edge (*v*, *u*) present in *E*. We also disallow loops such that ∀*u* : (*u*, *u*) ∈/ *E*. Two specific nodes in the graph are now chosen out of *V*, one vertex being the source *s* and one vertex being the sink *<sup>t</sup>*. We now want to send a flow from *<sup>s</sup>* to *<sup>t</sup>*, a flow being a total function *<sup>f</sup>* : *<sup>V</sup>* <sup>×</sup> *<sup>V</sup>* <sup>→</sup> <sup>R</sup> subject to the following constraints:

• A capacity constraint which expresses that flow can never exceed capacity:

$$\forall u, v \in V: 0 \le f(u, v) \le c(u, v) \tag{1}$$

• A flow conservation constraint which expresses that for all nodes (except source and sink) the incoming flow equals the outgoing flow:

$$\forall u \in V - \{s, t\}: \sum\_{v \in V} f(v, u) = \sum\_{v \in V} f(u, v) \tag{2}$$

Note that *f*(*u*, *v*) = 0 if (*u*, *v*) ∈/ *E*.

The value | *f* | of a flow is defined as the flow leaving the source minus the flow arriving at the source: | *<sup>f</sup>* | = <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*s*, *<sup>v</sup>*) − <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*v*,*s*). Equivalently, we can define the value of the flow just as well as | *<sup>f</sup>* | = <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*v*, *<sup>t</sup>*) − <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*t*, *<sup>v</sup>*). Note that one can define the value of the flow using a formula containing only *s* but not *t*, or a formula containing only *t* but not *s*; one never needs both *s* and *t* in one formula at the same time, see the book by Cormen et al. [4]. Also note that we allow <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*v*,*s*) as well as <sup>∑</sup>*v*∈*<sup>V</sup> <sup>f</sup>*(*t*, *<sup>v</sup>*) to be strictly larger than 0. The classic problem, given a flow network *<sup>G</sup>*, is to maximize the flow | *f* | from *s* to *t* under the constraints above. (The |.|-notation denotes flow value, not absolute value or cardinality.)

#### *3.2. Residual Graphs and Augmenting Paths*

One line of solving this so-called maximum flow problem is via the introduction of residual graphs and the calculation of augmenting paths. In short, we iteratively increase the value of the flow, starting with ∀*u*, *v* ∈ *V* : *f*(*u*, *v*) = 0. At each step, we calculate an augmenting path in the residual graph, which is nothing but the original graph, extended with return-edges. We will now formally define the necessary concepts.

*Appl. Sci.* **2019**, *9*, 863

First, we define the residual capacity *<sup>c</sup> <sup>f</sup>* : *<sup>V</sup>* <sup>×</sup> *<sup>V</sup>* <sup>→</sup> <sup>R</sup> in the following way. We iterate over all edges (*u*, *v*) ∈ *E* and for each such edge we define the following:

$$\begin{cases} c\_f(u,v) = c(u,v) - f(u,v) \\ c\_f(v,u) = f(u,v) \end{cases} \tag{3}$$

If there is no edge between *u* and *v* in any direction, i.e., (*u*, *v*) ∈/ *E* and (*v*, *u*) ∈/ *E*, then both *c <sup>f</sup>*(*u*, *v*) = 0 and *c <sup>f</sup>*(*v*, *u*) = 0.

Next, we define the residual network *Gf* = (*V*, *Ef*), where *Ef* = {(*u*, *v*) ∈ *V* × *V* : *c <sup>f</sup>*(*u*, *v*) > 0}. So, for each edge in the original graph we have a residual edge if the capacity of that edge is not filled up yet, and moreover we add a reverse edge when it carries any flow at all.

This residual graph will now be used to find augmenting paths. Given a flow network *G* = (*V*, *E*) and a flow *f* , an augmenting path *p* is defined as a path from *s* to *t* in the residual network, that is *p* = [*e*1,*e*2,...,*en*], such that

• edges belong to *Ef* , thus

$$\forall e\_i : e\_i \in E\_{f'} \tag{4}$$

• we start in *s*, thus


$$\forall i \in \{1, \ldots, n-1\}: \exists \mu, \upsilon\_i \, w: \mathfrak{e}\_i = (\mathfrak{u}, \upsilon) \land \mathfrak{e}\_{i+1} = (\upsilon, w) \tag{6}$$

• we end in *t*, thus

$$
\exists \upsilon : \varepsilon\_n = (\upsilon, t). \tag{7}
$$

In order to increase the flow *f* , we will change the flow along the edges of the augmenting path *p* with the maximum amount possible, namely the residual capacity of the path *p* which is

$$\mathfrak{c}\_f(p) = \min \{ \mathfrak{c}\_f(\mathfrak{u}, \mathfrak{v}) : (\mathfrak{u}, \mathfrak{v}) \in p \}. \tag{8}$$

The flow *f* can now be increased to a new flow *f* as follows:

$$f'(u,v) = \begin{cases} f(u,v) + c\_f(p) & \text{if } (u,v) \in p\_\prime \\ f(u,v) - c\_f(p) & \text{if } (v,u) \in p\_\prime \\ f(u,v) & \text{otherwise.} \end{cases} \tag{9}$$

Iteratively finding augmenting paths in a residual graph and increasing the flow will finally come to a stop when there are no augmenting paths anymore. Then, the maximum flow has been reached. Depending on the procedure used to find augmenting paths, the above algorithm is called Ford–Fulkerson or Edmonds–Karp.

#### *3.3. Realization of a Flow*

Apart from calculating the value of the maximum flow in a flow network, we are also interested in the actual paths followed by the flow; this set of paths is called the "realization" of the flow. Augmenting paths are a device to calculate the maximum flow, however, augmenting paths are not real paths in the flow network as they might comprise edges from *Ef* that are not in *E*. Define a path *p*

from *s* to *t* in *G* in exactly the same way as an augmenting path, except that edges in the path should be contained in *E* instead of *Ef* . Thus, in addition to (5)–(7), we replace (4) with

$$
\forall \mathcal{e}\_i \, : \, \mathcal{e}\_i \in E.
$$

Given a flow network *G* and accompanying maximum flow *f* , it is possible to determine a set

$$S = \{ (p\_1, c\_1), (p\_2, c\_2), \dots, (p\_n, c\_n) \} \tag{11}$$

that realizes that flow. Intuitively, we start with no flow in the network at all. Then, for each of the paths *pi* we fill all edges of that path with an additional amount *ci*. After all paths are processed, we obtain the maximum flow *f* .

Construction of this set *S* can be done simultaneously with the calculation of the maximum flow. During execution of the maximum flow algorithm, we perform the following steps every time a new augmenting path *p* with additional flow amount *c* has been found, starting with an empty set *S* = {}.


Thus, each iteration of the maximum flow algorithm results in an augmenting path that is possibly cut in sections to remove reversed edges from the residual graph. Forward sections are glued together with sections from previously found paths. When no augmenting paths can be found, the set *S* will contain paths and associated flow amounts, that together sum up to the maximum flow. That is, if upon termination *S* = {(*p*1, *c*1),...,(*pn*, *cn*)}, we have

$$f(u,v) = \sum\_{\{(p\_i,c\_i) : (u,v) \in p\_i\}} c\_i \tag{12}$$

**Figure 2.** Three sections: forward, backward and forward.

Thanks to the fact that every flow has a realization and the fact that it is easy to derive the flow from a realization, from now on we identify a realization of a flow and the flow itself, and use the two concepts as one.

#### **4. Policy-Compliant Paths and Algorithms**

We now continue to define the concept of policy-compliant paths. Then we want to find the maximum flow in the network, allowing only the use of such policy-compliant paths. We will provide algorithms that calculate lower bounds, as well as exact solutions.

#### *4.1. Finite State Automata*

We now proceed to define finite state automata [5,6], which are the simplest of all classic automata. Note that we will only consider deterministic automata, although for finite state automata we have that non-determinism does not add expressive power, i.e., the class of deterministic finite state automata is exactly the same as the class of non-deterministic finite state automata.

We define a deterministic finite state automaton as *M* = (*Q*, Σ, *δ*, *q*0, *F*), where


Such an automaton works as follows. Initially, we start with an input word *w* ∈ Σ∗, where Σ<sup>∗</sup> is the set of strings obtained by concatenating zero or more symbols from Σ. Each of these symbols of the word *w* gets consumed, from left to right, while the automaton changes its state accordingly to its transition function, taking as input the current state and the symbol just read. If the automaton starts in state *q*0, and finishes reading the word *w* while reaching state *q*, then we say that the word *w* is accepted by the automaton if and only if *q* ∈ *F*. For convenience, we introduce the extended transition function *δ*<sup>∗</sup> : *Q* × Σ<sup>∗</sup> → *Q*, where the second argument is a string rather than a single symbol. We define *δ*∗ recursively as follows:

$$\delta^\*(q, w) = \begin{cases} q & \text{if } w = \lambda, \\ \delta^\*\left(\delta(q, a), v\right) & \text{if } w = a v, \end{cases} \tag{13}$$

where *a* ∈ Σ and *v*, *w* ∈ Σ∗. Note that we also introduced the empty word *λ* out of convenience for this definition. We won't need it any further in this work.

#### *4.2. Labeled Networks and Policy-Compliant Paths*

Given a certain flow network *G* = (*V*, *E*) and an FSA *M*, we define a labeling function *l* : *E* → Σ that attaches one symbol to every edge of the graph, which becomes a labeled flow network. Before, we were interested in all paths from source to sink, in order to maximize the total flow value of the network. Now, we constrain the allowed paths to only those that are accepted by the FSA *M*. Formally, we define a policy-compliant path *p* from *s* to *t* in *G* to be a sequence of edges from *E*, such that the symbols that are encountered, in order, constitute a word that is accepted by *M*. Thus, in addition to (5)–(7) and (10), we add the following constraint:

$$\delta^\*(q\_0, l(e\_1)l(e\_2)\dots l(e\_n)) \in F \tag{14}$$

Thus, we are looking for paths from source to sink that are accepted by the automaton. Now we can define the policy-compliant maximum flow problem (PCMF) as calculating the maximum flow | *f* | from *s* to *t* in a labelled flow network *Gf* , subject to the constraint that all paths that realize the flow *f* are accepted by an FSA *M*, i.e., they are policy-compliant paths.

#### *4.3. Brute Force Approach*

A basic, inefficient way of solving the policy-compliant maximum flow problem is by way of brute force. In its essence, we first calculate all policy-compliant paths from source to sink, and afterward try out all different orders in which these paths can be filled up. The maximum flow will be reached by at least one of these orderings. As such, this approach will result in the exact solution, however, at the expense of large computational complexity and memory requirements. Indeed, as we employ our search via e.g., a depth-first-search, it might be the case that we need to visit a certain vertex *u* ∈ *E* more than once, because the internal state of the FSA might be different. Effectively, we need to build a new graph *GBF* that contains tuples of vertices and states. After appropriately defining the edges between these tuples, a typical depth-first-search might be used to enumerate successively all policy-compliant paths. More formally, we perform the following steps.


$$f\_s'(u,v) = \begin{cases} f\_s(u,v) + m\dot{m}\_i^s & \text{if } (u,v) \in p\_{i\prime}^s\\ f\_s(u,v) & \text{otherwise.} \end{cases} \tag{15}$$

5. After all paths *p<sup>s</sup> <sup>i</sup>* are processed, the final flow *fs* will be a lower bound for the policy-compliant maximum flow for the network *G*. Thus, *max*{| *fs*|} will be the value of the policy-compliant maximum flow.

It is clear that this approach results in the exact value of the solution to the policy-compliant maximum flow problem. However, except for small networks and automata it will quickly take too much space and time to execute.

#### *4.4. First Lower Bound*

To cut down on the needed computational resources, we start from the classic Ford–Fulkerson or Edmonds–Karp algorithm. In its essence, these algorithms start with an empty flow, after which a new augmenting path is looked for in the residual graph. If found, this augmenting path gives rise to an augmented flow, after which the residual graph is updated. Then, a new iteration is started. Once no new augmented paths can be found, the algorithm has reached the maximum flow.

In order to apply this approach to the policy-compliant maximum flow problem, it will be necessary to keep track of how the flow actually is realized, throughout the execution of the algorithm. This relates to the fact that flows cannot merge and split anonymously anymore, as discussed previously.

As a first step, it will be necessary to keep track of the different state-transitions that occur when flow crosses an edge *e* = (*u*, *v*) ∈ *E* in a path *p* and changes state during that transition from state *qu* ∈ *Q* in vertex *u* to state *qv* ∈ *Q* in vertex *v*. To that aim, we define transitions *t* to be members of the set *Q* × *Q* = {(*qu*, *qv*) : *qu* ∈ *Q* ∧ *qv* ∈ *Q*}. For every edge *e* ∈ *E* we will now keep track of all transitions that cross that edge, and thus define the total function *<sup>T</sup>* : *<sup>V</sup>* <sup>×</sup> *<sup>V</sup>* <sup>×</sup> *<sup>Q</sup>* <sup>×</sup> *<sup>Q</sup>* <sup>→</sup> <sup>R</sup>. This function keeps track of the amount of flow that crosses the edge (*u*, *v*) while changing state from *qu* to *qv* as the value *T*(*u*, *v*, *qu*, *qv*). As a direct consequence, we have at all times ∀*u*, *v* ∈ *V* : <sup>∑</sup>*qu*,*qv*∈*<sup>Q</sup> <sup>T</sup>*(*u*, *<sup>v</sup>*, *qu*, *qv*) = *<sup>f</sup>*(*u*, *<sup>v</sup>*) ≤ *<sup>c</sup>*(*u*, *<sup>v</sup>*) which expresses the fact that a flow over an edge can be split up in different transitions, the amount of which add up exactly to the amount of flow crossing that very edge.

Keeping track of this function *T* allows the algorithm to discriminate between the different states in one single vertex, and the amount of the flow that is actually passing through that vertex having that state at that moment. This information is then used to glue together already found policy-compliant paths and a newly found augmenting path. In detail, this procedure works as follows.


$$(u\_{j}, u\_{y}) \in E \land T(u\_{j}, u\_{y}, q\_{u\_{j}}, q\_{u\_{y}}) > 0. \tag{16}$$

Choose such an edge *e* = (*uj*, *uy*) ∈ *E* and state *quy* ∈ *Q*. When the augmenting path *p* is cut at vertex *uj*, we can glue the new flow that comes in via the forward section [*e*1, ... ,*ei*] to the flow that leaves vertex *uj* in state *quj* via edge *e* = (*uj*, *uy*) and reaches state *quy* at the other end in vertex *uy*. This flow is already available and does not concern us anymore.

6. Now, we need to cancel flow that used to cross the edges [(*vk*, *uk*), ... ,(*vj*, *uj*)], that is the backward section of the augmenting path *p*, crossed in the opposite direction. This gives rise to a second condition for this algorithm to find augmenting paths. This section contains a flow, reaching state *quj* at the end. Tracing back where this flow comes from, we require that there must exist some flow passing vertex *vk* in state *qvk* ∈ *Q* which crosses the edges [(*vk*, *uk*), ... ,(*vj*, *uj*)] finally reaching vertex *uj* in state *quj* . This observation immediately translates to the following condition:

$$\exists q\_{\mathbb{D}\_k} \in \mathbb{Q}: \delta^\*(q\_{\mathbb{D}\_k}, l((v\_k, u\_k)) \dots l((v\_j, u\_j))) = q\_{u\_j}. \tag{17}$$

Choose such a *qvk* ∈ *Q*. Then, the flow that is already present in the network and passes vertex *vk* in state *qvk* will be canceled over the section [(*vk*, *uk*),...,(*vj*, *uj*)].

7. We now need to consider the last forward section [*el*, ... ,*en*]. This section starts in vertex *vk* and ends in the sink *t*, thus ∃*v* ∈ *V* : *el* = (*vk*, *v*) and ∃*v* ∈ *V* : *en* = (*v*, *t*). We know that, having reached state *qvk* in *vk*, we will cross edges [*el*, ... ,*en*] and need to reach some state *qt* ∈ *F* to obtain a policy-compliant augmenting path. This immediately gives rise to the third and last condition, which can be stated as

$$\delta^\*(q\_{\mathbb{U}\_{k'}}l(e\_l)\dots l(e\_{\mathbb{U}})) = q\_l \in F. \tag{18}$$

8. Suppose all conditions (16)–(18) for this augmenting path are fulfilled, which means that we have found an augmenting path that can be joined with the already available policy-compliant paths to increase the total policy-compliant flow. We now need to calculate the amount of flow that can be sent over this augmenting path. First, define *cuj* = *min*{*c <sup>f</sup>*(*u*, *v*) : (*u*, *v*) ∈ [*e*1, ... ,*ei*]} to be the capacity that is at most available along the first forward section of the augmenting path and thus at the endpoint *uj* of this section. Likewise, define *cvk* = *min*{*c <sup>f</sup>*(*u*, *v*) : (*u*, *v*) ∈ [*el*, ... ,*en*]} to be the amount of flow that is at most available along the second forward section of the augmenting path, and thus the maximum amount of flow that can be diverted at the starting point *vk* of this section. For the backward section of the augmenting path, define *cvkuk* = *T*(*vk*, *uk*, *qvk* , *quk* ), ... , *cvjuj* = *T*(*vj*, *uj*, *qvj* , *quj* ), to determine the already crossing flow over the backward edges, taking into account the state-transitions which we derived above. Taking everything together, we have

$$c\_f(p) = \min\{c\_{u\_{j'}}c\_{v\_k u\_{k'}}, \dots, c\_{v\_{j'} u\_{j'}} c\_{v\_k}\} \tag{19}$$

(cf. (8)). This denotes the additional flow *c <sup>f</sup>*(*p*), under the assumption that *c <sup>f</sup>*(*p*) > 0, that can now be sent over the newly found augmenting path *p*, such that the resulting flow *f* can be realized via only policy-compliant paths.

9. We now want to update the value of the flow *f* . This can be done via (9). However, we also need to update the specific values of the function *T*, for all edges involved in the augmenting path. First, we update *T* for all edges *e* = (*ux*, *uy*) in the forward section [*e*1, ... ,*ei*] via *T* (*ux*, *uy*, *qux* , *quy* ) = *T*(*ux*, *uy*, *qux* , *quy* ) + *c <sup>f</sup>*(*p*). Thus, for such an edge *e* = (*ux*, *uy*) we add an additional *c <sup>f</sup>*(*p*) units of flow that change state from *qux* to *quy* while crossing edge *e*. For all edges [(*vk*, *uk*), ... ,(*vj*, *uj*)], which are crossed in the reverse direction in the backward section of the augmenting path, we set *T* (*ux*, *uy*, *qux* , *quy* ) = *T*(*ux*, *uy*, *qux* , *quy* ) − *c <sup>f</sup>*(*p*), actually canceling out the pre-existent transitions in that section. Finally, for the edges in the final forward section [*el*, ... ,*en*] we update the transitions via *T* (*ux*, *uy*, *qux* , *quy* ) = *T*(*ux*, *uy*, *qux* , *quy* ) + *c <sup>f</sup>*(*p*). This concludes the updates of

the transitions-function *T* and the work that needs to be done when an augmenting path has been found.


Although the approach described above results in correct lower bounds for the policy-compliant maximum flow problem, we point out that the conditions (16)–(18) that are imposed upon augmenting paths are very restrictive. Indeed, it is easy to describe specific cases where the approach misses additional augmenting paths due to these conditions, which can easily be relaxed at the expense of additional computational work. The next approach, which will result in a second lower bound that is tighter than the first bound, starts by relaxing condition (16).

In its essence, this condition says that for an augmenting path, containing a backward section, it is necessary to glue the immediately previous forward section to a pre-existent policy-compliant path at the cut-vertex, via an exact match of the state in that cut-vertex. Clearly, this is too restrictive a condition as it would be sufficient if the augmenting flow arriving at that vertex could be continued via the same edges, however taking different intermediate states in between.


$$\exists p\_i \in P : \delta^\*(q\_{u\_j}, l(e\_l^i) \dots l(e\_n^i)) \in F. \tag{20}$$

Choose such a *pi* ∈ *P*. That way, we are confident that the prefix *l*(*e*1)... *l*(*ei*) can be joined together with a suffix *l*(*e<sup>i</sup> l* )... *l*(*e<sup>i</sup> <sup>n</sup>*) such that for the whole word holds that *δ*∗(*q*0, *l*(*e*1)... *l*(*ei*)*l*(*e<sup>i</sup> l* )... *l*(*e<sup>i</sup> <sup>n</sup>*)) ∈ *F*.

4. At this point, we have to be careful which policy-compliant flows crossing the backward section [*ej*, ... ,*ek*] we are going to cancel. Indeed, as the first section [*e*1, ... ,*ei*] is glued together with some policy-compliant path *pi* ∈ *P* at vertex *uj* via condition (20), we are looking for states *qvk* ∈ *Q* at vertex *vk* such that, instead of (17) the following formula holds:

$$\exists q\_{\mathbb{D}\_k} \in Q : \delta^\*(q\_{\mathbb{D}\_{k'}} l((v\_{k'}, u\_k)) \dots l((v\_{j'}, u\_{\bar{j}}))) = q\_{u\_{\bar{j}}}^i. \tag{21}$$

Indeed, after crossing section [*ej*, ... ,*ek*] in reverse direction we want to end in state *q<sup>i</sup> uj* instead of *quj* .


Note that, in comparison to the approach for obtaining the first lower bound, condition (20) not only has merits for being more relaxed, it also carries a serious drawback as checking whether condition (20) holds involves more computational work than (16). One possible approach to check whether (20) actually holds, is by a depth-first walkthrough of the realization associated with the flow *f* . However, this is complicated by the fact that the flow might have multiple valid realizations, some of which are compatible with (20), while some are not. We discuss this now in detail.

1. Consider the augmenting path *p* = [[*e*1, ... ,*ei*], [(*uj*, *vj*), ... ,(*uk*, *vk*)], [*el*, ... ,*en*]] from above, where the state at *uj* is defined as *quj* = *δ*∗(*q*0, *l*(*e*1)... *l*(*ei*)). Suppose there exist two policy-compliant paths *pi* ∈ *<sup>P</sup>* and *pj* ∈ *<sup>P</sup>*, such that their states in vertex *uj* are *<sup>q</sup><sup>i</sup> uj* and *q j uj* respectively. Suppose that the realization of the flow that gave rise to these paths *pi* and *pj* is such that *pi* = [... ,(*uj*, *u*1),(*u*1, *u*3),(*u*3, *u*4),(*u*4, *t*)] and *pj* = [... ,(*uj*, *u*2),(*u*2, *u*3),(*u*3, *u*5),(*u*5, *t*)] and moreover that *l*((*uj*, *u*1)) = *a*, *l*((*u*1, *u*3)) = *b*, *l*((*u*3, *u*4)) = *c*, *l*((*u*4, *t*)) = *d*, *l*((*uj*, *u*2)) = *e*, *l*((*u*2, *u*3)) = *f* , *l*((*u*3, *u*5)) = *g* and finally *l*((*u*5, *t*)) = *h* with {*a*, *b*, *c*, *d*,*e*, *f* , *g*, *h*} ⊆ Σ. Then, as *pi* and *pj* are policy-compliant paths in the set *P*, it directly follows that *δ*∗(*q<sup>i</sup> uj* , *abcd*) ∈ *F* and *δ*∗(*q j uj* ,*ef gh*) ∈ *F*, see Figure 3.

**Figure 3.** Fulfilling condition (20).

2. Suppose that moreover it holds that *δ*∗(*q<sup>i</sup> uj* , *abgh*) ∈ *F* and *δ*∗(*q j uj* ,*ef cd*) ∈ *F*, that is, for the paths *pi* and *pj* it doesn't really matter which way they take in the second diamond. However, suppose that *δ*∗(*quj* , *abgh*) ∈ *F* as well as *δ*∗(*quj* ,*ef cd*) ∈ *F* but *δ*∗(*quj* , *abcd*) ∈/ *F* and also *δ*∗(*quj* ,*ef gh*) ∈/ *F*. 3. Now, there is no way to fulfill condition (20), as no policy-compliant path through *uj* that is part of the realization of *f* allows the prefix *l*(*e*1)... *l*(*ei*) to be extended by a suffix to arrive at a final state in the sink. However, if the realization of the flow *f* was such that it contained the policy-compliant paths [... ,(*uj*, *u*1),(*u*1, *u*3),(*u*3, *u*5),(*u*5, *t*)] and [... ,(*uj*, *u*2),(*u*2, *u*3),(*u*3, *u*4),(*u*4, *t*)], then the prefix of the augmenting path *p* could have been extended via a matching suffix to a final state. In short, condition (20) would have been fulfilled.

Thus, although the relaxed condition (20) can be used to generate augmenting paths that would go unnoticed if applying (16), it suffers from the fact that previous decisions (like how to route policy-compliant paths in case of multiple possibilities) can have a large impact on the number of augmenting paths that can be found later on. As such, we have no guarantee that this second approach will always calculate the exact value for the policy-compliant maximum flow problem, but merely a lower bound. However, as (20) and (21) are fulfilled everytime (16) and (17) are fulfilled, we can immediately derive that the second lower bound will be tighter than the first lower bound.

#### *4.6. Third Lower Bound*

Considering the approaches taken to calculate the first and the second lower bound, observe that at several times choices have to be made, e.g., when fulfilling condition (17). Apart from the fact that some of the conditions might be overly restrictive, we face the fact that we might miss augmenting paths during the search when we make a bad choice, as well as the fact that we might need additional iterations to reach the point where no more augmenting paths can be found. For example, consider (17) and suppose that multiple states exist that fulfill the condition. By choosing one of them, an already present policy-compliant path is chosen too, which implies that the amount of flow that can be canceled has an upper bound that is equal to the minimum value of residual capacity across the backward section. In turn, this restricts the total amount of flow that can be sent over the augmenting path due to (19).

Thus, in order to allow as much flow as possible to be sent over the augmenting path that is constructed, one can check at the same time whether other transitions can be applied in parallel, such that the total amount of flow that is canceled across the backward section is maximized. First, we want to allow multiple pre-existing policy-compliant paths that cross the cut-vertex to be involved in the process, such that we do not have to choose at all which path we want to work with in (20). Multiple policy-compliant paths might fulfill the condition, even reaching the same state in the cut-vertex but following another way to the sink. Also, we want to allow multiple states at the end of the backward section to be involved in the process, such that we do not have to choose at all which state we want to work with in (21). Multiple states might fulfill the condition, and they even might be equal to each other even though they are part of different policy-compliant paths. We need to match the possibilities in both sets to maximize the flow that we can cancel along the backward section in the middle of the augmenting path. This being, of course, under the additional condition (18).

We start again with augmenting path *p*, set *P* and set *S*. Now select a subset *P* ⊆ *P* such that

$$\forall p\_{\boldsymbol{i}} \in P^{\boldsymbol{\ell}}: \boldsymbol{\delta}^\*(q\_{\boldsymbol{u}\_{\boldsymbol{j}}}, \boldsymbol{l}(\boldsymbol{e}\_{\boldsymbol{l}}^{\boldsymbol{i}}) \dots \boldsymbol{l}(\boldsymbol{e}\_{\boldsymbol{n}}^{\boldsymbol{i}})) \in F \tag{22}$$

$$\forall p\_i \in P': \delta^\*(q\_{v\_k}^i, l((v\_k, u\_k)) \dots l((v\_j, u\_j))) = q\_{u\_j}^i \tag{23}$$

$$\forall p\_i \in P': \delta^\*(q\_{v\_k}^i, l(\mathfrak{e}\_l) \dots l(\mathfrak{e}\_n)) \in F \tag{24}$$

$$
\Sigma\_{p\_i \in P^\nu \mathcal{C}\_i} \subsetneq \min(c\_{u\_{j'}}, c\_{v\_k}).\tag{25}
$$

This set *P* can be iteratively built, starting from the empty set and choosing and adding *pi*'s while (25) holds. Once the iteration comes to a stop, augment the flow along the path *<sup>p</sup>* with the value <sup>Σ</sup>*pi*∈*P ci* and update *f* and *T*. As the conditions (18), (20) and (21) imply that (22)–(24) will also hold, we obtain once more a tighter bound. However, there is no guarantee that the policy-compliant maximum flow will be reached. This is the most generic specification of a solution to the policy-compliant maximum flow problem that we derived.

#### *4.7. Implementation and Experimental Results*

Implementing classical maximum flow algorithms is not difficult as they are well-studied. Typical approaches include augmenting path algorithms, push-relabel algorithms and others. They have varying complexities in solving the problem. In practice, execution times are also dependent on the structure and density of the network. We chose an augmenting path approach. Finite state automata are typically implemented in a table-driven way, and it is quite straightforward to implement them efficiently, like we did. However, combining a maxflow algorithm with an FSA multiplies the complexities of both. Indeed, sending a flow over an edge multiple times was now allowed, as the states at the beginning and the end of the edge when crossing it for the first time might be different from the states while crossing it a second time. Thus, utmost care is needed to implement these algorithms correctly. All code is freely available in Supplementary Materials.

Although this paper is mainly intended as a theoretical study, we performed a large number of experiments validating our approach and formulas, in order to build some insights in the practical use of our approach, execution times and the accuracy reached. However, these values must be interpreted as rough estimates, and can only serve as rules-of-thumb, as much is dependent on the actual networks and finite state automata used in the experiments. Indeed, for evaluating the work, a network generator was implemented which allowed for the use of random networks having differing numbers of vertices and edges. Of course, the higher the density of edges, the more paths can be found and the higher the calculation time for the varying parts of the algorithms. We experimented with graphs having up to 256 vertices, with a random choice of edges (from one to the maximum number, i.e., *n*(*n* − 1)/2). Also, we implemented an FSA generator, which allowed experimenting with different policies, all expressed as regular expressions. The maximal number of states allowed was maximally of the same order as the number of vertices, and clearly it holds that allowing more states leads to higher computational times. Indeed, a flow can reach a vertex in about *O*(|*V*|) possible states, which also means that an algorithm looking for augmenting paths can be in about *O*(|*V*| <sup>2</sup>) states itself. Also, the flow over any edge is now parameterized by the states of both endpoints of the edge; as we use up to |*V*| possible states this results in the memory-use for function *T* scaling in the order *O*(|*V*| 4).

For our experiments with networks with up to 32 vertices and FSAs with up to 32 states, we found that running the brute-force algorithm either finished within a second, or did not finish even after 30 min, after which we always stopped the experiments. Careful investigation shows that the density of the network, combined with the structure of the FSA, is of prime importance to whether or not the algorithm finishes fast. For all the randomly generated experiments where the brute-force approach finished and provided us with the maximal value, the first lower bound also reached that same maximal value, but even in the worst observed cases in less than a millisecond, which is three orders of magnitude faster than the brute-force approach. The first lower bound algorithm could not find the maximal value only for specifically designed networks and FSAs, such as the network in Figure 3, aimed specifically at fooling the algorithm. For random experiments, the second lower bound algorithm always gives the same result as the first algorithm, taking only a little bit more time to execute. The third lower bound tries out much more combinations and comes closer to the brute-force approach in computational time, and in our random experiments always obtained the same maximal flow values as the second lower bound. We also executed random experiments with up to 256 vertices without having the maximal flow value, because the brute-force algorithm did not terminate. For the subset of these experiments where the three lower bound algorithms did terminate, also no differences were obtained in the resulting values. Thus, differences in the accuracy of results were never detected in any random experiment, and could only be obtained with specifically designed and contrived networks and automata. Runtimes for the first lower bound algorithm, for random networks having up to 256 vertices, are shown in Figure 4. If the algorithm finished, it typically finished fast (i.e., within

about 10 s), as can be seen in the figure. If it didn't finish, we aborted the execution after 30 min (not depicted).

Concluding, we advise to use the first lower bound algorithm, as it terminates very fast and for our random experiments always reached the maximal value when known, i.e., when the brute-force approach also finished. It couldn't find the maximal value only in contrived examples, specifically designed to fool the algorithm. Also, we did not encounter random experiments where the three algorithms gave different results, so there is no reasonable situation where one might recommend using the second or third algorithm. We would like to stress again that the execution times we obtained varied wildly, being highly dependent on the structure of the network, the accompanying FSA and the combinatorial combination of the two.

**Figure 4.** Runtimes for the first lower bound algorithm, for random networks with up to 256 vertices.

#### **5. Conclusions and Future Work**

Knowing the maximum throughput between two nodes in a network is key knowledge for a network operator. In practice, network policies severely complicate this task and often leave the question open-ended. We put forward a formal model with which it is possible to model real-life policy constraints, and we analyzed the impact of policies on throughput. As exact solutions are difficult to obtain, we defined a series of conditions and algorithms that allow us to calculate a sequence of increasingly tighter bounds on the exact value. These algorithms are built upon classic algorithms like Ford–Fulkerson to solve the generic maximum flow problem, which are adapted to our needs and augmented with specific functions to ease bookkeeping of additional data like transitions or the specific realization associated with the flow.

Although the approach to specify the conditions on the augmenting paths in order to enforce their policy-compliance is reasonable in combination with Ford-Fulkerson, it requires complex formulae that express conditions on the augmenting path, end-to-end. Future work in the direction of push-relabel algorithms would solve this problem: on the basis of local conditions, it might be possible to express specific constraints that guarantee path-compliancy upon arrival of the flow.

**Supplementary Materials:** All code is freely available and downloadable at http://www.dna.idlab.ugent.be.

**Author Contributions:** Conceptualization and methodology: P.A., D.C. and M.P.; software, validation, formal analysis, resources, visualization, writing: P.A.; investigation, writing—review and editing: P.A., D.C. and M.P.

**Funding:** This research was funded by Ghent University—imec, IOP 2016 Derudder-Pickavet.

**Acknowledgments:** The authors wish to thank Ghent University—imec for the support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

FSA Finite state automaton

PCMF Policy-compliant maximum flow

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
