A mathematical model for antibiotic control of bacteria in peritoneal dialysis associated peritonitis

1.   Introduction

It is well-known that the optimal feedback control can be synthesized by the dynamic programming (DP) approach, while finding the optimal control in feedback form is generally argued as the Holy Grail of control theory ([16]). However, before we arrive at the feedback control by this approach, there are two difficulties that need to face. The notion of viscosity solution to partial differential equation (PDE) overcomes the first one that the HJB equation usually has no classical solution regardless of the smoothness of its coefficients ([7,8,9]). Under this framework, such kind of weak solution simultaneously holds the existence and uniqueness. The value function of the optimal control problem is the unique viscosity solution to the HJB equation. Nevertheless, for most of optimal control problems, it is usually impossible to find the analytical viscosity solution, which is the second difficulty to hinder the utilization of the DP approach. The numerical solution is almost the only practically significant way for finding the optimal control via solving the HJB equation. Luckily enough, nowadays there are some literatures available devoted to these topics where the analysis of DP methods for deterministic and stochastic control problems is presented in detail ([3,6,11,12]).

Besides the Bellman DP principle, another milestone work in control theory is contributed to the Pontryagin maximum principle, which establishes a necessary optimality condition satisfied by the investigational optimal control systems consisted of the state equation and its adjoint system ([23,26]). Nevertheless, in essence, the control presented by the Pontryagin maximum principle is of open-loop, which is usually less than satisfactory. Moreover, it is notorious for its sophistication and the fact that it rarely gives an insight into the structure of the optimal controls ([24]). Generally speaking, a feedback control system is considered superior to an open-loop system. Many unnecessary disturbances and noise signals from outside the system can be rejected by using feedback ([4]). Furthermore, even from the numerical solution point of view, for the two-point boundary problem obtained by the necessary condition of optimality, the recognized shooting method has to overcome the difficulty of "guess" for the initial data to start the iterative numerical process ([14,27]), which is generally not an easy job.

People need to make great efforts to develop more efficient algorithms to find the optimal control. Two real reasons for this could not be clearer. One is due to the lack of the analytical solution to the practical optimal control problems. We must numerically solve the investigational problem. Furthermore, the ultimate aim of the control theory is to use it in the real world. It also coincides with the requirement of the information age. The other reason is caused by the curse-of-dimensionality, the inherent defect of the DP approach ([22]). Thanks to many model reduction techniques ([1,3,13,20,19]), this problem can be partially circumvented. It makes some effective algorithms for low dimensional problems, thus initially confined to only "toy'' problems, can be generalized to deal with much more complicated practical optimal control problems. The upwind finite-difference scheme is such a well-adapted algorithm and has been successfully applied to many examples ([13,14,15]). And most of all, its convergence has been rigorously proven in [28].

Certainly it is not easy to find an optimal feedback control by the HJB-based approach and lots of complicated computations are usually involved. Even so, on this topic, there are still many interesting results available in literature, such as [1,5,17,19,22,29], to name just a few. We can unfortunately select only some more relevant research to mention here. For two practical optimal control problems, [14,15] obtain the optimal feedback controls by repeatedly calling the algorithm to the HJB equation. In [2,10,18], the HJB equation is firstly solved to gain the value function. And then the convex combination of numerical value function is applied to obtain the approximation of optimal feedback control. Note that the obtained numerical solution of control function is completely ignored and not directly applied to obtain optimal feedback control. A class of explicit Markov chain approximation methods are introduced and studied in detail in [21]. And these methods are tailored to solving the HJBs for both stochastic and deterministic continuous time nonlinear optimal control problems on finite and infinite horizons. The paper [20] considers a distributed volume control problem in infinite time region and obtains the optimal control by the convex combination of the control values at each point of the polyhedron, which can be interpreted as a special interpolation. The idea of using the interpolation to find the optimal control is also mentioned in [13] but it is kept mostly to one-dimensional linear interpolation for one-dimensional problem. The multi-dimensional problems are not involved in the analysis.

