{'fraction_non_alphanumeric': 0.0463615215753088, 'fraction_numerical': 0.009305068641263393, 'mean_word_length': 4.437733832175307, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, '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}}
Optimal policy design to mitigate epidemics on networks using an SIS model

Understanding how to effectively control an epidemic spreading on a network is a problem of paramount importance for the scientific community. The ongoing COVID-19 pandemic has highlighted the need for policies that mitigate the spread, without relying on pharmaceutical interventions, that is, without the medical assurance of the recovery process. These policies typically entail lockdowns and mobility restrictions, having thus nonnegligible socio-economic consequences for the population. In this paper, we focus on the problem of finding the optimum policies that"flatten the epidemic curve"while limiting the negative consequences for the society, and formulate it as a nonlinear control problem over a finite prediction horizon. We utilize the model predictive control theory to design a strategy to effectively control the disease, balancing safety and normalcy. An explicit formalization of the control scheme is provided for the susceptible--infected--susceptible epidemic model over a network. Its performance and flexibility are demonstrated by means of numerical simulations.


I. INTRODUCTION
The ongoing COVID-19 pandemic has highlighted the key role played by public health authorities in enacting nonpharmaceutical interventions (NPIs) to "flatten the epidemic curve" when no effective pharmaceutical treatments such as vaccines are available [1], [2]. However, NPIs typically entail the implementation of harsh measures, including lockdowns and restrictions of personal freedom of movement, which may yield severe socio-psychological and economic consequences [3]. Thus, they should be implemented keeping a reasonable balance between safety and normalcy. To this aim, the development of accurate tools to predict the course of an epidemic and evaluate the impact of different interventions has become a task of paramount importance for the scientific community, aiming at assisting public health authorities in their decisions.
The mathematical modeling of epidemics has emerged as a valuable framework to perform such a task [4], [5], [6], [7], [8]. Relevant examples can be found in the useful insights provided into the ongoing COVID-19 pandemic [9], [10], [11], [12], [13], [14], [15], [16], [17], [18]. Within this framework, network models have become popular as they allow us to capture the relation between human mobility and the spatial spread of a disease. Of particular interest is the problem of understanding how to effectively control the spread of an epidemic disease on a network by acting on the nodal dynamics and on the network structure.
In the literature, such problems have often been addressed assuming limited changes in the network structures, that is, by studying how to re-arrange the network structure or distribute antidote in order to increase the population's resistance against a possible epidemic. Important results have been found by using geometric programming [19]; distributed algorithms have been recently proposed to address this problem [20], [21]. However, the ongoing health crisis has highlighted the importance of having control schemes, which take into account the dynamic evolution of the epidemic process, and can thus be updated online, as the outbreak evolves. While a considerable body of literature has been proposed to study dynamical control strategies for vaccination campaigns and antidote distribution using optimal control theory [22], [23] and Model Predictive Control (MPC) [24], [25], [26], limited results are available in the absence of effective treatments, that is, when the control action has to focus on the contagion mechanisms rather than on the recovery mechanism. Recently, motivated by the ongoing COVID-19 pandemic, some efforts have been devoted to bridging this gap by proposing feedback control interventions and leveraging MPC [13], [14], [16], [18].
Here, inspired by these works, we propose an optimal control approach in order to address the problem of mitigating an epidemic process spreading by means of regional policies that entail activity reductions and targeted mobility restrictions. This holistic approach to NPIs constitutes a key contribution of this paper and differs from other works in the literature that typically focus on optimizing a specific intervention policy (e.g., social distancing in [27]) only. Specifically, we consider a discrete-time deterministic Susceptible-Infected-Susceptible (SIS) epidemic model on a network in which two types of control actions are included to mitigate the contagion: actions to reduce the social activity in some regions of the network (modeling, for instance, lockdown measures, closures of activities, and curfews), and policies to limit or ban travels between specific nodes (e.g. cities, provinces, regions, countries). Then, an optimal control problem is formulated to find a strategy that mitigates the spread of the disease in the network, limiting the negative consequences of NPIs. The proposed optimal control problem takes into account important features, such as the balance between safety and normalcy, the need of keeping the epidemics under control, and the increase of socioeconomic costs associated to the implementation of NPIs, caused by accumulating factors. Different from many control schemes proposed in the literature (see for example [4], [7] for an overview), our formalization does not (necessarily) have the objective to eradicate the disease (which may be extremely costly -and practically unfeasible-using only NPIs). In contrast, it allows the controller to set an acceptable prevalence of the disease (which may depend on the hospital capacity and may vary across the nodes), providing thus a framework that might realistically be adopted to assist public health authorities in their decision toward mitigating epidemic outbreaks. Finally, we propose several simulations in which we highlight the performance of the proposed control strategy, the benefits and disadvantages of farsighted policies in comparison to myopic ones, and perform a sensitivity analysis on how the objectives of the policy maker affect the final solution.
The rest of the paper is organized as follows. In Section II, we present the SIS epidemic model. In Section III, we discuss the control policies. In Section IV, we propose our optimal control policy. In Section V, we illustrate the results of our numerical simulations. Section VI concludes the paper and outlines the avenues of future research.

A. Notation
We gather here the notation used throughout the paper. The set of nonnegative integer numbers is denoted by N, and the set of nonnegative real numbers is denoted by R + . Given a vector x, we denote by x ⊤ its transpose and, for a positive definite matrix S ≻ 0, we denote by x S := √ x ⊤ Sx the norm of x weighted by S. With 0 and 1, we denote the all-0 and all-1 vectors, respectively. Given a matrix S, we denote by σ(S) its maximum singular value. S ij and S i denote the element in the i-th row and j-th column and the i-th row of S, respectively. We use the following notation col ((x i ) i∈ [1,..

II. SIS EPIDEMICS MODEL OVER NETWORKS
In this section, we introduce the SIS epidemic model, illustrate its main properties, and discuss the assumptions needed to formulate the control problem.
We consider a discrete-time deterministic SIS epidemic model on a network [6]. The network structure arises from the local interactions among several nodes (communities), which represent geographical entities, such as countries, regions, or even cities, thanks to the high flexibility of the model. We assume there are N communities interconnected by P links, which can represent for instance roads, air routes or even simple geographic adjacencies. In general, a link between two communities means that there is a flow of people between them. The network is hereafter formalized via an undirected and connected graph G, where the communities correspond to the set V := {1, . . . , N } of the nodes. A connection between two communities i, j ∈ V is denoted via an edge e ℓ connecting the two corresponding nodes, defined as the unordered couple e ℓ := {i, j}. We assume that all the self-loops {i, i}, i ∈ V are present. The edge set is the collection of all the edges of the graph, i.e., E := {e 1 , . . . , e P }. Therefore, the graph is defined as G := (V, E). The fraction of infected individuals in community i ∈ V at time t ∈ N is denoted by x i (t) ∈ [0, 1]. This quantity describes the temporal evolution of the health state of community i. We assume that the health state evolves according to the dynamics of a discrete-time SIS model [6], i.e., for all time-steps t ∈ N, we have where x + i := x i (t + 1) and x i := x i (t). Next, we discuss the role and physical interpretation of each parameter appearing in (1).
1) Recovery rate µ ∈ [0, 1]: It is the rate at which the individuals manage to recover from the disease. Here, we assume that the recovery rate is constant across the entire population and cannot be increased by the controller, capturing those epidemics for which a cure is not yet developed, e.g., the early spread of COVID-19.
Assumption 1 (Positive recovery rate): For all i ∈ V, the recovery rate µ > 0 is constant and strictly positive.
2) Infection rate β i > 0: This represents the rate at which the individuals in community i ∈ V become infected when they get in touch with others. The higher this value is, the easier people become infected. Such a value can differ among the communities, e.g., due to the implementation of different NPIs. Infection rates are gathered in the n-dimensional vector β. The trajectory of x(t) := col ((x i (t)) i∈V ) can naturally have two possible behaviors, depending on the model parameters. Either it converges to the disease-free equilibrium, i.e., x(∞) = 0, or the disease-free equilibrium becomes unstable, and the trajectory converges to a (unique) endemic equilibrium with x i (∞) > 0, for all i ∈ V [28], [29]. The threshold between these two regimes depends on whether σ(diag (β)Ā)/µ is smaller or greater than 1. A simplernetwork-independent-sufficient (but nonnecessary) condition for the trajectory to converge to an endemic equilibrium can be established by requiring that β i /µ > 1 for all i ∈ V. In this work, we focus on the case in which the disease does not die out naturally, and thus NPIs have to be put in place toward mitigating the epidemic outbreak. Hence, we make the following assumption.
Assumption 2 (Disease spreading): For all i ∈ V, the infection rate satisfies β i > µ.
3) Communities interaction A i,j ∈ [0, 1]: People in j ∈ V can move to the neighboring communities, denoted by N j := {ℓ ∈ V : {j, ℓ} ∈ E}, and interact with people there. The population that flows from j to i ∈ N j influences the health state of community i, i.e., x i . This is modeled in (1) via the weighted adjacency matrix A : where A ij is the weight associated to the edge (i, j) of G. The diagonals of A represent the part of population that remains in the same community. We now introduce the following assumption on this matrix. Fig. 1: Schematic of the controlled network SIS model. The timevarying weighted adjacency matrix A(t) determines the interaction that community i has with other communities. The time-varying infection rate βi(t) determines the new infections, while the constant recovery rate µ governs the recovery process.

