On the dynamics of a viral marketing model with optimal control using indirect and direct methods

The complexity of optimal control problems requires the use of numerical methods to compute control and optimal state trajectories for a dynamical system, aiming to optimize a particular performance index. Considering a real viral advertisement, this article compares the dynamics of a viral marketing epidemic model with optimal control under different cost scenarios and from two perspectives: using numerical methods based on the Pontryagin’s Maximum Principle (indirect methods) and methods that treat the optimal control problem as a nonlinear constrained optimization problem (direct methods). Based on the trade–off between the maximization of information spreading and the minimization of the costs associated with it, an optimal control problem is formulated and studied. The existence and uniqueness of the solution are proved. Our results show not only that the cost of implementing control policies is a crucial parameter for the spreading of marketing messages, but also that low investment costs in control strategies fulfill the proposed trade–off without compromising the financial capacity of a company.


Introduction
Viral Marketing (VM) is a marketing strategy that takes advantage of network effects to fostering customers to promote information product exchange among their social network [1].
The impact of network effects on consumers can be quite erratic, and marketing companies cannot always guarantee the success of an advertising campaign.In this regard, research studies have been applied epidemiological models beyond their conventional use, particularly on studying the dissemination of information content within a certain population (see, e.g., [2,3,4,5] and references cited therein).Other investigations have developed models to study viral propagation of memes, via social communication, on a target population [6], spreading schemes for VM within social networks [7], algorithms that minimize the number of initial seeds to generate maximum awareness of products [8], or even approaches that try to better identify the seeds that can broadly spread a marketing message [9].Since companies generally face trade-off situations related to the conception of strategies that maximize the spreading of information among consumers in a cost-effective way, it is relevant to conceive mathematical techniques able to model these particular features.In this regard, in addition to applications in the area of health (e.g., Human Respiratory Syncytial Virus (HRSV) [10] or Dengue [11]), Optimal Control Theory [12] has been explored to assess the impact of control measures on information epidemics (see, e.g., [3,4,13,14]).

SIR Epidemiological Model and Properties
Proposed in [16], the classical SIR model subdivides the population into three mutually-exclusive compartments: susceptible individuals (S), infected individuals (I) and recovered individuals (R).This model is represented by the following system of ordinary differential equations with initial values: However, these compartments can be described under a marketing context.In fact, the marketing definition of these compartments and model parameters was already discussed in a similar context (see, e.g., [4,5,6]).The state variables S(t), I(t) and R(t) thus represent, respectively, the target population who can spread the marketing message, the individuals who spread it into their social networks and the ones who stop diffusing it, over time t.Regarding the model parameters, β and γ traduce the predisposition to share the marketing message by targeted individuals and the interruption of the spreading process by infected individuals, respectively.Note that, when the marketing message is appealing, individuals in the class S move to the class I at a rate β.Naturally, over time, individuals in the class I stop sharing the marketing message and move to the class R at a rate γ.Moreover, the total population size, N , is assumed to be constant over time t, that is, N = S(t) + I(t) + R(t).The value N was estimated based on the number of YouTube users in 2013 [17].
After modeling the system (1) according to the maximum peak of infection reached in the campaign's data referred in [15] (first seven days of the campaign), the parameters β = 67.5244and γ = 65.0751 were estimated based on the FMINSEARCH function that implements the Nelder-Mead algorithm (see [18]), from MATLAB optimization toolbox.Additionally, considering that at the beginning of the campaign almost everyone is susceptible to have contact with the marketing message, the basic reproduction number for system (1) is given by message is not widespread.On the contrary, if R 0 > 1, the message is widely spread within the target population (see, e.g., [19] for more details on the dynamics of R 0 ).In this regard, for the estimated parameters β and γ, we get R 0 > 1 -which corroborates that Dove's campaign was a viral marketing epidemic.

