{'fraction_non_alphanumeric': 0.05730570406044272, 'fraction_numerical': 0.015464145274599037, 'mean_word_length': 3.9625082218811665, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A dynamic programming approach for controlled fractional SIS models

We investigate a susceptible-infected-susceptible (SIS) epidemic model based on the Caputo-Fabrizio operator. After performing an asymptotic analysis of the system, we study a related finite horizon optimal control problem with state constraints. We prove that the corresponding value function is a viscosity solution of a dynamic programming equation. We then turn to the asymptotic behavior of the value function, proving its convergence to the solution of a stationary problem, as the planning horizon tends to infinity. Finally, we present some numerical simulations providing a qualitative description of the optimal dynamics and the value functions involved.

We study the fractionary SIS system where β, γ, N > 0, and I 0 , S 0 ≥ 0 satisfy I 0 + S 0 = N . Moreover, we make the technical assumption α + (1 − α)(β − γ) ≤ M (α). (2) Note that above condition is trivially satisfied if α = 1, since M (1) = 1 by definition. The functions S(t) and I(t) represent the size of susceptible and infective individuals, respectively. N is assumed to be the size of the total population at initial time. β is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject, and γ is the recovery rate. Fractional epidemic models attracted the interest of researchers, see [CLYL21] for a recent review, due to the possibility of tuning the derivative order α for applications to real data fitting, see for instance [LWLT19]. An important feature is the ability to incorporate memory effects into the model: in particular, we refer [KSAQB19] (and the references therein) for a detailed study of the inclusion of memory in epidemic models by using the Caputo-Fabrizio operator. The peculiarity of this fractional operator, introduced in [CF15], is the presence of a non-singular, exponential kernel. We show below that (1) rewrites as Dynamics as (3) underlie the class of SIS models with saturated incidence rate H α (I) := λαI 1+kαI , and saturated treatment function T α (I) := rαI 1+kαI , see [ZZL14,ZL08]. The system (3) models a stable population in which the natural recovery rate is zero: individuals recover from the disease only if they are treated, and they are healed at the rate T α (I). In particular, r is the cure rate, whereas 1/(1+k α I) measures the reverse effect of the delayed treatment, due, for instance, to the limited health system capacity. Also the incidence rate H α (I) is assumed to saturate as I increases: this models the psychological effect of the awareness in the susceptible population about the existence of a large size of infected individuals, inducing a more cautious behavior [CS78]. It is worth noting that (3) is also related to equations arising in the dynamics between tumor cells, immune-effector cells, and immunotherapy [KP98]. We refer to [Pic12] and the references therein, for an application of control theory to the optimization of cancer therapies in a general, nonlinear setting.
Using the conservation of the population, i.e., the identity S = N −I, the system (1) reduces to a non-linear ordinary differential equation depending on the infected population I(t) only, see Theorem 1. We then complete this equation with a linear control term, and we address the problem of minimizing the size of infected individuals plus a quadratic cost on the control. To this end, we study the associated dynamic programming equation, namely an evolutive Hamilton-Jacobi equation for which the value function u α (x, t) is proved to be a viscosity solution, see Theorem 2. Our main result is Theorem 4, in which we prove, as t → +∞, the convergence of u α (x, t) to the value function v α (x) of an associated stationary problem. Theorem 2 and Theorem 4, and the techniques adopted for their proofs, are inspired by the paper [FIL06]. For a general introduction on the topic, we refer to [Lio82].
Finally, we introduce a suitable finite difference scheme for solving the Hamilton-Jacobi equation and building the corresponding optimal trajectories. Some numerical tests complete the presentation, validating our results and providing a qualitative analysis of the solutions.

The SIS model with Caputo Fabrizio derivative
We begin by noting the indentity for f ∈ H 1 ((a, b)): We have the following result, relating (1) to an algebraic identity and to an ordinary differential equation.
In the remaining case α = 1, we recover the classical logistic equation whose unique solution in [0, +∞) is explicitly given by and satisfies I(t) ≥ 0 for all t > 0 if I 0 > 0.

Equlibria and asymptotic behavior
We introduce the reproduction factor ρ := β/γ and we describe the asymptotic behavior of (5) according to the cases of ρ > 1 and ρ ≤ 1.
Proposition 1. If ρ > 1 (respectively, ρ ≤ 1) then the endemic population E := N (1− 1 ρ ) (resp, the equilibrium 0) is an asymptotically stable equilibrium point for (5). In particular, for I 0 > 0, the corresponding solution I satisfies Proof. We rewrite the second equation in (5) in terms of ρ: We set I e := E if ρ > 1, I e := 0 otherwise, and we define V (x) := 1 2 (x − I e ) 2 . We prove that V is a Lyapunov function for (6). To this end, note that V is a positive definite function in R \ {I e } and that V is smooth in R. Moreover, by a direct computation, we get and this concludes the proof.
Remark 1. This result was earlier proved in [ZZL14] in the framework of saturated SIS models, here we propose an alternative proof for the fractional context under exam.
3. An optimal control problem for fractional SIS Fix β, γ, N > 0, ρ = β/γ as before. We consider a controlled version of the infected population dynamics (5): where x ≥ 0 and ξ ∈ L 1 ((0, T )) is the control function. We consider the following finite horizon optimal control problem: minimize with respect to ξ T 0 I 2 (t) 2 + ξ 2 (t) 2 dt + φ(I(T )) subject to (7) and where T > 0, [0, T ] is the given planning horizon, and φ(I(T )) is the expected future cost, depending on the size of the individuals remaining infected at time T . We assume that φ is a continuous, nonnegative function attaining its global minimum at x = 0. The state constraint I ≥ 0 is introduced for modeling reasons. We prove that the value function u α (x, T ) associated to the problem (8) is a viscosity solution of a related dynamic programming equation.
To this end, we denote by I(t; ξ, x) the absolutely continuous solution of (7). We say that a trajectory-control pair (I(t; ξ, x), ξ(t)) is admissible if ξ ∈ L 1 ((0, T )) and I(t; ξ, x) ≥ 0 for all t ∈ (0, T ), and we denote by A T ⊂ L 1 ((0, T )) the set of admissible controls. Then, we define the value function associated to (8) and we consider the following Hamilton-Jacobi equation in which the Hamiltonian is provided by the Legendre transform where the supremum is achieved by the optimal control ξ = −Du.
Remark 2. Equation (10) is the dynamic programming equation associated to the problem of minimizing with respect to ξ T 0 I 2 (t) 2 + ξ 2 (t) 2 dt + φ(I(T )) subject to (7) , namely a version of (8) in which the state constraint I(t) ≥ 0 is removed. The possibility to neglect the state constraint, combined with the unboundedness of the control set and the nonlinearity of the drift b α (which is not Lipschitz when α = 1) are elements of novelty that require an ad hoc analysis.
In agreement with the classical theory, we have the following result.
Proof. The proof is based on showing that u α simultaneously fulfills the definition of viscosity sub-solution and of viscosity super-solution. Technical computations can be easily adapted from [FIL06, Theorem 10] (see also the proof of Theorem 3 below) and we omit them for brevity. The main ingredients are the Dynamic Programming Principle (Proposition 2) and the continuity of u α (Proposition 3), that are proved in detail the following subsections.
3.1. Auxiliary results for the proof of Theorem 2 Let R + := (0, +∞). For (x, T ) ∈ [0, +∞) × R + , let C + (x, T ) denote the space of non-negative, absolutely continuous functions X : [0, T ] → [0, +∞) such that X(0) = x. Then, the value function u α defined in (9) can be rewritten in a form suitable for our purposes: The proof of Theorem 2 relies on several preliminary results, that can be summarized in a set of estimates for u α , proved in Section 3.1.1, the Dynamic Programming Principle, and a continuity result for u α , proved in Section 3.1.2 below.
Proof. The lower estimate readily follows by the assumption that φ attains its global minimum at x = 0. Indeed, for all T ≥ 0 and ε > 0 there exists X ∈ C + (x, T ) such that To prove the upper estimate, consider the curve X ∈ C + (x, T ) defined by Proof. Let X ∈ C + (x, T ) satisfy (16). Then X(t) ≥ 0 by the definition of C + (x, T ). To prove X(t) ≤ C α 2 (R), note that by Lemma 1 In particular, also recalling that we assumed φ to be non negative, for all From this we deduce Integrating the left handside of above expression, we get Since 3.1.2. Continuity of u α . The next two results investigate the dependence of u α on the initial datum φ and the local uniform continuity of u α with respect to the space variable x, respectively. For any continuous function We denote by ω φ,R the modulus of continuity of φ restricted to [0, R].
Proof. Assume that X ∈ C + (x, T ) satisfies By Lemma 1 and Lemma 2, we respectively get Moreover, it follows from the latter inequality -see also the definition of Moreover, since for all A > 0 and y ∈ R and this, together with (20) implies (18). Now, by Lemma 3, we readily get In view of the arguments above, X verifies (18). This, together with the arbitrariness of ε, implies and completes the proof.
Proof. We first let T ≤ |x−y|. Taking ν R as in Lemma 4, and by enlarging it if necessary, we may assume without loss of generality |φ( Assume now T > |x − y| and, by swapping the variable's names x and y if necessary, assume that u α (x, T ) ≤ u α (y, T ). Fix ε ∈ (0, 1) and select X ∈ C + (x, T ) such that Since ε < 1, we know by Lemma 2 that On the other hand, replacing ν R by ν C α 2 (R) if necessary, we have Applying above inequality tox = X(T − |x − y|), we get one gets, also using the last inequality in (21), and this concludes the proof of the local uniform continuity of u α in x. Sincẽ C α (R) is independent from φ, we deduce from Lemma 4 and from the defini- Proposition 2 (Dynamic Programming Principle). Let α ∈ [0, 1]. For all x ≥ 0 and for all S, T > 0 Proof. Set for brevitỹ Let Y ∈ C + (X(T ), S) be such that Then the map By the arbitrariness of ε, one deduces u α (x, T + S) ≤ũ α (x, T, S). Now, to prove the inverse inequality, let ε > 0 and X ∈ C + (x, T + S) be such that It follows that Thenũ α (x, T, S) ≤ u α (x, T + S) and this concludes the proof.
Proof. Let R > 0. By Lemma 5 there exists a modulus of continuity γ R such that for every In other words, settingφ(·) := u α (·, S) we have thatφ is a locally uniformly continuous map, and by Lemma 1 it attains its global minimum at x = 0. By Proposition 2, for every T > 0 Applying Lemma 4 toū α (x, T ) := u α (x, T + S) and toφ, we deduce that there exists a modulus of continuityν R satisfying Note thatν R depends onφ (hence on u α (·, S)) only via ||φ|| L ∞ ([0,R]) and via the modulus of continuity ofφ. In particular, by Lemma 1, we have ||φ|| L ∞ ([0,R]) ≤ C α 1 (R). On the other hand, the modulus of continuity ofφ is simply γ R , i.e. the modulus of continuity of u α (·, S), which is independent from S by Lemma 5. We then conclude that and, consequently, we deduce the continuity of u α (x, ·) for all x ∈ [0, R]. Since R is arbitrary, the proof is complete.
In particular The proof is postponed in Section 4.1 below. We finally state our main result, whose proof is showed in Section 4.2.

Proof of Theorem 3.
It is left to prove that for all α ∈ [0, 1], v α is a viscosity solution of (25).
Letting ε → 0 + , we finally conclude We deduce, by the arbitrariness of ϕ, that v α is a also viscosity supersolution of (25) and this concludes the proof.

Preliminary results and proof of Theorem 4
We begin with the following regularity result for ψ.

Numerical approximation and simulations
In this section, we introduce a numerical scheme for solving the Hamilton-Jacobi equation (10), and we explore some properties of the numerical solution to confirm the results presented in the previous sections. In particular, we analyze its asymptotic behavior in time, in comparison with the solution of the stationary problem (25). Then, we employ a simple Euler integrator to build optimal trajectories for the control problem (9), and we show the results in different scenarios. For the reader's convenience, we recall here the Hamilton-Jacobi equation (10) together with the boundary condition in space at x = 0, provided by Lemma 1: where the Hamiltonian H : [0, +∞) × R → R, given by is rewritten in a form suitable for the discretization. We approximate the unbounded space-time domain by means of a rectangle [0, x max ] × [0, T ], for sufficiently large real numbers x max > 0 and T > 0, and we introduce a uniform grid with nodes (x i , t n ) = (i∆x, n∆t) for i = 0, . . . , N x and n = 0, . . . , N t , where ∆x = xmax Nx and ∆t = T Nt denote space and time steps respectively, and N x , N t are given integers.
Moreover, we denote the approximations of u, b α , φ on the grid respectively by U n , and we collect them in the vectors U n = (U n 0 , . . . , U n Nx ), B α = (B α,0 , . . . , B α,Nx ), Φ = (Φ 0 , . . . , Φ Nx ). Then, we introduce finite differences for approximating space and time derivatives. In particular, discretization in space requires some care, in order to define a numerical Hamiltonian H which correctly approximates viscosity solutions of the equation. More precisely, using forward/backward differences, we introduce the following two-sided approximation of Du, for i = 1, . . . , N x − 1, and we set which selects the gradient components in an upwind fashion, according to the sign of −B α + 1 2 DU , where (·) + = max{·, 0}, (·) − = min{·, 0} denote respectively the positive and negative parts of their arguments (for further details, we refer the interested reader to [Set99,FF13]).
Finally, we employ a forward difference in time, and we end up with the following explicit time-marching scheme: We remark that the forward difference D R U i is not defined for i = N x , corresponding to the point x = x max at the right boundary. Since the feedback control for the underlying control problem is given by ξ = −Du, we can argue that, for x max sufficiently large (e.g. greater than the endemic population E = N (1 − 1 ρ ) if ρ > 1), the optimal dynamics in (7) pushes to the left, and towards the origin, the trajectory starting from x max . This implies that Du(x max , t) ≥ 0 for all t ∈ (0, T ), and we also have b α (x max ) ≤ 0. We conclude that (−b α (x max ) + 1 2 Du(x max , t)) − = 0, hence we can always select for the node i = N x just the contribution given by the backward approximation D L U i .
We also remark that, due to the definition of b α , the right hand side of the above scheme becomes larger and larger as x max increases. This requires a severe CFL restriction on the discretization steps ∆t and ∆x, in order to preserve stability.
Once the solution has been computed, we reconstruct the following semidiscrete feedback control, by merging the two components of DU and interpolating their values in space, via a linear interpolation operator I: Then, we build the optimal trajectories by integrating (7), using a simple forward Euler scheme y n+1 = y n + ∆t (b α (y n ) + ξ n (y n )) y 0 = x . Now, let us setup the parameters for the numerical experiments. We choose the space domain size x max = 4, the final time T = 5, and N x = 200, N t = 4000 nodes in space and time respectively. Moreover, we choose M (α) ≡ 1, N = 9 4 and γ = 1 as model parameters defining the advection speed b α in (5). On the other hand, α and ρ (and accordingly β = γρ) will be set differently for each test.
We start by choosing the exit cost φ(x) = x and a reproduction factor ρ = 3 2 such that the corresponding endemic population E = 3 4 is a stable equilibrium (attractive in the domain [0, x max ]) for the uncontrolled system, see Proposition 1. In Figure 1, we show the results obtained for α = 1 and α = 1 2 at different times. In the left panels, we report the value functions compared to φ, in the right ones the corresponding optimal controls.
Note that the time horizon T = 5 has been set large enough to reveal the asymptotic behavior in time of the solutions, namely their convergence, up to machine error, to stationary regimes. Similarly, the space boundary x max =4 is large enough to distinguish the growth of the solutions for x → +∞. We observe a linear behavior for the case α = 1, and a quadratic behavior for the case α = 1 2 . This can be better appreciated looking at the corresponding optimal controls, and it is confirmed by the simulation in Figure 2, in which we show, for 0 < α ≤ 1, the asymptotic behavior of the L ∞ norms in space of U and DU at the final time (achieved by monotonicity at x max ) as both x max , T → +∞. In particular, we find out that α = 1 is the only value that produces a globally Lipschitz continuous solution. A rigorous proof of this statement is still under investigation. In Figure 3, we compare some optimal trajectories obtained for α = 1 (top panels) and α = 1 2 (bottom panels). In each plot, we report the endemic population E (dashed line), the uncontrolled/controlled trajectories (bold lines), and the corresponding optimal controls (thin lines). Moreover, we choose two different initial data for the dynamics (7), x = 0.5 and x = 1.25, respectively below and above E. As discussed in the introduction, the fractionary SIS system can be recasted in the model (3), with ordinary derivatives and saturated growth rates. In particular, the growth of infective individuals is softened as α decreases. This effect is apparent in the uncontrolled trajectories. In the same time horizon, we observe that the uncontrolled trajectory approaches E for α = 1, while for α = 1 2 it is "lazier" and still far from the endemic value at the final time. On the other hand, we observe that the optimal control always succeeds in steering the system to the origin (the unstable equilibrium in this case). Nevertheless, while the controlled trajectories are quite similar (as their optimal controls) when the evolution starts from x = 0.5 < E, the case x = 1.25 > E for α = 1 2 requires an additional effort to compensate the slower decay of the corresponding dynamics. Indeed, we observe an optimal control with a larger amplitude in the fragment [0, 0.5] of the time interval.
The numerical results for the case with a reproduction factor ρ ≤ 1 are quite similar to the previous ones, and we omit them for brevity. We just remark that the endemic value now falls out of the space domain (E ≤ 0), while the state x = 0 is a stable equilibrium for the system, see again Proposition 1. This implies that, for all the initial data in (0, x max ], the corresponding uncontrolled trajectories eventually converge to x = 0, whereas the controlled ones have a faster decay, in order to optimize the cost functional (9) for the optimal control problem.
Let us now consider an example with a non smooth exit cost, namely we choose φ(x) = min{2x + 1 2 , 6x 2 }, so to produce a kink in the solution.  Moreover, we choose the parameters as in the previous tests, with the exception of the space domain size, that we set to x max = 2 in order to achieve a sharper CFL condition and mitigate the numerical diffusion of the scheme.
In Figure 4, we show the results for the case ρ = 3 2 and α = 1. In each plot we report, at different times, the value function compared to φ and also the optimal control. We clearly observe that the kink in the solution moves and eventually exits the domain as the time increases. Asymptotically, we obtain a smooth solution as in the previous tests.
We finally consider the case of a smooth exit cost φ(x) = x+exp(−40(x− 1 2 ) 2 ), which corresponds to penalize the final distribution of infective individuals around the point x = 1 2 more than for larger values (up to about x = 3 2 ). The results for the case ρ = 3 2 and α = 1 are reported in Figure 5. We observe that, in the first part of the evolution, the point x = 1 2 acts as a barrier, preventing some states of the system to be steered to the desired one x = 0. More precisely, the local minimizer of φ (around about x = 0.8) is more favorable for states beyond this barrier, where the optimal control has a change of sign. This creates a kink in the solution, which starts moving towards the right boundary of the domain only at a later time.
In Figure 6, we compare the corresponding optimal trajectories, obtained for the initial data x = 0.48 and x = 0.52, respectively slightly below and above the barrier. In the first case, we obtain a controlled trajectory similar to the previous tests, with just a larger amplitude in the control due to the choice of φ. On the other hand, the second case confirms the scenario discussed above. Indeed, the optimal control acts in the positive direction for a small amount of time, pushing the controlled trajectory close to the endemic population, then readily jumps to a negative value, and starts steering the system to the origin. We computev α on our numerical grid, approximating the integral by a simple trapezoidal quadrature rule, using the same space step ∆x. Moreover, we set x max = T = 10, and we choose the same exit cost φ of the previous test (note that this affects the convergence of u in time, whilev α only depends on φ(0)).
In Table 1, we report the results of the comparison under grid refinement, for different choices of α and ρ, evaluating the difference u(·, T ) −v α (·) both  in L ∞ and L 2 space norms. As ∆x → 0, we clearly observe a decay of the errors, respectively of order O(∆x) and O(∆x 2 ) for the two norms, and also a slowdown in convergence as α and ρ decrease. This numerical experiment is in agreement with the result proved in Theorem 4. In particular, convergence is obtained on bounded space intervals, and we have observed in all the experiments that possible irregularities of u are pushed out of the domain towards infinity, before approaching the stationary smooth solutionv α in the limit T → +∞.  Table 1. Comparison between u andv α under grid refinement for different model parameters.