Assumption 3 (Stochastic and positive diagonals):
The weighted adjacency matrix A associated with the undirected graph G is stochastic, i.e., A1 = 1, and with strictly positive diagonal entries, i.e., A ii > 0, for all i ∈ V.
In the following section, we discuss the control policies.

III. CONTROL POLICIES
To mitigate the spread of a disease, a policy maker can apply endogenous or exogenous NPIs to each local community. The former category includes all those measures put in place inside the local community i, e.g., lock-downs, usage of face masks or encouraging social distancing. The latter instead concerns those measures that limit the inflow of people from neighboring communities, i.e., implementation of travel bans or requesting quarantine upon entrance. It is clear that these policies directly affect the dynamics in (1). Specifically, for each community i, the endogenous measures reduce the infection rate β i , while the exogenous ones act on A. The static values of β i and A refer to the uncontrolled evolution of the epidemic, while the intervention of the policy makers transforms these parameters into time-varying functions β i (t) and A(t), which can be seen as inputs to be designed in a (possibly) optimal way to control the epidemic spreading. Thus, the dynamics of the controlled evolution of the SIS epidemic become as illustrated in the schematic in Fig. 1.
If the policies put in place are effective, the value of the infection rate should decrease, i.e., We consider policies that affect the uncontrolled coefficients in (1) linearly. So, we define time-varying functions β i (t) and A(t) in (2) as where v i (t) ∈ [0, β i ]. Furthermore, it is reasonable to assume that the policy cannot implement new connections among the community, this translates into the following set of constraints on W (t) : From (3b) and (4), it follows that the time-varying weighted adjacency matrix A(t) defines a directed subgraph G(t) ⊆ G, which is composed of the same nodes and a subset of the edges of G, i.e., G(t) := (V, E(t)) with E(t) ⊆ E, so it may not be connected at every t ∈ N.