In this paper, we are concerned with a feedback design for numerical solution to the optimal control problem in finite time horizon under the framework of the DPVS approach. To synthesize the optimal feedback control, we solve the HJB equation by the upwind finite-difference scheme. Different from the usual existing algorithms, the numerical control function is interpolated in turn to gain the approximation of optimal feedback control-trajectory pair, which avoid solving the HJB equation repeatedly, thus efficaciously promote the computation efficiency. Moreover, the design can deal with the optimal control problem equipped with multiple controls. And more generally, here, the adopted interpolation is not limited to one-dimensional linear interpolation. Other interpolation techniques including multi-dimensional linear interpolation, cubic spline interpolation, nearest-neighbor interpolation and cubic-Hermite interpolation go well likewise in this new design.

The remainder of this paper is organized as follows. In Section 2, we revisit the procedure of the DPVS approach and present the upwind finite-difference scheme to the HJB equation. An interpolation algorithm then follows, which consists of the core of this paper. The effectiveness of the algorithm is proven in Section 3 by solving five quite distinct examples from each other before some concluding comments are made in Section 4.

2.   DPVS approach and interpolation algorithm

Set $ T>0 $ and $ U\subset \mathbb{R}^m $ be a closed, bounded and convex set. Denote by $ \mathcal{U}[0,T] = L^\infty([0,T]; U) $ the admissible control set. We consider the following control system

where $ y:[0,T]\rightarrow \mathbb{R}^n $ is the state, the control $ u(\cdot)\in \mathcal{U}[0,T] $ and $ y_0\in \mathbb{R}^n $ is a given initial value. Assume that the mapping $ f:[0,T]\times \mathbb{R}^n\times U \rightarrow \mathbb{R}^n $ is continuous. Take the cost functional as

where both the running cost $ L:[0,T]\times\mathbb{R}^n\times\mathbb{R}^m\rightarrow \mathbb{R} $ and the terminal cost $ h:\mathbb{R}^n\rightarrow \mathbb{R} $ are continuous. The optimal control problem is to minimize (2.2) over $ \mathcal{U}[0, T] $, i.e.,

subject to the control system (2.1).

In order to present the DP principle, we consider a family of dynamic optimal control problems. Namely, over $ \mathcal{U}[t, T] = L^\infty([t,T]; U) $, we minimize the cost functional

subject to

where $ (t,x)\in [0,T)\times \mathbb{R}^n $ is the initial time and state. Define the value function by

which is the associated minimal cost functional. It satisfies the DP principle ([3])

for any $ \tau\in [t,T) $.

Then, if the value function $ v $ is smooth enough, such as $ v \in C^1([0,T]\times \mathbb{R}^n) $, the following HJB equation can be derived from the DP principle

where $ \nabla_x $ denotes the gradient with respect to $ x $. By solving the HJB equation (2.5), we next turn to the feedback synthesis, i.e., the construction of the optimal control $ u^*(t) $ with the optimal state $ y^*(t) $. For any mapping $ \bar{u}(t,x) $ satisfying

we set

for almost all $ t\in[0,T] $, where $ y^*(t) $ satisfies

with $ y^*(0) = y_0 $. Then we have

for almost all $ t \in[0,T] $. Thus $ u^*(\cdot) $ is an optimal control of (2.3) from [3]. Here, $ \bar{u} $ in (2.6) is called as the control function or the optimal feedback mapping. When the mapping $ \bar{u} $ is known in (2.7), the system (2.7) is a closed loop system.

The DPVS approach tells us that once a viscosity solution of the HJB equation (2.5) is obtained numerically, we are able to find the numerical solution of the feedback law by the feedback synthesis above.

Next we revisit the upwind finite-difference scheme to the HJB equation (2.5). Actually, it is such a well-adapted algorithm and has been successfully applied to many examples ([14,13,15]). And most of all, its convergence has been rigorously proven in [28].