Optimal Control Problem
Based on the SIR model (1), let us now consider the fractions s(t), i(t) and r(t) such that s(t) = S(t) N , i(t) = and t f is the considered campaign deadline.Considering these fractions, the system (1) is mathematically and epidemiologically well posed (see [19] and the references cited therein), ensuring the boundedness of state trajectories.
To the system (1), one control function u(•) is now added, resulting in a controlled SIR epidemiological model given by the following system of ordinary differential equations: ( The rate of change of the total population is given by ds(t) dt + di(t) dt + dr(t) dt = 0 and the control function u(•) represents the recruitment of susceptible individuals to act as spreaders within social networks (e.g., through advertisements in mass media [3]) by marketing companies, purposing to maximize the levels of information diffusion.The admissible control set Ω is defined by For the bounded control u(•), the application of the control policy is maximum if u = u max and null if u = 0.The main goal is to find the optimal value u * for the control u, in such a way that the state trajectories s, i and r are the solution of the system (2) over [0, t f ], and maximize the objective functional (3).At this point, following the work of [3], the optimal control problem consists in maximizing the fraction of individuals who had contact with the marketing message at t f , by minimizing the intervention costs related to the implementation of the control strategy u, i.e., subject to (2), where the nonnegative constant B represent the weight of the investment costs associated with the control strategy u.Under a marketing perspective and assuming that the costs do not take a linear form, we considered a quadratic structure for the weighted control function instead of a linear form as in [4].Henceforth, let s * (t), i * (t), r * (t) and λ * 1 (t), λ * 2 (t), λ * 3 (t) denote the optimal state and adjoint functions, respectively.

Proof
This proof is based on the Cesari Theorem [20].Hence, the following conditions must be satisfied: 1.The set of solutions of the controlled system with initial conditions (2) with u ∈ Ω is not empty.2. The set Ω is closed and convex.3. The right hand side of the state system (2) is continuous, bounded above by a sum of the bounded function u and the state variables, and can be written as a linear function in u with state variables as coefficients.4. The integrand of the objective functional ( 3) is concave on Ω with respect to u and there exists constants Firstly, preserving the initial conditions and noting that r(t) = 1 − s(t) − i(t), the system (2) can be reduced to The Condition 2 is satisfied since Ω is closed and convex by definition.Condition 3 is satisfied as follows (see, e.g., [21,22] for similar arguments): The reduced system (4) can be rewritten as where the matrix A, the state vector x and the nonlinear function B(x) are defined as: In addition, notice that where by taking is uniformly Lipschitz continuous.Hence, also attending to the definition of Ω and the non-negativity of the state solutions, Condition 1 is satisfied.Since the state variables s(t), i(t) and r(t) and the cost functional (3) take values in a compact interval [0, 1], and 0 u(t) u max , ∀t ∈ [0, t f ], then V (x) is bounded above by a sum of the bounded function u and the state variables.Moreover, it is easy to show that V (x) can be written as a linear function in u(t) with state variables as coefficients, thus proving the Condition 3. Finally, Condition 4 is satisfied since the integrand of the functional (3) is concave on Ω and According to the Pontryagin's Maximum Principle (PMP) [12] and considering the optimal control problem (2) and (3), if u * (•) is a control that is optimal for the problem, then there exists a Lipschitz continuous map called covector and where the function H defined by is called the Hamiltonian.
Based on the above, it is possible to derive the optimality system, which consists on the state system with initial conditions (2) coupled with the adjoint system (11) and transversality conditions (12) together with the characterization (13).
Additionally, since the state and adjoint functions are bounded and the systems ( 2) and ( 11) preserve the Lipschitz structure, u * is unique for t f sufficiently small (see, e.g., [23]).However, the uniqueness of optimality system holds for any value of t f since the state system (2) is autonomous.
The optimal control and state can be obtained by solving the optimality system using numerical methods for optimal control problems.

Numerical Implementation and Methods
To compute the solution of the optimality system and maximize the objective functional (3), two different implementations are considered: indirect and direct methods.Under Dove's campaign real data [15], these implementations were compared in terms of model dynamics for the first seven days (t ∈ [0, 6]) using both the parameter values β and γ estimated in Section 2 and the initial conditions Indirect methods to solve optimal control problems are based on the Pontryagin's Maximum Principle to compute the optimal solution.Therefore, these methods require the state system with initial values (2), and necessary conditions (10), (11), (12) to find the solution of a given optimal control problem, by converting it into a boundary value problem.On the other hand, direct methods use discretization with respect to time to solve optimal control problems, in which the cost functional is directly optimized [24], treating the optimal control problem as a nonlinear optimization problem (NLP) (see, e.g., [25]).Some advantages and drawbacks can be pointed out to both indirect and direct methods.According to Trélat [26], direct methods are more robust and less sensitive to the choice of the initial conditions than indirect methods, being more easier to implement.However, in comparison with indirect methods, reaching to a desirable precision is not so easy when direct methods are employed.In addition, the author states not only the possibility of obtain local minima when the direct discretization of an optimal control problem is employed, but also the necessity of a large amount of memory, which in turn can lead to inefficiencies when, for instance, a large dimension problem is considered.On the other hand, indirect methods provide high levels of numerical accuracy, but their implementation can be quite difficult due to the necessity of computing derivatives and necessary conditions related to the PMP.
Based on the foregoing, in order to assess the dynamics of the model using both indirect and direct methods, three distinct weight factors B = 10 −1 , B = 10 and B = 10 2 were considered to generate numerical simulations.These factors traduce, respectively, a low, moderate and high investment cost scenario in recruiting susceptible population to act as spreaders on the target population, by marketers.The choice of these values is related to the formulation of investment cost scenarios that best describe and illustrate the marketing reality.
At this point, as an indirect numerical approach, Forward-Backward Sweep method (FBS) is considered.Generally, this method uses a fourth order Runge-Kutta scheme to solve forward and backward in time the state and adjoint equations in compliance with the differential equations in the optimality system (see, e.g.[27] for more detailed information).
In contrast, to use direct methods, the problem (3) must be discretized with respect to time t.For that, a first order Euler's method is considered where the time t = kh, for k ∈ {0, 1, ..., M } moves forward in uniform steps of length h and M is the the last index of the Euler's scheme, corresponding to t f = 6.Note that higher order discretization schemes can also be used despite of no notorious advantages arise [28].
Let f : [0, 1] 3 → R 3 be a C 1 map.According to the Euler's scheme, considering where the update is given by The approximation x n+1 of x(t) at the point t n+1 has an error committed proportional to h 2 .Moreover, it is considered h = 0.006 since it gives a proper balance between accuracy and computational efficiency.
KNITRO, IPOPT and LOQO were used as software packages to solve the discretized problem above.A short description of these packages is given as follows.
KNITRO [29] is a nonlinear solver for large-scale linear, quadratic, and general smooth nonlinear optimization problems with several variables, continuous and integer.For continuous problems, several algorithms are available such as: direct interior-point, conjugate gradient interior-point and sequential linear-quadratic active-set algorithms.
IPOPT [30], for Interior Point OPTimizer, is a software package for large-scale nonlinear optimization, written in FORTRAN and C, that implements a primal-dual interior-point method using line search methodologies based on filter methods.LOQO [31], is an optimization solver for smooth constrained problems, supported to linear, quadratic, nonlinear, convex or nonconvex objectives and constraints in continuous variables.
Using these solvers, the discretized problem (17) was modeled in the algebraic modeling language AMPL [32] and numerically solved in NEOS Server platform [33].The indirect method was implemented and solved in MATLAB.It should be noted that, within the direct method, the CPU time can vary depending on the machine used in NEOS server platform.Thus, this variable is not taken into account when analyzing the results.

Numerical Results and Discussion
This section reports and compares the computational results not only regarding the state and control variables of the state system (2), but also for the cost functional (3), using both the indirect and direct numerical methods and the three cost scenarios previously described.
Henceforth, each figure illustrates the solution obtained by each numerical method adopted for both control and state variables.In the computational measurements, a local optimal solution was attained for all the solvers and cost variations.The number of iterations and function evaluations spent by each direct solver for each weight factor B is presented in Table 1.Table 2 reports the values obtained for the proposed objective functional using the four considered solvers.Starting with the weight factor B = 10 −1 (Fig. 1), all the solvers suggest that for a small investment cost in the control u, the fraction of susceptible population decreases over time t, which means that the control measure adopted is effective on encouraging individuals to spread the marketing message within their social network, inducing, therefore, high levels of infection in a short period of time (see Fig. 1b).All the methods also show that virtually all the target population tends to spread the message, since, for t = 6, the fraction of susceptible population is almost null (see Fig. 1a) and the individuals tend to fully recover from previous infection states (see Fig. 1c), meaning that they had previous contact with the marketing message.
In terms of the application of the control policy (Fig. 1d), mainly LOQO's curve allows inferring that, when the investment cost is low, the control u should be applied gradually from the beginning of the campaign in order to attain the maximum peak of infection.In fact, the higher levels of cost functional (3) are obtained for B = 10 −1 (see Table 2).It should be emphasized that despite of the dynamics of LOQO be different than the obtained with the other solvers for the control signal u, its behavior on the state variables as well as on the cost functional is consistent with the dynamics of all numerical solvers tested.The differences in the dynamics of the solvers can be explained by the numerical steps that promote global convergence.
On the other hand, for a reasonable investment cost in implementing u (Fig. 2), all the solutions obtained show that the levels of information spreading are lower than the ones obtained for B = 10 −1 (see, Fig. 2b).In fact,  whereas for B = 10 −1 the majority of the individuals spreads or has contact with the marketing message, Fig. 2a reveals that, for B = 10, about half of the target population is not stimulated to share the intended message, which in turn induces lower levels of recovery (see Fig. 2c). Figure 2d suggests that facing a reasonable investment cost in control actions, the control measures should be applied with higher intensity at the beginning of the campaign and not over time as in the case of B = 10 −1 .In terms of the proposed objective functional, for B = 10, all the solvers conducted to an approximated solution, which is roughly half than the obtained for B = 10 −1 (see Table 2).This fact is natural, since the adoption of this control weight represents a bigger effort in terms of cost minimization.Regarding the factor B = 10 2 , illustrated in Fig. 3, the simulations related to the fraction of susceptible and recovered population were omitted in order to avoid redundancy.However, for a high investment cost in u, the recorded levels of infection are the lowest (Fig. 3a).In terms of the application of the control policy (Fig. 3b), the results are analogous to the ones obtained for B = 10 (Fig. 2d), but with a lower magnitude.
Transversally, whatever the factor B considered, it is possible to observe that different numerical solvers lead to an analogous dynamics.These dynamics are exactly the same for KNITRO and IPOPT, which corroborates the robustness of the obtained results, notwithstanding KNITRO leads to a smaller number of iterations and function evaluations than IPOPT (see Table 1

Conclusion
Using real numerical data, this paper formulates an optimal control problem to maximize the spreading of information by minimizing costs.Indirect and direct numerical techniques to solve the problem are compared, and the dynamics of the state and control variables under different cost scenarios are studied.
Our results support the conclusion that both indirect and direct methods lead to analogous behaviors in terms of state variables, confirming the reliability of the obtained results.However, in what concerns to the control variable u, the dynamics vary inasmuch as the investment costs in control policies have a direct influence on the objective functional.At this point, numerical simulations showed that when the investment cost related to the adoption of u decreases, the magnitude of u increases.In contrast, when the investment cost associated with the implementation of u increases, the magnitude of u decreases.In addition, for all considered approaches, the cost of implementation of control policies has a significative impact on the dissemination of the information content.Whenever this cost increases, the fraction of population who has contact with the marketing message decreases.As a consequence, the number of infected individuals (who has contact with the marketing message) is smaller.Therefore, the greater the factor B, the smaller the objective functional.
On the other hand, we observed that the numerical estimation of β and γ may give rise to difficulties when real data is considered.It should also be emphasized that despite of FBS be very intuitive, its flexibility might be a problem since it requires explicit derivatives of necessary conditions.On the contrary, direct methods are easier to implement but induce smaller levels of numerical accuracy.
All in all, optimal control strategies combined with indirect and direct numerical methods can be a fruitful tool to study consumer's response to several marketing strategies.As future work, it would be interesting to test other solvers with additional optimal control measures or even perform a cost-effectiveness analysis considering delay differential equations, in order to better describe the dynamics of information epidemics.

Table 1 .
Computational performance indicators of direct solvers.
). Model dynamics for FBS, KNITRO, IPOPT and LOQO under a reasonable investment cost in u.Parameter values: B = 10 and umax = 0.1 .