IV. OPTIMAL CONTROL POLICY
In this section, we present the main result of this paper on optimal control policies.

A. Problem formulation
The adoption of stringent policies carries costs including both monetary (e.g. recession) and social (e.g. personal restrictions) aspects [3]. In this section, we formalize the problem of choosing the optimal NPIs that minimize the cost while keeping under control the epidemic. Therefore, at each time instant t ∈ N, we have to design the values of the control actions W (t) and v i (t) for all time instants k in the prediction horizon T := {t, t + 1, . . . , t + T h }, the length of which (T h ∈ N) may vary due to the specific epidemic.
In the remainder, we use the index k when we refer to an instant belonging to T and t otherwise.
Each population i ∈ V is assumed to establish a desired trajectoryx i (k), for all k ∈ T \ {t}, that is the fraction of infected population x i (k) considered to be acceptable at time k. An interesting case is the one in which the value ofx i (k) is greater than 0 at the beginning of T , and it decreases over time, i.e.,x(k) ≥x(k+1) for all k, k+1 ∈ T \{t}. Its slope suggests how aggressive the desired community's policies should be. Notice that the terminal valuex i (t + T h + 1) is not necessarily constrained to be 0. In fact, it is reasonable that some communities accept a small fraction of infected in exchange for relaxed NPIs.
Next, let us define the vector of all policies put in place at time instant t by agent i as the (N + 1)-dimensional vector u i (t) := col (v i (t), W i1 (t), . . . , W iN (t)). The cost for each community i consists of two antagonizing components. The first is the health-care cost c HC i (x i (t);x i (t)) due to the presence of more infected than the desired onesx i (t), which is individually decided by each community and communicated to the health authority. If the optimal policies leads to x i (t) ≤x i (t), then the cost is 0; otherwise it is assumed quadratic in the difference with respect to the desired value, and thus it reads as 0 being a (possibly time-varying) weight. On the other hand, we consider a (quadratic) control cost associated with the implementation of control policies. So, the global cost that the N communities face over the prediction horizon reads as where S i (k) ≻ 0 is the diagonal matrix of the N + 1 weights associated to the NPIs u i (k). The fraction of infected and the adopted policies of all the populations at t are x(t) and u(t) := col ((u i (t)) i∈V ), respectively. The pair x = col ((x(t + 1)) t∈T ) and u = col ((u(t)) t∈T ) denotes in compact form all the variables involved. In our general formulation, the weights associated to the health-care cost (i.e., q i (t)) and to the control cost (i.e., S i (t)) are timevarying, since the same fraction of infected individuals or the same level of NPIs may yield different costs depending on the timing, for example due to the increasing hospital preparedness or the accumulation of socio-economic costs. It is worth noticing that only 2P − 2N values of W (k) can be freely designed at each time instant k ∈ T . In fact, for each one of the P −N off-diagonal edges, we can set W ij (k) and W ji (k); then the weight associated with the self-loop is constrained by (4a). Thus, the number of elements of u(t) to be optimized is 2P − N . Finally, we cast the problem of designing the optimal policies u * for the control dynamics (2) where This problem belongs to the class of Nonlinear Model Predictive Control (NL-MPC), due to the highly nonlinear and nonconvex controlled SIS dynamics in (2). It is well known that these problems are hard to solve in their original form. Nevertheless, we establish that a feasible solution to (P) always exists, and thus that the problem is worth to be studied.