Take a polygonal region $ [0,T]\times \prod_{k = 1}^n[a_k,b_k] $ as the computational domain. Set the time step $ \Delta t = - \frac{T}{N} $ and denote $ t_j = T+j\Delta t,j = 0,1,2,\ldots,N $, which divide the time domain $ [0, T] $ into $ N $ equal parts. Similarly, for each $ k = 1,2,\ldots,n $, divide $ [a_k,b_k] $ into $ M_k $ equal parts. Denote $ \Delta x_k = \frac{b_k - a_k}{M_k} $, $ k = 1,2,\cdots,n, $ the space steps. And let $ i = \left(i_1,i_2,\cdots,i_n \right) $, $ ik^+ = \left(i_1,\ldots,i_{k-1},i_k+1,i_{k+1},\ldots,i_n \right) $, $ ik^- = \left(i_1,\cdots,i_{k-1},i_k-1,i_{k+1},\cdots,i_n \right) $, where $ i_k = 0,1,\cdots, M_k $, $ k = 1,2,\ldots, n $. Denote the point $ x_i = \left(a_1+i_1\Delta x_1,a_2+i_2 \cdot \Delta x_2,\ldots,a_n+i_n\Delta x_n \right) $. Then the upwind finite-difference scheme of the HJB equation (2.5) can be written as

for $ i_k = j+1,j+2,\ldots,M_k-j-1 $, $ k = 1,2,\ldots,n $ and $ j = 0,1,\ldots,N-1 $. Here, the sign$ \; z $ denotes the sign of $ z $. And $ u_i^j, v_i^j $ are, respectively, the approximations of $ \bar{u}(t_j,x_i) $ and $ v(t_j,x_i) $. Then, denote

In the following, we construct a HJB-based feedback algorithm for (2.3). Its detailed steps are listed as follows.

Step 1. Initial partition on time and space. Divide the calculation domain into mesh as mentioned in Section 2 above. Denote $ \left(u^j, y^j \right) $ by the optimal control-trajectory pair at $ t^j = -j\Delta t $ for $ j = 0,1,\ldots, N $.

Step 2. Compute the approximations of value function and control function. Based on the numerical scheme (2.8), set $ \gamma_k = \frac{\Delta t}{\Delta x_k} $ for $ k = 1,2,\ldots,n $.

Sub-step 2.1. By the same space steps $ \Delta x_k $, $ k = 1,2,\ldots,n $, extend the initial solution region $ \prod_{k = 1}^n [a_k,b_k] $ to a new one $ \prod_{k = 1}^n \left[a_k-(N+1)\Delta x_k,b_k+(N+1)\Delta x_k \right] $ by adding the $ 2(N+1) $ nodes. (This is a rough expansion and more precise extension steps can be found in [29]). Renumber the points on the new region $ \prod_{k = 1}^n \left[a_k-(N+1)\Delta x_k,b_k+(N+1)\Delta x_k \right] $ like in Section 2. Denote $ ii = (ii_1,ii_2, $ $ \ldots,ii_n) $, $ iik^+ = (ii_1,\ldots,ii_{k-1},ii_k+1,ii_{k+1},\ldots,ii_n) $, $ iik^- = (i_1,\ldots,ii_{k-1},ii_k-1,ii_{k+1},\ldots,ii_n) $, where $ ii_k = 0,1,\ldots, M_k+2(N+1) $ and $ k = 1,2,\ldots, n $. Define $ x_{ii} = (a_1+(ii_1-N-1)\Delta x_1, \; a_2+(ii_2-N-1)\Delta x_2,\ldots, \; a_n+(ii_n-N-1)\Delta x_n) $.

Sub-step 2.2. Initialize $ { } \left\{v_{ii}^0 \right\}_{ii_k = 0}^{M_k+2(N+1)} $ and $ { } \left\{u_{ii}^0 \right\}_{ii_k = 0}^{M_k+2(N+1)} $ as follows.

Sub-step 2.3. For $ j = 0 $ to $ N-1 $, compute $ v_{ii}^{j+1} $ and $ u_{ii}^{j+1} $ as follows. For $ ii_k = j+1 $ to $ M_k+2N-j+1 $, $ k = 1,2,\ldots ,n $, set

