Solving Optimal Control Problem of Monodomain Model Using Hybrid Conjugate Gradient Methods

Abstract

We present the numerical solutions for the PDE-constrained 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 Hestenes-Stiefel-Dai-Yuan (HS-DY) method and the Liu-Storey-Conjugate-Descent (LS-CD) 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 PDE-constrained optimization problem. PDE-constrained 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 PDE-constrained 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 PDE-constrained 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, Polak-Ribière-Polyak (PRP) method [15, 16], the Hager-Zhang (HZ) method [17], and the Dai-Yuan (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 Rogers-modified FitzHugh-Nagumo 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 𝐻=Ω×[0,𝑇] and 𝜕𝐻=𝜕Ω×[0,𝑇]. The optimal control problem of monodomain model with Rogers-modified FitzHugh-Nagumo ion kinetic is therefore given by the following: min𝐽𝑉,𝑤,𝐼_𝑒=12^𝑇_0Ω𝑜|𝑉|²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐𝜆𝑑𝑡s.t.𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒=0,in𝐻𝜕𝑤()𝐷𝜕𝑡−𝑓𝑉,𝑤=0,in𝐻_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻𝑥,0=𝑉⁰,𝑤(𝑥,0)=𝑤⁰,onΩ, (2.1)(2.2) where 𝐼_ion(𝑉,𝑤)=𝑐₁𝑉𝑉1−𝑉_th𝑉1−𝑉_𝑝+𝑐₂𝑓()𝑤𝑉,𝑉,𝑤=𝑐₃𝑉𝑉_𝑝−𝑐₄𝑤. (2.3)(2.4) 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 surface-to-volume ratio of the cell membrane, 𝐶_𝑚 is the membrane capacitance, 𝐼_ion(𝑉,𝑤) is the current density flowing through the ionic channels, 𝑤(𝑥,𝑡) is the ionic current variable, 𝑓(𝑉,𝑤) is the prescribed vector-value function, 𝑉_th 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 Rogers-modified FitzHugh-Nagumo model [19] for representing ion kinetics.

According to Kunisch and Wagner [20], the control-to-state 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: min𝐽̂𝐼_𝑒=12^𝑇_0Ω𝑜𝑉𝐼_𝑒²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐𝑑𝑡, (2.5) where (2.5) is known as the reduced cost functional.

3. Discretization of the Optimal Control Problem

We adopt the optimize-then-discretize 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. First-Order Optimality System

For deriving the infinite dimensional optimality system, the Lagrange functional ℒ is formed as follows: 1ℒ=2^𝑇_0Ω𝑜𝑉𝐼_𝑒²𝑑Ω_𝑜+𝛼_Ω𝑐𝐼_𝑒²𝑑Ω_𝑐+𝑑𝑡_𝑇0Ω𝜆𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒+𝑝𝑑Ω𝑑𝑡_𝑇0Ω𝜕𝑤()𝜕𝑡−𝑓𝑉,𝑤𝑞𝑑Ω𝑑𝑡, (3.1) where 𝑝(𝑥,𝑡) and 𝑞(𝑥,𝑡) are the adjoint variables used to adjoin the constraints in (2.2) to the reduced cost functional in (2.5). The first-order 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: ℒ_𝑉∶𝑉|_𝑜+𝜆𝐷1+𝜆∇⋅_𝑖∇𝑝+𝛽𝐶_𝑚𝜕𝑝𝐼𝜕𝑡−𝛽_ion𝑉𝑝−𝑓_𝑉ℒ𝑞=0_𝑤∶−𝛽𝐼_ion𝑤𝑝−𝜕𝑞−𝑓𝜕𝑡^𝑤ℒ𝑞=0_𝑝∶𝜆𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐶_𝑚𝜕𝑉𝜕𝑡−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒ℒ=0_𝑞∶𝜕𝑤()ℒ𝜕𝑡−𝑓𝑉,𝑤=0_𝐼𝑒∶𝛼𝐼_𝑒−11+𝜆𝑝=0, (3.2)(3.3)(3.4)(3.5)(3.6) 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: 𝐷_𝑖𝑝()()∇𝑝⋅𝜂=0,on𝜕𝐻𝑥,𝑇=0,𝑞𝑥,𝑇=0,onΩ. (3.7)(3.8)

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: 𝛽𝐶_𝑚𝜕𝑉=𝜆𝜕𝑡𝐷1+𝜆∇⋅_𝑖∇𝑉−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒,in𝐻𝜕𝑤()𝐷𝜕𝑡=𝑓𝑉,𝑤,in𝐻_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻𝑥,0=𝑉⁰,𝑤(𝑥,0)=𝑤⁰,onΩ, (3.9) and the adjoint system is given by the following: 𝛽𝐶_𝑚𝜕𝑝𝜆𝜕𝑡=−𝐷1+𝜆∇⋅_𝑖∇𝑝+𝛽𝐼_ion𝑉𝑝+𝑓_𝑉𝑞−𝑉|_𝑜,in𝐻𝜕𝑞𝐼𝜕𝑡=−𝛽_ion𝑤𝑝−𝑓_𝑤𝐷𝑞,in𝐻_𝑖𝑝()()∇𝑝⋅𝜂=0,on𝜕𝐻𝑥,𝑇=0,𝑞𝑥,𝑇=0,onΩ. (3.10) Also, by utilizing (3.6), the reduced gradient is given by the following: ∇𝐽̂𝐼_𝑒=𝛼𝐼_𝑒−11+𝜆𝑝. (3.11) 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 first-order 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 𝛽𝐶_𝑚𝜕𝑉=𝜆𝜕𝑡𝐷1+𝜆∇⋅_𝑖∇𝑉,in𝐻𝛽𝐶_𝑚𝜕𝑉𝜕𝑡=−𝛽𝐼_ion()−1𝑉,𝑤𝐼1+𝜆^𝑒,in𝐻𝜕𝑤()𝐷𝜕𝑡=𝑓𝑉,𝑤,in𝐻_𝑖𝑉()∇𝑉⋅𝜂=0,on𝜕𝐻𝑥,0=𝑉⁰𝑤(),onΩ𝑥,0=𝑤⁰,onΩ, (3.12) while the adjoint system in (3.10) becomes 𝛽𝐶_𝑚𝜕𝑝𝜆𝜕𝑡=−𝐷1+𝜆∇⋅_𝑖∇𝑝,in𝐻𝛽𝐶_𝑚𝜕𝑝𝐼𝜕𝑡=𝛽_ion𝑉𝑝+𝑓_𝑉𝑞−𝑉|_𝑜,in𝐻𝜕𝑞𝐼𝜕𝑡=−𝛽_ion𝑤𝑝−𝑓_𝑤𝐷𝑞,in𝐻_𝑖𝑝()𝑞()∇𝑝⋅𝜂=0,on𝜕𝐻𝑥,𝑇=0,onΩ𝑥,𝑇=0,onΩ. (3.13)

For the discretization procedure, the linear PDEs in (3.12) and (3.13) are discretized with the Crank-Nicolson 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: 𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛+1=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐕1+𝜆^𝑛,𝐕^𝑛+1=𝐕^𝑛+Δ𝑡₂−𝐈^𝑛_𝐢𝐨𝐧𝐶_𝑚−𝐈^𝑛_𝐞𝛽𝐶_𝑚(),𝐰1+𝜆_𝑛+1=𝐰^𝑛+Δ𝑡₂𝐟^𝑛,𝐕(𝑥,0)=𝐕⁰,𝐰(𝑥,0)=𝐰⁰. (3.14) Here Δ𝑡₁ and Δ𝑡₂ are the local time-steps, 𝐌 is the mass matrix, and 𝐊 is the stiffness matrix. On the other hand, the discretized adjoint system is given by the following: 𝛽𝐶_𝑚𝐌+Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛=𝛽𝐶_𝑚𝐌−Δ𝑡₁2𝜆𝐊𝐩1+𝜆^𝑛+1,𝐩^𝑛=𝐩^𝑛+1+Δ𝑡₂𝐕^𝑛+1_𝑜𝛽𝐶_𝑚−𝐈^𝑛+1_ion𝑉𝐶_𝑚𝐩^𝑛+1−𝐟^𝑛+1_𝑉𝛽𝐶_𝑚𝐪^𝑛+1,𝐪^𝑛=𝐪^𝑛+1+Δ𝑡₂𝛽𝐈^𝑛+1_ion𝑤𝐩^𝑛+1+𝐟^𝑛+1_𝑤𝐪^𝑛+1,𝐩(𝑥,𝑇)=𝟎,𝐪(𝑥,𝑇)=𝟎. (3.15)

4. Optimization Procedure

Two hybrid conjugate gradient methods are used to solve the discretized optimal control problem of monodomain model, namely, the Hestenes-Stiefel-Dai-Yuan (HS-DY) method [22] and the Liu-Storey-Conjugate-Descent (LS-CD) method [23]. These hybrid conjugate gradient methods combine the good numerical performance of the Hestenes-Stiefel (HS) and Liu-Storey (LS) methods with the strong convergence properties of the Dai-Yuan (DY) and Conjugate-Descent (CD) methods.

4.1. Hybrid Conjugate Gradient Methods

Starting from an initial guess 𝐈⁰_𝐞, our control is updated using the following recurrence: 𝐈_𝐞^𝑘+1=𝐈^𝑘_𝐞+𝛿^𝑘𝐝^𝑘,𝑘=0,1,…, (4.1) where 𝛿^𝑘>0 is the step-length computed by the Armijo line search and 𝐝^𝑘 is the search direction. Given an initial step-length 𝛿>0 and 𝜇,𝜌∈(0,1), the Armijo line search chooses 𝛿^𝑘=max{𝛿,𝛿𝜌,𝛿𝜌²,…} such that 𝐽̂𝐈^𝑘_𝐞+𝛿^𝑘𝐝^𝑘≤𝐽̂𝐈^𝑘_𝐞+𝜇𝛿^𝑘∇𝐽̂𝐈^𝑘_𝐞^𝑇𝐝^𝑘. (4.2) On the other hand, the search direction is defined by 𝐝^𝑘=−∇𝐽̂𝐈^𝑘_𝐞𝐽̂𝐈,if𝑘=0−∇^𝑘_𝐞+𝜃^𝑘_∗𝐝^𝑘−1,if𝑘>0, (4.3) where 𝜃^𝑘_∗∈ℜ is the conjugate gradient update parameter. The conjugate gradient update parameters for the HS-DY and LS-CD methods are given as follows: 𝜃^𝑘_HS-DY=max0,min𝜃^𝑘_HS,𝜃^𝑘_DY𝜃^𝑘_LS-CD=max0,min𝜃^𝑘_LS,𝜃^𝑘_CD, (4.4) where 𝜃^𝑘_HS=∇𝐽̂𝐈^𝑘_𝐞^𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1𝐝^𝑘−1𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1𝜃^𝑘_DY=∇𝐽̂𝐈^𝑘_𝐞²𝐝^𝑘−1𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1𝜃^𝑘_LS∇𝐽̂𝐈=−^𝑘𝐞𝑇∇𝐽̂𝐈^𝑘_𝐞−∇𝐽̂𝐈_𝐞^𝑘−1𝐝^𝑘−1𝑇∇𝐽̂𝐈_𝐞^𝑘−1𝜃^𝑘_CD∇𝐽̂𝐈=−^𝑘𝐞2𝐝^𝑘−1𝑇∇𝐽̂𝐈_𝐞^𝑘−1. (4.5)

4.2. Optimization Algorithm

In this section, we present the optimization algorithm for solving the discretized optimal control problem of monodomain model using the HS-DY and LS-CD methods. The optimization algorithm for these two hybrid conjugate gradient methods is therefore given as follows.

Step 1. Provide an initial guess 𝐈⁰_𝐞 and set 𝑘=0.

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 𝑘≥1, check the following stopping criteria: 𝐽̂𝐈^𝑘_𝐞−𝐽̂𝐈_𝐞^𝑘−1≤10⁻³∇𝐽̂𝐈^𝑘_𝐞≤10⁻³1+𝐽̂𝐈^𝑘_𝐞. (4.6) If one of them is met, stop.

Step 7. Compute the conjugate gradient update parameters 𝜃^𝑘_HS-DY and 𝜃^𝑘_LS-CD using (4.4).

Step 8. Compute the search direction 𝐝^𝑘 using (4.3).

Step 9. Compute the step-length 𝛿^𝑘 that satisfies condition in (4.2).

Step 10. Update the control variable 𝐈_𝐞^𝑘+1 using (4.1). Set 𝑘=𝑘+1 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 HS-DY and LS-CD methods.

5.1. Experiment Setup

The numerical experiment is carried out on a two-dimensional computational domain of size 1×1 cm², that is, Ω=[0,1]×[0,1] and the final simulation time is set to be 𝑇=2ms. Figure 1 displays the positions of the subdomains in the computational domain Ω. From Figure 1, Ω_𝑐1 and Ω_𝑐2 are the control domains, Ω_𝑐1 and Ω_𝑐2 are the neighborhoods of the control domains, Ω_𝑜=Ω⧵(Ω_𝑐1∪Ω_𝑐2) is the observation domain, and Ω_exi⊂Ω_𝑜 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 Ω_𝑐1 and Ω_𝑐2, and the rest of them compose the neighborhoods of the control domains Ω_𝑐1 and Ω_𝑐2.

Figure 1 

Computational domain Ω and its subdomains.

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: 𝑉(𝑥,0)=𝑉⁰=105mV,𝑥∈Ω_exi𝑤()0mV,otherwise𝑥,0=𝑤⁰𝐼=0,𝑥∈Ω_𝑒(𝑥,0)=𝐼⁰_𝑒=0mAcm⁻³,𝑥∈Ω_𝑐0mAcm⁻³,otherwise. (5.1)

Table 1 

Parameters used in the numerical experiment.

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 HS-DY and LS-CD 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.

Figure 2 

Minimum values of reduced cost functional 𝐽̂ for HS-DY and HS methods.

Figure 3 

Minimum values of reduced cost functional 𝐽̂ for LS-CD and LS methods.

As shown in Figure 2, the HS-DY 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 HS-DY method outperforms the HS method. On the other hand, both LS-CD 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 (LS-CD 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 HS-DY and LS-CD methods. From the figure, the gradient for the HS-DY method decreased sharply at the 8th iteration, followed by a smooth decrease till the end of optimization iterations. In contrast, the gradient for the LS-CD method decreased less sharply than the HS-DY method at the 8th iteration, and it finally approaches zero with a smooth decrease. Observe that the gradients for both HS-DY and LS-CD methods are almost the same starting from the 100th iteration to the end of optimization iterations.

Figure 4 

Norm of reduced gradient ‖∇𝐽‖̂ for HS-DY and LS-CD methods.

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 5 

The uncontrolled solutions (𝑉) at (a) 0.2 ms; (b) 0.8 ms; (c) 1.5 ms; and (d) 2.0 ms.

Figure 6 illustrates the optimally controlled solutions at times 0.2 ms, 0.8 ms, 1.5 ms, and 2.0 ms using the HS-DY and LS-CD methods. For the optimally controlled case, the excitation wavefront is successfully dampened out by the optimal applied extracellular current 𝐈^opt_𝐞. Also, observe that the excitation wavefront is almost completely dampened out at time 1.5 ms.

(a)

(b)

(c)

(d)

Figure 6 

The controlled solutions (𝑉opt) at (a) 0.2 ms; (b) 0.8 ms; (c) 1.5 ms; and (d) 2.0 ms.

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 HS-DY and LS-CD methods. Our experiment shows that both HS-DY and LS-CD methods successfully converge to the optimal solution. By comparison to the classical conjugate gradient methods, both HS-DY and LS-CD 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.