Proposition 1 (Solution existence):
For every initial condition x(t) ∈ [0, 1] N , there exists at least one pair (x * , u * ) that is a solution to (P).
Proof: From the constraints' definition (2) and (4), it follows that the pair u = 0 and x, obtained via (1), always lies in Ω, hence Ω = ∅. The box constraints on the control inputs and the fact that the dynamics (2) are positively invariant with respect to [0, 1] N , see [28,Lem. 1], implying that Ω is bounded. For every constraint, the set of points satisfying it is closed. Therefore, Ω is the intersection of closed sets, and hence it is closed itself. Since J is continuous, we invoke the Weierstrass theorem to conclude the existence of (x * , u * ) ∈ Ω that minimizes J, and, in turns, solves (P).

B. NL-MPC solution algorithm
In the literature, there are several approaches to solve NL-MPC. Specifically, in order to solve nonlinear constrained optimization problems with differentiable cost functions, the most popular approaches are based on the Sequential Quadratic Programming (SQP). As shown in [30,Ch. 18], [31] and reference there in, these methods provide excellent convergence properties and ensure a fast convergence to a (local) optimum of the original nonlinear problem. The SQP is an iterative algorithm in which, during each iteration p, a candidate optimal trajectory (x p , u p ) is computed as the solution of a Quadratic Programming (QP), obtained by linearizing the constraints and approximating the cost via a quadratic funciton, see [30,Alg. 18.3]. The linearization is performed with respect to the trajectory (x p−1 ,ũ p−1 ) defined as where (x p−1 , u p−1 ) is the solution of the QP solved at the previous iteration, and α p−1 can be computed for example via line search as in [30,Eq. 18.28]. If the candidate solution satisfies some convergence conditions then (x * , u * ) = (x p , u p ); otherwise a new iteration is performed.
The sole nonlinearity in the constraints of (P) is associated with the dynamics in (2). The effectiveness of the NPIs depends on the predictive accuracy of the linearized dynamics. This difference is studied via numerical simulations in Section V, where it is shown that the linearized dynamics allow to obtain a very high prediction accuracy for a sufficiently long prediction horizon. Thus, the used NL-MPC algorithm generates an effective solution for the control problem.
Finally, in Figure 2 we depict the complete control scheme to solve the problem of designing optimal NPIs to control an SIS type dynamics. Specifically, at each time instant t an instance of (P) is cast. The optimal trajectory (x * , u * ) is computed over T via SQP. Finally, a receding horizon approach is implemented by applying only NPIs associated to the first instant, i.e., u * (t), and then the loop starts again. This implementation allows us to minimise the prediction error inherently present in the prediction of the model.