Sub-step 2.4. Choose the numerical solution of value function and control function over $ [0, T]\times \prod_{k = 1}^n [a_k,b_k] $. Denote them by $ \left\{\{v_i^j\}_{i_k = 0}^{M_k} \right\}_{j = 0}^N $ and $ \left\{\{u_i^j\}_{i_k = 0}^{M_k} \right\}_{j = 0}^N $, $ k = 1,2,\ldots,n $, respectively. Here, $ v_i^j $ is the approximation of value function $ v(t_j,x_i) $ and $ u_i^j $ is the approximation of control function $ \bar{u}(t_j,x_i) $.

Step 3. Compute the optimal pairs $ \mathit{\boldsymbol{\{(u^j,y^j)\}_{j = 0}^N}} $.} For initial setting, let $ y^0 = y_0 $. Set $ j = 0 $.

Sub-step 3.1. Calculate the $ j $-th optimal feedback control $ u^j $ by applying the numerical solution $ \left\{u_i^{N-j} \right\}_{i_k = 0}^{M_k} $, $ k = 1,2,\ldots,n $, obtained in Step 2. Choose the linear interpolation (or cubic spline interpolation, nearest-neighbor interpolation and cubic-Hermite interpolation) for point $ y^j $ by using the value $ \left\{u_i^{N-j} \right\}_{i_k = 0}^{M_k} $ at known points $ \{x_i\}_{i_k = 0}^{M_k} $, $ k = 1,2,\ldots,n $. Actually, the value obtained by interpolation is an approximation of $ \bar{u}(t^j,y^j) $, which is denoted as $ \tilde{u}(t^j,y^j) $. Then, we set

Sub-step 3.2. Calculate the $ (j+1) $-th optimal trajectory $ y^{j+1} $. Discrete the state equation by the forward difference:

With the known $ u^j $ and $ y^j $, we solve the finite difference equation above to obtain the $ (j+1) $-th optimal trajectory $ y^{j+1} $.

Sub-step 3.3. Iterate for the next time instant. Set $ j = j+1 $. If $ j>N $, end the procedure. Otherwise, go to Sub-step 3.1. Continue the iteration until we get all optimal feedback control-trajectory pairs $ \left\{(u^j,y^j) \right\}_{j = 0}^N $.

Remark 1. In Sub-step 3.1, we can also adopt other interpolation methods, such as polynomial interpolation, average interpolation, the Newton interpolation, Lagrange's interpolation and so on ([25]). Similarly, other difference method can be used to discrete state equation in Sub-step 3.2.

3.   Simulations

This section is devoted to five numerical examples to verify the effectiveness of the algorithm. The simulations are implemented by using MATLAB programme in a desktop computer with a core process (i5-6600M) and 8GB of random access memory.

Example 1. Consider the following control problem ([31])

where $ \mathcal{U} = L^\infty([0,1];U) $, $ U = [0,1] $, $ y_0\in \mathbb{R} $ is a given initial value. Its corresponding HJB equation is

Its exact value function $ v $ is

and the control function is

Here, we choose the region $ D = [0,1] \times \left[-\frac{3}{2},\frac{3}{2}\right] $, i.e., $ 0\leq t\leq 1 $ and $ -\frac{3}{2}\leq x\leq\frac{3}{2} $. The number of nodes $ N $ on the $ t $ axis is chosen as $ 41 $. And the number of nodes $ M $ on the $ x $ axis is set as $ 31 $. Then the time step $ \Delta t = - \frac{1}{40} $ and the space step $ \Delta x = \frac{1}{10} $. Adopting the linear interpolation, we numerically solve this example by the interpolation algorithm presented in Section 2 above. And the results are plotted in Figures 1 and 2.

Moreover, in order to show the performance of the presented interpolation algorithm, two existing algorithms are additionally chosen to solve Example 1, simultaneously. The corresponding results are also given in Figures 1 and 2. And the comparisons on numerical simulations are listed in Table 1.

Actually, one chosen algorithm is the minimisation of the Hamiltonian. People can refer to Algorithm 3.2 in [13] for details. For convenience, we use Algorithm 1 to denote it. In the other one algorithm, the optimal feedback control is also gained by minimizing the Hamiltonian function. But it is different from Algorithm 1, in which the value function is obtained by the convex combination of numerical value function ([18]). We name it as Algorithm 2.

