Reconsider Example 23.3, where an invasive plant species was modeled. The species consists of a finite number of satellites with a main focus. We wish to minimize the total area covered by the invader at the end of a finite time period. A discrete time model is used. The optimal control problem we considered before was

min u

[ r2j,T +B


subject to rj,t+1 = ( rj,t +


ε+ rj,t

) (1− uj,t), rj,0 = ρj ,

0 ≤ uj,t ≤ 1 for j = 0, 1, . . . , T − 1, (24.1)

whereB is a weight parameter, and ρj are the known initial conditions. Recall, the state rj,t is the radius length of focus j at time step t, the control uj,t is the amount of radius decrease for focus j at time step t, and k is the spread rate. To simplify matters, we will fix some of the constants for this lab. The

number of foci n will be set to five, so that there is one main focus and four satellites. The initial values will be set to ρ1 = 0.5, ρ2 = 1, ρ3 = 1.5, ρ4 = 2, ρ5 = 10, so that the main focus is denoted by j = 5. Further, the arbitrary constant ε will be fixed at ε = 0.01. Further, we wish to avoid the transversality condition λj,T = 2rj,T . This

type of boundary condition cannot be handled by the standard forwardbackward sweep (the adapted sweep can be used, of course; see Exercise 24.1). Therefore, we change to

φ(x1,T , . . . , x5,T ) = 5∑

xj,T ,

which gives the transversality conditions λj,T = 1. Therefore, our new optimal control problem is

min u

[ rj,T +B





ε+ rj,t

) (1− uj,t), rj,0 = ρj ,

0 ≤ uj,t ≤ 1. (24.2)

To be clear, we are not claiming (24.1) and (24.2) are equivalent optimal control problems. They will no doubt produce different results. However, by only removing the square, we have not drastically altered the dynamics of the problem. This new problem should have solutions of the same flavor as the original. Before beginning, we discuss the slight changes in the code, as compared

to previous forward-backward sweep routines. The first thing to note, is that while our vectors have been indexed from 0 to T , MATLAB indexes vectors starting with 1. Therefore, the states and adjoints, which should run from 0 to T , actually start at 1 and end with T + 1. Similarly, the controls should end with T −1, but end with T . This is only a superficial difference, occurring entirely in the code. In fact, the t variable, as defined, gives the correct index values. Namely, t(1) = 0, t(2) = 1, . . . , t(T + 1) = T . This t will be used to graph the solutions correctly indexed. The other minor changes involve the differing vector lengths. Before, we

updated the characterization of the control all at once, via vectors. Here, we must be more careful, as the controls uj have one less term than the states and adjoints. Therefore, it is done term-by-term in a for loop. Then, the standard convex combination is used. Second, instead of outputting the solutions as the dummy matrix y, we must use two matrices (y and z). The y variable contains the states (running from 0 to T ), and the controls (running from 0 to T − 1) are given by z. Everything else in the code should be familiar. As we noted in the previous

chapter, instead of a Runge-Kutta 4 routine, the difference equations are solved exactly as stated. In MATLAB run the program lab14. To begin, try the values

B = 1 k = 1 T = 10 , (24.3)

and do not vary any parameters. The optimal solutions are displayed in a different manner than normal. All five state variables are plotted together on the left. For emphasis, the state variable representing the main foci, r5, is plotted in red, with the values denoted by circles. The remaining four states are in blue, marked by x’s. Also, while the x’s and circles denote the values of the states, dotted lines connect the points. This is done in order to make the plots easier to see and read. The five controls are plotted together on the right, with the same color and marker scheme. This first simulation illustrates one of the major differences between con-

tinuous time and discrete time models. Here, choosing the control to be 1 at any time immediately forces the corresponding radius to 0. The optimal strategy depicted is to use no control for the first 8 years, then maximum

radii. This causes the radii to grow naturally for the majority of the time interval, only to be completely eradicated in the last year. We note, with these parameters, any controls with value 1 for one year, and 0 at all others, would produce the same final state values, and are therefore optimal as well. Enter (24.3), varying with k = 0.5. The radii for the first system (with

k = 1) are depicted as before. The main focus for the second system (with k = 0.5) is plotted in purple, with circle markers. The remaining four foci are plotted in green with x markers. The controls are similar. The control strategy is the same for both simulations. Namely, use no control until the last year, then use maximum control. The foci grow naturally until the last year, when they are pushed to 0. The only difference occurs during the natural growth phase, where the result of the different k values is clear. This control strategy, that of maximum control at the end, is optimal in

these cases primarily because of the weight parameter B = 1. Minimizing the final radii values and the total control are of equal importance, which allows this harsh control strategy. Suppose we increase the importance of keeping the control low. Enter (24.3), varying with B = 2. The optimal solutions are identical. We still are neglecting the control too much. Vary with B = 10. This causes a change in strategy. The second system employs controls with less total square sum. As such, the radii do not decrease all the way to 0. Note, all controls in the second simulation steadily increase, and the (second) main focus steadily decreases. However, the other four foci increase at first, and actually end at values higher than their initial positions. We also point out that while the five foci begin at varying values, they end at exactly the same value (0) in the first simulation, and very near the same value in the second. Now examine

B = 5 k = 1 T = 10 . (24.4)

Here, the optimal strategy is a mixed one. For the four minor foci, we should use a steadily increasing control, as we just saw. On the other hand, the main focus should be eradicated via a max control at the end of the time interval. With this medium B value, it is still advantageous to eliminate the main focus, which begins much larger than the others, but to only maintain the others. With

B = 4.19 k = 1 T = 10 (24.5)

the main focus and the largest minor focus are eliminated, while the others are maintained. Try

B = 4.1 k = 1 T = 10 (24.6)


FIGURE 24.1: Optimal radius value r∗5 with B = 10, k = 1, and T = 10 plotted in solid circles, together with r5 with no control applied, in open circles.