Numerical methods for the simulation
Many biological systems of interest exhibit a nonlinear dynamic behavior which makes the analytical solution of models representing such systems rather complicated, if not impossible, for most of the realistic situations. In addition to nonlinearity, these processes may present a spatially distributed nature. As a consequence they must be described using PDEs which, in turns, makes the analytical approach even more difficult. Numerical techniques must be, therefore, employed to solve the model equations.
Most of numerical methods employed for solving PDEs, in particular those employed in this work, belong to the family of methods of weighted residuals
in which the solution of the PDE system (2) is approximated by a truncated Fourier series of the forma
Depending on the selection of the basis functions ϕ
(ξ) different methodologies arise. In this work, two groups will be considered: those using locally defined basis functions as it is the case in classical techniques like the finite difference method or the finite element method and those using globally defined basis functions.
Methods using local basis functions
The underlying idea is to discretize the domain of interest into a (usually large) number N of smaller subdomains. In these subdomains local basis functions, for instance low order polynomials, are defined and the original PDE is approximated by N ordinary differential equations (ODE). The shape of the elements and the type of local functions allow distinguishing among different alternatives.
Probably the most widely used approaches for this transformation are the finite difference and the finite element methods. The reader interested on an extensive description of these techniques is referred to the literature[29–31]. Both of these methods have been successfully applied in the context of dynamic optimization[19, 32].
However it must be highlighted that in many biological models, especially those in 2D and 3D, the number of discretization points (N) to obtain a good solution might be too large for their application in optimization. Methods using global basis functions, which will reduce the computational effort, constitute an efficient alternative.
Methods using global basis functions
Different techniques like the eigenfunctions obtained from the Laplacian operator, Chevyshev or Legendre polynomials, among others have been considered over the last decades - see[33
] and references therein for a detailed discussion -. Probably the most efficient order reduction technique is the proper orthogonal decomposition
] and because of this, it will be chosen in this work to obtain the reduced order models. In this approach each element ϕ
) of the set of basis functions in (8) is computed off-line as the solution of the following integral eigenvalue problem[34
corresponds with the eigenvalue associated with each global eigenfunction ϕ
. The kernel R
( ξ ξ
) in equation (9
) corresponds with the two point spatial correlation function, defined as follows:
with x( ξ t
) denoting the value of the field at each instant t
and the summation extends over a sufficiently rich collection of uncorrelated snapshots at j = 1,⋯, ℓ. The basis functions obtained by means of the POD technique are also known as empirical basis functions or POD basis.
The dissipative nature of this kind of systems makes that the eigenvalues obtained from Eqn (9) can be ordered so that λ
. This property allows us to define a finite (usually low) dimensional subset ϕ
which captures the relevant features of the system[35
]. The number of elements (N
) in this subset is usually chosen using a criteria based on the energy captured by the POD basis. Such energy is connected to the eigenspectrum
or, to be more precise, to the inverse of the eigenvalues
In order to compute the time dependent coefficients in Eqn (8), the original PDE system (2) is projected onto each element of the POD basis set. In this particular case, such projection is carried out by multiplying the original PDE by each ϕ
and integrating the result over the spatial domain, this is:
Substituting the Fourier series approximation (8) into Eqn (12) leads to:
The basis functions obtained from (9) are orthogonal and can be normalized so that:
Therefore, Eqn (13) can be rewritten as:
is a row vector of the form
corresponds with the following column vector
. Expression (14) can be rewritten in matrix form as follows:
], and. Initial conditions for solving Eq (15) are obtained by projecting the original initial conditions x( ξ,0) over the basis functions, this is. At this point the basis functions ϕ
are known from Eq (9) while time dependent coefficients are computed by solving Eq (15), therefore the approximation of the original field x can be recovered by applying Eqn (8), this is. It is important to highlight that the number of elements N in the basis subset ϕA can be increased to approximate the original state x with an arbitrary degree of accuracy.
Dynamic optimization methods
There are several alternatives for the solution of dynamic optimization problems from which the direct methods are the most widely used. These methods transform the original problem into a non-linear programming (NLP) problem by means of complete parameterization, multiple shooting or control vector parameterization methods. Basically, all of them are based on the use of some type of discretization and approximation of either the control variables or both the control and state variables. The three alternatives basically differ in: the resulting number of decision variables, the presence or absence of parameterization related constraints and the necessity of using an initial value problem solver.
While the complete parameterization or the multiple shooting approaches may become prohibitively expensive in computational terms, the control vector parameterization approach allows handling large-scale dynamic optimization problems, such as those related to PDE systems, without solving very large NLPs and without dealing with extra junction constraints.
The control vector parameterization proceeds dividing the process duration into a number of elements and approximating the control functions typically using low order polynomials. The polynomial coefficients become the new decision variables and the solution of the resulting NLP problem (outer iteration) involves the system dynamics simulation (inner iteration).
Nonlinear programming methods may be largely classified in two groups: local and global methods. Local methods are designed to generate a sequence of solutions, using some type of pattern search or gradient and Hessian information that will converge to a local optimum. However the NLP arising from the application of the control vector parameterization method are frequently multimodal (i.e. presenting multiple local optima), due to the highly nonlinear nature of the dynamics. In this scenario, the initial guess may be located in the basin of attraction of a local minimum. This may be easily assessed by solving the problem from different initial guesses (multistart). In fact, this may be regarded as the first global optimization strategy. However experience demonstrates that there is no guarantee of arriving to the global solution, even starting from a large number of different initial guesses, and becomes computationally too expensive as illustrated in the examples considered in[10, 42] and later in this work.
Over the last decade a number of researchers have proposed different techniques for the solution of multimodal optimization problems. Depending on how the search is performed and which information they are exploiting the alternatives may be classified in two major groups: deterministic and stochastic.
Global deterministic methods in general take advantage of the problem’s structure and guarantee global convergence for some particular problems that verify specific smoothness and differentiability conditions. A number of works have recently approached the solution of dynamic optimization problems using convex relaxations or branch-and-bound strategies[42, 44, 45]. Although very promising, the necessary conditions for these methods to be applicable may not be guaranteed for the cases of interest and the computational cost may become prohibitive, particularly as the number of decision variables and the simulation cost increase.
The main drawbacks of global deterministic methods have motivated the use of stochastic methods that do not require any assumptions about the problem’s structure. They make use of pseudo-random sequences to determine search directions toward the global optimum. This leads to an increasing probability of finding the global optimum during the run time of the algorithm, although convergence may not be guaranteed. The main advantage of these methods is that, in practice, they rapidly arrive to the proximity of the solution.
The most successful approaches lie in one (or more) of the following groups: pure random search and adaptive sequential methods, clustering methods or metaheuristics. Metaheuristics are a special class of stochastic methods which have proved to be very efficient in recent years. They include both population (e.g., genetic algorithms) or trajectory-based (e.g., simulated annealing) methods. They can be defined as guided heuristics and many of them try to imitate the behavior of natural or social processes that seek for any kind of optimality. Some of these strategies have been successfully applied to the dynamic optimization of bioprocesses.
Despite the fact that many stochastic methods can locate the vicinity of global solutions very rapidly, the computational cost associated to the refinement of the solution is usually very large. In order to surmount this difficulty, hybrid methods and metaheuristics that have been recently developed which combine global stochastic methods with local gradient based methods in two phases or in several phases as in the scatter search based method eSS[23, 48].
Finally, knowing that global optimization methods become prohibitively expensive with an increasing number of decision variables, a control refining technique has been used so as to obtain smoother control profiles. This technique consists of performing successive re-optimizations with increasing control discretization level. A detailed description of the mesh refining approach used is presented in. The main steps are the following:
Step 1: The problem is solved using a coarse control discretization level (for example, 5−10) with the hybrid optimization method.
Step 2: The best solution found is transformed by multiplying the discretization level by for example 2−4 and the result is employed as the starting point for the local method.
Step 3: Step 2 is repeated until the established number of refinements has been achieved.