When $ y_0 = -\frac{1}{2} $, it is easy to obtain the optimal control $ u^*(t)\equiv 0 $, $ t\in[0,1] $, and the optimal trajectory $ y^*(t)\equiv -\frac{1}{2} $, $ t\in [0,1] $. Figure 1(a) shows the numerical solutions of optimal feedback control agrees well with the exact solution for these three algorithms. In Figure 1(b), the numerical solution of optimal trajectory also agrees well with its exact solution. When $ y_0 = \frac{1}{2} $, the optimal control is $ u^*(t)\equiv 1 $, $ t\in[0,1] $, and the optimal trajectory is $ y^*(t) = \frac{1}{2}e^t $, $ t\in [0,1] $. Figure 2 shows both the numerical solutions of optimal feedback control and trajectory agree well with the exact solutions when initial state $ y_0 = \frac{1}{2} $.

For these three algorithms, Table 1 lists the computed errors in the maximum norm between the numerical solutions and exact solutions including the optimal feedback control $ u^* $ and optimal trajectory $ y^* $ with initial state $ y_0 = - \frac{1}{2} $ and $ y_0 = \frac{1}{2} $. So small error estimates tell us that the numerical solutions obtained by the interpolation algorithm, Algorithms 1 and 2 are very close to the exact solutions.

Moreover, the CPU times are listed in Table 1. They are the computing time of these three algorithms after we obtain the numerical control function and numerical value function by upwind finite-difference method. By comparing the CPU times, we find that the interpolation algorithm takes the least time than Algorithms 1 and 2 with the similar error estimation. It illustrates that linear interpolation is meaningful for finding an approximation of the optimal feedback control, although the control function (3.1) is discontinuous. In all, numerical results shows that the interpolation algorithm is not only feasible but also efficient for gaining the optimal feedback control.

Example 2. Consider the following control problem ([30])

where $ \mathcal{U} = L^\infty([0,1];U) $, $ U = \mathbb{R} $. Its corresponding HJB equation can be read as

The exact solution to the HJB equation (3.2) is

and the control function is

With (2.6), (2.7) and (3.3), we obtain the optimal feedback control $ u^*(t) = \frac{-2e^{-t}}{1+e^{-2}} $ and the optimal trajectory $ y^*(t) = \frac{e^{-t}+e^{t-2}}{1+e^{-2}} $.

Now we numerically solve Example 2. Choose the computational region $ [0,1]\times[0,1] $, i.e., $ 0\leq t\leq 1 $ and $ 0\leq x\leq 1 $. Here we choose a uniform mesh along $ t $ axis with $ N = 41 $ nodes, i.e., $ \Delta t = - \frac{1}{40} $, and a uniform mesh along $ x $ axis with $ M = 21 $ nodes, i.e., $ \Delta x = \frac{1}{20} $. In this example, the chosen interpolation is also the linear interpolation. The obtained numerical solutions are plotted in Figure 3. By this figure, we see that both the numerical solutions of optimal feedback control and trajectory agree well with the exact solutions. Moreover, in Table 2, we list the maximum errors between the numerical solutions and exact solutions of optimal feedback control and trajectory. They are, respectively, $ {\rm{0.025896}} $ and $ {\rm{0.017297}} $, which are very small. It also shows that the algorithm is effective.

In addition, in Examples 1 and 2, the control problems are all one-dimensional problems. Next, we present other two optimal control examples with two space dimensions.

Example 3. Consider the following two-dimensional control problem ([29]) with two state variables and one control.

where $ \mathcal{U} = L^\infty([0,1];U) $, $ U = [0,1] $, $ (x_0,y_0) \in \mathbb{R}^2 $. Its corresponding HJB equation is

In this example, the exact solution of the value function is

and the control function is

From (2.6), (2.7) and (3.5), when initial states $ (x_0,y_0) = \left(-\frac{1}{2}, - \frac{1}{2} \right) $, the optimal feedback control is $ u^*(t)\equiv 0, \; t \in [0,1] $ and the optimal trajectory is $ (x^*(t),y^*(t)) \equiv \left( - \frac{1}{2}, -\frac{1}{2} \right), \; t \in [0,1] $. When the initial states $ (x_0,y_0) = \left(\frac{1}{2}, \frac{1}{2}\right) $, the optimal feedback control is $ u^*(t)\equiv 1 $ and the optimal trajectory is $ (x^*(t),y^*(t)) = \left(\frac{1}{2}e^t,\frac{1}{2}e^t\right), \; t\in[0,1] $.