V. SIMULATIONS
In this section, we present several numerical simulations that validate the procedure proposed to design optimal NPIs and provide insightful information on how different parameter choices (e.g., the weights on the control and the length of the prediction horizon) influence the final optimal policies.

A. Reference tracking and reproduction number
The simulations are performed over a randomly generated connected network with N = 4 that is described by a weighted adjacency matrix with off-diagonal nonzero entries chosen from a normal distribution with mean 0.1 and variance of 0.2, and diagonal entries such thatĀ1 = 1, obtaining  Figure 2 (red solid curve) for all i ∈ V. In (b), we plot the amount of endogenous (blue) and exogenous (green) control applied to every community i ∈ V to obtain xi(t), respectively βi(t) and Wi(t) 2 2 .
The uncontrolled infection rate β i of each node is selected uniformly and randomly in the interval [0.3, 0.6], obtaining β = [0.3, 0.59, 0.3 , 0.45]. The recovery rate is chosen to be µ = 0.15. So, the uncontrolled dynamics converge over time to an endemic equilibrium [28]. Note that, since µ is the inverse of the time for recovering, if one time-step corresponds to a day, then this setup may represent a disease like seasonal flu, for which (approximately) 6 days are enough to recover.
On the other hand, if one time-step is interpreted as one week, the setup may capture diseases with longer recovery time, e.g., gonorrhea. The initial fraction of infected for each population is set as x(0) = [0.65, 0.55, 0.75, 0.40], that is, a scenario of an endemic disease. We denote the trajectory of the fraction of infected individuals in the uncontrolled dynamics (1) by x uc (t). The desired reference trajectorŷ x i (t) of each community is assumed to start close to x i (0) and linearly decrease until it reaches its terminal value at t = 20 that isx(20) = [0.1168, 0.0548, 0.0856, 0.1175].
The length of the prediction horizon is set equal to T h = 10 and we will discuss its optimal value in Section V-B. The weights q(t) and S(t) are chosen constant in time, uniform across the different nodes, and equal to q i (t) = 1, for all i ∈ V, while S is so that all off-diagonal terms are equal to 0, the diagonal terms that correspond to entries of v are equal to 0.2 and those that correspond to entries of W are equal to 0.05. Note that the weights are such that applying the same amount of exogenous and endogenous control respectively have similar costs.
In Figure 3a, we present the trajectory of x i (t), obtained from (2), by following the control scheme in Figure 2 (red solid curves), compared with the value of the uncontrolled SIS dynamics (red dashed curves), and the reference (blue solid curve). As expected from having a higher weight on the healthcare cost, the value of x i (t) remains relatively close to the desired onex i (t). From Figure 3b, it is clear that, in this scenario, acting on local restrictions as lock-down and social distancing is more effective than implementing travel bans. This is consistent with empirical observations during the ongoing COVID-19 pandemic, suggesting that travel bans are more effective in the early stages of an epidemic outbreak, i.e., when x(0) is close it 0 [15]. In the first instances, the control action is more intense due to the decreasing reference and, after t = 20, it is relaxed. This can also be assessed via the ratio β i /µ depicted in Figure 4a. In fact, during the initial time instants, all this ratios decrease reaching values lower than 1. This suggests that, by adopting the NPIs implemented in the first phase for the entire duration of the simulation, the trajectory would converge to the disease-free equilibrium, at the cost of dramatically increasing the social and economical impact of the control.
In Figure 4a, we also report the evolution of the singular value σ(B(t)A(t))/µ where B(t) := diag (β(t)), which is known to determine the behavior of the epidemic process [5],  [28]. It is also strictly related to the effective reproduction number R t , which can be obtained by multiplying such a quantity by the fraction of susceptible individuals in the population. As expected, the value at which σ(B(t)A(t))/µ settles is above 1, since the terminal value of the referencê x(t) is greater than the disease-free equilibrium. It is interesting to notice that the optimal control W (t) tends to isolate the communities by increasing the weight of each self-loop and decreasing the off-diagonal entries. Analyzing more complex graphs may lead to the emergence of different patterns in the exogenous control actions, and this is left for future research.
In Figure 4b, the accuracy of the linearized model used to solve (P) at time t = 4 is shown. By solving the optimization problem, we obtain the trajectories x * (t) and u * (t) for t ∈ T = {4, . . . , 14}, which we denote by x l (t). We compute x nl (t) by applying the whole u * to (2) and letting it evolve in an open loop. The estimation of the nonlinear behavior of the dynamics is almost perfect during the first six time instants of T , while it becomes less accurate during the second half. An even longer prediction horizon would exacerbate this difference affecting the performance, as we show in the next section.
Overall, the policies computed by the proposed control scheme generate noteworthy effects in controlling the epidemics by achieving the two main objectives of the proposed approach, i.e., an acceptable level of infected and at a low social and economical price. Nevertheless, it is evident that the optimal solution to the problem (P) highly depends on the choice of the weights q i (t) and S i (t) in the cost function J(x, u). To better understand the effect of the endogenous and exogenous measures on the (steady-state) health state of the overall population, we have performed a sensitivity analysis of the weights applied to the control action. For all i and t, we set q i (t) = 1 and S i (t) equal to the block diagonal matrix that has on the diagonal s v and sw N I, which are the weights on v i (t) and W i (t), respectively. In order to compare the performance obtained with different weights, we introduce the following index: As shown in Figure 5, we obtain larger values of π (i.e., the system converges close to the endemic equilibrium x uc ) when the weights of the control actions are large, while smaller values of π (i.e., the system converges close to the desired referencex) arise if s v and s w are small. Furthermore, from Figure 5 one can observe that changing the weight s v makes the index π change more significantly than changing the weight of the exogenous actions. Then, we can conclude that the endogenous actions are more effective than the exogenous ones. However, we can also observe that when the cost of implementing for instance lock-downs is very high (i.e., larger values of s v ), putting in place travel bans is very beneficial.

