Research Article  Open Access
Kin Wei Ng, Ahmad Rohanin, "Solving Optimal Control Problem of Monodomain Model Using Hybrid Conjugate Gradient Methods", Mathematical Problems in Engineering, vol. 2012, Article ID 734070, 14 pages, 2012. https://doi.org/10.1155/2012/734070
Solving Optimal Control Problem of Monodomain Model Using Hybrid Conjugate Gradient Methods
Abstract
We present the numerical solutions for the PDEconstrained optimization problem arising in cardiac electrophysiology, that is, the optimal control problem of monodomain model. The optimal control problem of monodomain model is a nonlinear optimization problem that is constrained by the monodomain model. The monodomain model consists of a parabolic partial differential equation coupled to a system of nonlinear ordinary differential equations, which has been widely used for simulating cardiac electrical activity. Our control objective is to dampen the excitation wavefront using optimal applied extracellular current. Two hybrid conjugate gradient methods are employed for computing the optimal applied extracellular current, namely, the HestenesStiefelDaiYuan (HSDY) method and the LiuStoreyConjugateDescent (LSCD) method. Our experiment results show that the excitation wavefronts are successfully dampened out when these methods are used. Our experiment results also show that the hybrid conjugate gradient methods are superior to the classical conjugate gradient methods when Armijo line search is used.
1. Introduction
Many science and engineering problems can be expressed in the form of optimization problems that are governed by partial differential equations (PDEs). This class of optimization problems is known as PDEconstrained optimization problem. PDEconstrained optimization problems arise widely in many areas such as environmental engineering [1], atmospheric science [2], biomedical engineering [3], and aerodynamics [4, 5]. However, solving such PDEconstrained optimization problems is a challenging task due to their size, complexity, and infinite dimensional nature. In order to deal with these numerical challenges, different approaches such as preconditioning [6–8] and parallel computing [9] have been proposed by researchers.
The specific PDEconstrained optimization problem considered in the present paper is the optimal control problem arising in cardiac electrophysiology, where the monodomain model appears as the governing equations. The monodomain model consists of a parabolic PDE coupled to a system of nonlinear ordinary differential equations (ODEs), which has been widely used for simulating cardiac electrical activity [10–13]. Thus, the above optimal control problem can be generally known as optimal control problem of monodomain model.
The optimal control problem of monodomain model was first proposed by Nagaiah et al. [14] with the control objective to dampen the excitation wavefront of the transmembrane potential using optimal applied extracellular current. Three classical nonlinear conjugate gradient methods have been applied by Nagaiah et al. [14] to solve the optimal control problem, namely, PolakRibièrePolyak (PRP) method [15, 16], the HagerZhang (HZ) method [17], and the DaiYuan (DY) method [18]. Later, Ng and Rohanin [11] employed the modified conjugate gradient method for solving the optimal control problem of monodomain model. For the present paper, we present the numerical solution for the optimal control problem of monodomain model using hybrid conjugate gradient methods.
The structure of the paper is organized as follows. Section 2 presents the optimal control problem of monodomain model with Rogersmodified FitzHughNagumo ion kinetic. In Section 3, we discuss the numerical approach used to discretize the optimal control problem. Next, the optimization procedure is presented in Section 4 while the numerical experiment results are given in Section 5. Finally, we conclude our paper with a short discussion of our work in Section 6.
2. The Optimal Control Problem of Monodomain Model
Let be the computational domain with Lipschitz boundary and be the final simulation time. We further set and . The optimal control problem of monodomain model with Rogersmodified FitzHughNagumo ion kinetic is therefore given by the following: where Here is the control domain, is the observation domain, is the regularization parameter, is the outer normal to , is the transmembrane potential, is the extracellular current density stimulus, is the constant scalar used to relate the intracellular and extracellular conductivity tensors, is the intracellular conductivity tensor, is the surfacetovolume ratio of the cell membrane, is the membrane capacitance, is the current density flowing through the ionic channels, is the ionic current variable, is the prescribed vectorvalue function, is the threshold potential, is the plateau potential, and , , , are positive parameters.
Notice that the optimal control problem of monodomain model consists of (2.1)–(2.4). Equation (2.1) is the cost functional that we need to minimize, with and as the state variables while as the control variable. The control variable is chosen such that it is nontrivial only on the control domain and extended by zero on . The monodomain model, as given by (2.2) appears as the only constraint for the optimal control problem of monodomain model. Lastly, (2.3) and (2.4) are obtained from the Rogersmodified FitzHughNagumo model [19] for representing ion kinetics.
According to Kunisch and Wagner [20], the controltostate mapping is well defined for the optimal control problem of monodomain model, that is, . Consequently, the cost functional in (2.1) can be rewritten as the follwong: where (2.5) is known as the reduced cost functional.
3. Discretization of the Optimal Control Problem
We adopt the optimizethendiscretize approach to discretize the optimal control problem of monodomain model. This classical approach first derives the infinite dimensional optimality system, and the resulting optimality system is then discretized.
3.1. FirstOrder Optimality System
For deriving the infinite dimensional optimality system, the Lagrange functional is formed as follows: where and are the adjoint variables used to adjoin the constraints in (2.2) to the reduced cost functional in (2.5). The firstorder optimality system (infinite dimensional) is obtained by equating the partial derivatives of (3.1) with respect to the state , adjoint , and control variables equal to zero: where denotes the transmembrane potential in the observation domain and denotes the partial derivative with respect to . We further obtain the boundary and terminal conditions:
Next, the state and adjoint systems can be formed using the boundary and initial conditions in (2.2) as well as (3.2)–(3.8). As a result, the state system is given by the following: and the adjoint system is given by the following: Also, by utilizing (3.6), the reduced gradient is given by the following: In order to compute the reduced gradient in (3.11), we are required to solve the state system in (3.9) to obtain the value of , and then the adjoint system in (3.10) to obtain the value of .
3.2. Numerical Discretization
Once the firstorder optimality system (as defined in (3.9)–(3.11)) has been derived, the numerical discretization needs to be carried out. However, the state system in (3.9) and the adjoint system in (3.10) are computationally demanding since they consist of a parabolic PDE coupled to a system of nonlinear ODEs. In order to reduce this computational demand, the operator splitting technique as proposed by Qu and Garfinkel [21] is applied to (3.9) and (3.10) for decomposing the systems into subsystems that are much easier to solve. As a result, the state system in (3.9) becomes while the adjoint system in (3.10) becomes
For the discretization procedure, the linear PDEs in (3.12) and (3.13) are discretized with the CrankNicolson method in time and the Galerkin finite element method in space. On the other hand, the nonlinear ODEs in (3.12) and (3.13) are discretized with forward Euler method in time. The discretized state system is therefore given by the following: Here and are the local timesteps, is the mass matrix, and is the stiffness matrix. On the other hand, the discretized adjoint system is given by the following:
4. Optimization Procedure
Two hybrid conjugate gradient methods are used to solve the discretized optimal control problem of monodomain model, namely, the HestenesStiefelDaiYuan (HSDY) method [22] and the LiuStoreyConjugateDescent (LSCD) method [23]. These hybrid conjugate gradient methods combine the good numerical performance of the HestenesStiefel (HS) and LiuStorey (LS) methods with the strong convergence properties of the DaiYuan (DY) and ConjugateDescent (CD) methods.
4.1. Hybrid Conjugate Gradient Methods
Starting from an initial guess , our control is updated using the following recurrence: where is the steplength computed by the Armijo line search and is the search direction. Given an initial steplength and , the Armijo line search chooses such that On the other hand, the search direction is defined by where is the conjugate gradient update parameter. The conjugate gradient update parameters for the HSDY and LSCD methods are given as follows: where
4.2. Optimization Algorithm
In this section, we present the optimization algorithm for solving the discretized optimal control problem of monodomain model using the HSDY and LSCD methods. The optimization algorithm for these two hybrid conjugate gradient methods is therefore given as follows.
Step 1. Provide an initial guess and set .
Step 2. Solve the discretized state system in (3.14).
Step 3. Evaluate the reduced cost functional in (2.5).
Step 4. Use the result obtained in Step 2 to solve the discretized adjoint system in (3.15).
Step 5. Update the reduced gradient using (3.11).
Step 6. For , check the following stopping criteria: If one of them is met, stop.
Step 7. Compute the conjugate gradient update parameters and using (4.4).
Step 8. Compute the search direction using (4.3).
Step 9. Compute the steplength that satisfies condition in (4.2).
Step 10. Update the control variable using (4.1). Set and go to Step 2.
5. Numerical Experiment
In this section, we present the numerical experiment for the optimal control problem of monodomain model. The experiment setup is presented first, followed by the experiment results which are solved by the HSDY and LSCD methods.
5.1. Experiment Setup
The numerical experiment is carried out on a twodimensional computational domain of size cm^{2}, that is, and the final simulation time is set to be . Figure 1 displays the positions of the subdomains in the computational domain . From Figure 1, and are the control domains, and are the neighborhoods of the control domains, is the observation domain, and is the excitation domain. For domain discretization, the computational domain is discretized into 8192 triangular elements, with 7936 of them compose the observation domain , 144 of them compose the control domains and , and the rest of them compose the neighborhoods of the control domains and .
Table 1 lists the parameters that are used in our numerical experiment, with some of them adopted from Colli Franzone et al. [24]. Lastly, the initial values for the state and control variables are given as the following:

5.2. Experiment Results
In this section, we present the experiment results for the optimal control problem of monodomain model. The minimum values of the reduced cost functional along the optimization process for the HSDY and LSCD methods are depicted in Figures 2 and 3, respectively. Notice that the logarithmic scales are used in Figures 2 and 3 for clear presentation on how the minimum values of are decreased during the optimization process.
As shown in Figure 2, the HSDY method successfully converges to the optimal solution by taking 704 optimization iterations. However, this is not the case for the HS method. For the HS method, it failed to converge to the optimal solution and stopped at the 2nd iteration. This phenomenon happens because the search direction generated by the HS method may not be a descent direction and its global convergence is not guaranteed when the Armijo line search is used. These experiment results indicate that the HSDY method outperforms the HS method. On the other hand, both LSCD and LS methods successfully located the optimal solution by taking 700 and 703 optimization iterations, as shown in Figure 3. Again, the hybrid conjugate gradient method (LSCD method) performs better than the classical conjugate gradient method (LS method). Thus, we can conclude that the hybrid conjugate gradient methods are not only globally convergent but also superior to the classical conjugate gradient methods.
Figure 4 illustrates the corresponding norm of reduced gradient for the HSDY and LSCD methods. From the figure, the gradient for the HSDY method decreased sharply at the 8th iteration, followed by a smooth decrease till the end of optimization iterations. In contrast, the gradient for the LSCD method decreased less sharply than the HSDY method at the 8th iteration, and it finally approaches zero with a smooth decrease. Observe that the gradients for both HSDY and LSCD methods are almost the same starting from the 100th iteration to the end of optimization iterations.
Next, the uncontrolled solutions at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms are illustrated in Figure 5. Note that the uncontrolled solutions are obtained where no extracellular current is applied to the computational domain. As shown in Figure 5, the excitation wavefront spreads from the inside to the outside of the computational domain during the time interval from 0 ms to 2 ms. These experiment results imply that the excitation wavefront will continue to spread to the computational domain if the control (the extracellular current) is not switched on.
(a)
(b)
(c)
(d)
Figure 6 illustrates the optimally controlled solutions at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms using the HSDY and LSCD methods. For the optimally controlled case, the excitation wavefront is successfully dampened out by the optimal applied extracellular current . Also, observe that the excitation wavefront is almost completely dampened out at time 1.5 ms.
(a)
(b)
(c)
(d)
6. Conclusion
In this paper, we have presented the numerical experiment results for the optimal control problem of monodomain model using the hybrid conjugate gradient methods, namely, the HSDY and LSCD methods. Our experiment shows that both HSDY and LSCD methods successfully converge to the optimal solution. By comparison to the classical conjugate gradient methods, both HSDY and LSCD methods show their superiority over HS and LS methods in terms of global convergence as well as number of optimization iterations. We therefore conclude that the hybrid conjugate gradient methods outperform the classical conjugate gradient methods and are suitable for solving the optimal control problem of monodomain model owing to their low memory requirement and simple computation.
Acknowledgments
The work is financed by Zamalah Scholarship and Research University Grant (GUP) Vote 06J26 provided by Universiti Teknologi Malaysia (UTM) and the Ministry of Higher Education (MOHE) of Malaysia.
References
 V. Akcelik, G. Biros, and O. Ghattas, “Parallel multiscale GaussNewtonKrylov methods for inverse wave propagation,” in Proceedings of the IEEE/ACM SC2002 Conference, pp. 1–15, November 2002. View at: Google Scholar
 M. Fisher, J. Nocedal, Y. Trémolet, and S. J. Wright, “Data assimilation in weather forecasting: a case study in PDEconstrained optimization,” Optimization and Engineering, vol. 10, no. 3, pp. 409–426, 2009. View at: Publisher Site  Google Scholar
 O. Schenk, M. Manguoglu, A. Sameh, M. Christen, and M. Sathe, “Parallel scalable PDEconstrained optimization: antenna identification in hyperthermia cancer treatment planning,” Computer Science, vol. 23, no. 34, pp. 177–183, 2009. View at: Publisher Site  Google Scholar
 C. E. Orozco and O. N. Ghattas, “Massively parallel aerodynamic shape optimization,” Computing Systems in Engineering, vol. 3, no. 1–4, pp. 311–320, 1992. View at: Google Scholar
 S. B. Hazra and V. Schulz, “Simultaneous pseudotimestepping for aerodynamic shape optimization problems with state constraints,” SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 1078–1099, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Benzi, E. Haber, and L. Taralli, “A preconditioning technique for a class of PDEconstrained optimization problems,” Advances in Computational Mathematics, vol. 35, no. 2–4, pp. 149–173, 2011. View at: Publisher Site  Google Scholar
 E. Haber and U. M. Ascher, “Preconditioned allatonce methods for large, sparse parameter estimation problems,” Inverse Problems, vol. 17, no. 6, pp. 1847–1864, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 T. Rees and M. Stoll, “Blocktriangular preconditioners for PDEconstrained optimization,” Numerical Linear Algebra with Applications, vol. 17, no. 6, pp. 977–996, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. Biros and O. Ghattas, “Parallel LagrangeNewtonKrylovSchur methods for PDEconstrained optimization. Part I: the KrylovSchur solver,” SIAM Journal on Scientific Computing, vol. 27, no. 2, pp. 687–713, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. M. Shuaiby, M. A. Hassan, and M. ElMelegy, “Modeling and simulation of the action potential in human cardiac tissues using finite element method,” Journal of Communications and Computer Engineering, vol. 2, no. 3, pp. 21–27, 2012. View at: Google Scholar
 K. W. Ng and A. Rohanin, “Numerical solution for PDEconstrained optimization problem in cardiac electrophysiology,” International Journal of Computer Applications, vol. 44, no. 12, pp. 11–15, 2012. View at: Google Scholar
 Y. Belhamadia, A. Fortin, and Y. Bourgault, “Towards accurate numerical method for monodomain models using a realistic heart geometry,” Mathematical Biosciences, vol. 220, no. 2, pp. 89–101, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 K. W. Ng and A. Rohanin, “The effects of control domain size on optimal control problem of monodomain model,” International Journal of Computer Applications, vol. 47, no. 10, pp. 6–11, 2012. View at: Google Scholar
 C. Nagaiah, K. Kunisch, and G. Plank, “Numerical solution for optimal control of the reactiondiffusion equations in cardiac electrophysiology,” Computational Optimization and Applications, vol. 49, no. 1, pp. 149–178, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 E. Polak and G. Ribière, “Note sur la convergence de méthodes de directions conjuguées,” vol. 16, pp. 35–43, 1969. View at: Google Scholar  Zentralblatt MATH
 B. T. Polyak, “The conjugate gradient method in extremal problems,” USSR Computational Mathematics and Mathematical Physics, vol. 9, no. 4, pp. 94–112, 1969. View at: Google Scholar
 W. W. Hager and H. Zhang, “A new conjugate gradient method with guaranteed descent and an efficient line search,” SIAM Journal on Optimization, vol. 16, no. 1, pp. 170–192, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Y. H. Dai and Y. Yuan, “A nonlinear conjugate gradient method with a strong global convergence property,” SIAM Journal on Optimization, vol. 10, no. 1, pp. 177–182, 1999. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. M. Rogers and A. D. McCulloch, “A collocationGalerkin finite element model of cardiac action potential propagation,” IEEE Transactions on Biomedical Engineering, vol. 41, no. 8, pp. 743–757, 1994. View at: Publisher Site  Google Scholar
 K. Kunisch and M. Wagner, “Optimal control of the bidomain system (I): the monodomain approximation with the RogersMcCulloch model,” Nonlinear Analysis: Real World Applications, vol. 13, no. 4, pp. 1525–1550, 2012. View at: Publisher Site  Google Scholar
 Z. Qu and A. Garfinkel, “An advanced algorithm for solving partial differential equation in cardiac conduction,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 9, pp. 1166–1168, 1999. View at: Publisher Site  Google Scholar
 Y. H. Dai and Y. Yuan, “An efficient hybrid conjugate gradient method for unconstrained optimization,” Annals of Operations Research, vol. 103, no. 1–4, pp. 33–47, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Zhou, Z. Zhu, H. Fan, and Q. Qing, “Three new hybrid conjugate gradient methods for optimization,” Applied Mathematics, vol. 2, no. 3, pp. 303–308, 2011. View at: Publisher Site  Google Scholar
 P. Colli Franzone, P. Deuflhard, B. Erdmann, J. Lang, and L. F. Pavarino, “Adaptivity in space and time for reactiondiffusion systems in electrocardiology,” SIAM Journal on Scientific Computing, vol. 28, no. 3, pp. 942–962, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
Copyright
Copyright © 2012 Kin Wei Ng and Ahmad Rohanin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.