Next we use the interpolation algorithm to solve the HJB equation (3.4) in the range $ [0,1] \times \left[-\frac{3}{2}, \frac{3}{2}\right]^2 $, i.e., $ 0\leq t \leq 1 $, $ - \frac{3}{2}\leq x\leq \frac{3}{2} $, $ - \frac{3}{2}\leq y\leq \frac{3}{2} $. Here we still choose a uniform mesh on $ t $ axis, $ x $ axis and $ y $ axis. The nodes $ N $ on $ t $ axis is set as $ 41 $, i.e., $ \Delta t = - \frac{1}{40} $. The nodes on $ x $ and $ y $ axis are set as $ M_1 = M_2 = 31 $, i.e., $ \Delta x = \Delta y = \frac{1}{10} $. Again, the linear interpolation is chosen to be applied in the algorithm.

From Figures 4 and 5, we see that, even with different initial values, the images of numerical solutions agree well with the exact solutions for the optimal feedback control and trajectory.

Example 3 shows that, even for two-dimensional control problem, the linear interpolation is meaningful for finding an approximation of the optimal feedback control. On the other hand, we find that the discontinuity of the control function (3.5) does not hinder us to find the optimal feedback control.

Furthermore, all three examples above are equipped with one control. And only the linear interpolation is adopted in the algorithm. And to show the effectiveness and efficiency of the interpolation algorithm in an even better fashion, we consider an optimal control problem with two states and two controls in Example 4. And four interpolations including linear interpolation, cubic spline interpolation, nearest-neighbor interpolation and cubic-Hermite interpolation are successively adopted in the investigation.

Example 4. Consider the following optimal control problem with two state variables and two controls.

where $ (x_{10},x_{20}) \in \mathbb{R}^2 $, $ \mathcal{U} = L^\infty([0,1];U) $, $ U = \{(p,q)\in \mathbb{R}^2 \;|\; |p|\leq 1, |q| \leq 1, |p+q^2|\leq 1\} $. Denote by $ |\cdot| $ the absolute value. The HJB equation of (3.6) can be written as

Here we choose the computational region $ [0,1]\times[-1,1]^2 $, i.e. $ 0\leq t\leq 1 $, $ -1\leq x_1\leq 1 $ and $ -1\leq x_2\leq 1 $. The number of nodes $ N $ on $ t $ axis is chosen as $ 41 $. And the number of nodes on $ x_1 $ and $ x_2 $ axis are set as $ M_1 = M_2 = 41 $. Then the time step $ \Delta t = - \frac{1}{40} $ and the space steps $ \Delta x_1 = \Delta x_2 = \frac{1}{20} $. To visualize the numerical result, we plot the value function $ v $, two control functions $ \bar{u}_1 $ and $ \bar{u}_2 $, as well as $ \bar{u}_1+\bar{u}_2^2 $ at different across-sections in Figures 6-8. From Figures 6(b-d), 7(b-d) and 8(b-d), we see that both the numerical controls $ \bar{u}_1 $ and $ \bar{u}_2 $, and $ \bar{u}_1+\bar{u}_2^2 $ satisfy the given constraints.

Moreover, in this example, the exact value function $ v $ and control function $ (\bar{u}_1,\bar{u}_2) $ of the HJB equation (3.7) are unknown in literature. Although the control variables have so many constraints including nonlinear constraint, when $ (x_{10}, x_{20}) = \left(\frac{1}{2}, - \frac{1}{2}\right) $, by calculation we find the optimal controls $ u_1^*(t) = \frac{-e^{-t}}{1+e^{-2}} $, $ u_2^*(t) = \frac{e^{-t}-e^{t-2}}{2(1+e^{-2})} $, $ t\in [0,1] $, and the optimal trajectories $ x_1^*(t) = \frac{e^{-t}+e^{t-2}}{2(1+e^{-2})} $, $ x_2^*(t) = -\frac{e^{-t}+e^{t-2}}{2(1+e^{-2})} $, $ t\in [0,1] $. Furthermore, in this example, we utilize the linear interpolation and the cubic spline interpolation to find the optimal feedback control. And their respective numerical results are plotted in Figures 9 and 10. From these two figures, we see that, no matter if it is the linear interpolation or the cubic spline interpolation, all numerical solutions have a good agreement with the exact solutions.

