## 3 Optimal Control Problem
### Implementation of the Optimal Control Problem  
We wish to find optimal controls numerically by applying a forward-backward iterative method [25]. Initially, a constant control is assumed, and the state equations are solved in the forward-time direction from our standard set of initial conditions. Given this solution to the state equations, the adjoint equations are then solved in the backwards-time direction, beginning with the final time \(T\) and final condition \(\lambda_{3}(T)=c_{5}\) and \(\lambda_{i}(T)=0\) for \(i\in\{1,2,4,5\}\). The controls to be used for the next forward run are then updated using the characterization of the optimal control given in (15), and the forward-backward solution process is repeated with the updated control functions. This process is iterated until the controls and the solutions to all of the differential equations converge to within acceptable numerical tolerances. See [28; 29; 40] for the background on this procedure and the approximation of delay equations.  
To use the forward-backward sweep, we rewrite the adjoint equations for our system as:\[\frac{d\lambda_{1}}{dt}= -\lambda_{1}\left(r\left(1-\frac{2x_{1}+x_{2}+x_{5}}{K}\right)-(1-b_ {1}u_{1})\beta x_{3}\right)-\lambda_{2}\left(-\frac{r}{K}x_{2}+(1-b_{1}u_{1}) \beta x_{3}\right)+\lambda_{5}\frac{r}{K}x_{5}, \tag{16}\] \[\frac{d\lambda_{2}}{dt}= -c_{2}+\lambda_{1}\frac{r}{K}x_{1}-\lambda_{2}\left(r\left(1- \frac{x_{1}+2x_{2}+x_{5}}{K}\right)-(1+a_{1}u_{2})\mu x_{4}-\rho x_{4}\right)\] \[\qquad-\lambda_{3}(1-(1-f)b_{2}u_{1}-fa_{2}u_{2})\pi-\lambda_{5} \left(\rho x_{4}-\frac{r}{K}x_{5}\right)-\lambda_{4}(t+\tau)ax_{4}\chi_{[0,T- \tau]}\] \[\frac{d\lambda_{3}}{dt}= -c_{1}+\lambda_{1}(1-b_{1}u_{1})\beta x_{1}-\lambda_{2}(1-b_{1}u _{1})\beta x_{1}+\lambda_{3}c\] \[\frac{d\lambda_{4}}{dt}= -\lambda_{2}\left(-(1+a_{1}u_{2})\mu x_{2}-\rho x_{2}\right)+ \lambda_{4}d(1-a_{3}u_{2})-\lambda_{5}(\rho x_{2}-vx_{5})-\lambda_{4}(t+\tau )ax_{2}\chi_{[0,T-\tau]}\] \[\frac{d\lambda_{5}}{dt}= -\lambda_{1}\left(-\frac{r}{K}x_{1}+q\right)+\lambda_{2}\frac{r} {K}x_{2}-\lambda_{5}\left(r\left(1-\frac{x_{1}+x_{2}+2x_{5}}{K}\right)-q-\
u x _{4}\right)\]  
Notice that in the time interval \([T-\tau,T]\), the advance terms (those with argument \(t+\tau\), found in the second and fourth equations) drop out, so we have five ordinary differential equations. On the interval \([0,T-\tau]\), we once again have advance equations, but the solutions to the ODEs on \([T-\tau,T]\) provide the required initial data to solve these equations. Thus, the adjoint equations are advance differential equations on \([0,T-\tau]\) and ordinary differential equations on \([T-\tau,T]\), subject to the final condition \(\lambda_{3}(T)=c_{5}\) and \(\lambda_{i}(T)=0\) for \(i\in\{1,2,4,5\}\).  
For our numerical simulations, we used the built-in MATLAB numerical delay differential equation solver, dde23. This tool is not capable of solving advance equations directly, so we made a change of variables to convert the advance differential equations to a system of delay differential equations. Specifically, we define a new reversed-time variable \(\sigma=T-t\) and new adjoint variables \(L_{i}(\sigma)=\lambda_{i}(T-\sigma)=\lambda_{i}(t)\). In terms of these new variables, adjoint equations for \(L_{i}(\sigma)\) are ordinary differential equations on \([0,\tau]\) and delay differential equations on \([\tau,T]\), subject to the initial condition \(L_{3}(0)=c_{5}\) and \(L_{i}(0)=0\) for \(i\in\{1,2,4,5\}\). The new system was solved numerically using dde23 as part of the implementation of the forward-backward method to find the optimal control.