B. Myopic vs predictive policies
Next, we consider the problem of choosing the optimal length for the prediction horizon. This is a critical choice that a policy maker has to perform. In fact, if T h is too small the policies will be myopic failing to prepare in time for the future evolution of the epidemics. On the other hand, longer prediction horizons will require an increased computational effort. Moreover, due to the difference between the real dynamics and the ones used in (P), a long prediction horizon may lead to inaccurate estimation and consequently to incorrect precautionary policies, see Figure 4b. Therefore, we believe that there is an optimal length for the prediction horizon.
To validate this claim via simulations, we generate 50 different scenarios in which A, β andx are randomly chosen as described in the previous section, while the rest of the coefficients of the problem do not change. As index of the performance we consider the sum over the whole simulation of the cost actually experienced by the communities, we denote it simply by J. Figure 6a clearly depicts a Pareto front in which, after an initial reduction of J due to the growing T h , there is a diminishing return. Moreover, for even greater values of T h the performance starts to slowly get worse. The computational time of each one of the simulation grows linearly with the value of T h , as shown in Figure 6b, this can be a serious issue in the case of large scale networks. Therefore, the optimal value of T h for the simulated cases is in the range 5 − 7. In fact, it is large enough to benefit from the initial steep increment in terms of performance but not that large to make the estimations less accurate.

VI. CONCLUSION AND OUTLOOK
In this paper, the problem of designing optimal NPIs to ensure containing an epidemic on a network while minimizing the costs due to restrictive policies has been cast as a nonlinear constrained optimization problem over a prediction horizon. We have proposed a solution via a linearization technique and a receding horizon strategy, in which the real evolution of the disease is used as a feedback to cast a new instance of the optimization problem and compute the upcoming policies that have to be implemented.
The performance of the proposed control scheme has been demonstrated via preliminary numerical simulations. Specifically, the good accuracy of the linearized dynamics for sufficiently long time horizons guarantees accurate farsighted forecasts, avoiding thus the risks of relying on myopic policies. Our numerical findings have shown nontrivial solutions for the optimal control strategies, depending on the costs associated with the implementation of the different types of NPIs and on the requirements on the desired evolution of the epidemic outbreak. We finally would like to stress that, while the current formulation has been developed using an SIS epidemic model, the same control scheme can be extended to more complex epidemic models, such as those used to capture the features of the ongoing COVID-19 pandemic [9]. Such an extension, and the analysis of the proposed control scheme in real-world scenarios is an important direction for the future research.