Then, in order to investigate the influence of different interpolation on the interpolation algorithm, we provide the error results of four interpolation cases for this example. In Table 4, we list the computational errors in the maximum norms between the numerical solutions and exact solutions including optimal control $ (u_1^*,u_2^*) $ and optimal trajectory $ (x_1^*,x_2^*) $ with initial state $ x_{10} = \frac{1}{2} $ and $ x_{20} = - \frac{1}{2} $.

In the linear interpolation case, they are $ {\rm{0.038706}} $ and $ {\rm{0.02681}} $, respectively. As for the cubic spline interpolation case, they are $ {\rm{0.038706}} $ and $ {\rm{0.026801}} $, respectively. In nearest-neighbor and cubic-Hermite interpolation cases, we also get very small errors as well. Although the errors in nearest-neighbor case are little bit bigger than others, all of them show a nice approximation to the exact solutions. Relatively speaking, by Table 4, we find that the nearest-neighbor interpolation is the worst among these four interpolations. And other three interpolations almost perform equally well. From these small error estimates, we know that the numerical solutions obtained by the interpolation algorithm are very close the exact solutions. The experimental results of Example 4 demonstrate again that interpolation algorithm is extremely effective.

Let me emphasize that, for four examples above, all optimal controls are continuous, while the case of optimal control with jumps is not considered. Next, we use the interpolation algorithm to solve a control problem in which optimal control has jumps. We attempt to show that the interpolation algorithm is also feasible for this kind of problems.

Example 5. Consider the following optimal control problem with two state variables and one control.

where $ \mathcal{U} = L^\infty([0,2\pi];U) $, $ U = [-1,1] $. The HJB equation of (3.8) can be written as

Here we choose the computational region $ (t, x_1, x_2) \in [0,2\pi]\times \left[-\frac{3}{2},\frac{3}{2}\right]^2 $. The number of nodes $ N $ on $ t $ axis is chosen as $ 41 $. And the number of nodes on $ x_1 $ and $ x_2 $ axis are set as $ M_1 = M_2 = 11 $. Then the time step $ \Delta t = - \frac{\pi}{20} $ and the space steps $ \Delta x_1 = \Delta x_2 = \frac{3}{10} $. In this case, for the given initial state $ (x_{10}, x_{20}) = (0, 0), $ the unique optimal control $ u^*(t) $ and the corresponding optimal trajectory $ (x^*_1(t), x^*_2(t)) $ of the system (3.8) can be obtained analytically by

and

From (3.9), we see that the optimal control is not continuous. It has two jumps in $ [0, 2\pi] $. For this example, we execute the computation by the interpolation algorithm above. We choose the linear interpolation and the Heun method in Sub-step 3.2. The simulation results are plotted in Figure 11, in which we find all numerical solutions have a good agreement with the exact solutions to optimal control and trajectory. The computational errors in the maximum norms between the numerical solutions and exact solutions to optimal control $ u^* $ and optimal trajectory $ (x_1^*,x_2^*) $ are, respectively, 4.8000e-05 and 0.0911. These small error estimates demonstrate that interpolation algorithm is likewise effective for finding the numerical optimal controls for solving control problem with optimal control with jumps.

4.   Conclusion

In this paper, we construct an interpolation algorithm for finding an approximation of optimal feedback control. And five experimental examples with different constraints (even nonlinear constraints) on control are successfully solved to demonstrate the effectiveness of the presented interpolation algorithm. The numerical results show that the algorithm with interpolation is very efficient for solving practical optimal control problems.

Acknowledgments

The authors would like to thank the editor and the anonymous referees for their very careful reading and constructive suggestions that improve the manuscript substantially.