A Regularization Homotopy Strategy for the Constrained Parameter Inversion of Partial Differential Equations

The main difficulty posed by the parameter inversion of partial differential equations lies in the presence of numerous local minima in the cost function. Inversion fails to converge to the global minimum point unless the initial estimate is close to the exact solution. Constraints can improve the convergence of the method, but ordinary iterative methods will still become trapped in local minima if the initial guess is far away from the exact solution. In order to overcome this drawback fully, this paper designs a homotopy strategy that makes natural use of constraints. Furthermore, due to the ill-posedness of inverse problem, the standard Tikhonov regularization is incorporated. The efficiency of the method is illustrated by solving the coefficient inversion of the saturation equation in the two-phase porous media.


Introduction
Inverse problems of partial differential equations arise in a variety of practical problems in science and engineering. These range from biomedical and geophysical imaging to groundwater flow modeling [1][2][3][4][5][6][7][8][9]. It is very significant to conduct research into the theory of inverse problems and their applications. Generally, it is very difficult to solve inverse problems. This is due in part to the ill-posedness of the problem itself, which causes the solution not to depend continuously on the measurement data. It is also due to the presence of numerous local minima in the cost function, which hinders the convergence of the ordinary numerical methods. In order to overcome these two difficulties, a homotopy strategy that makes natural use of constraints is designed for the general parameter inversion of partial differential equations.
The homotopy method is a novel and effective method which has been successfully applied to find solutions of various nonlinear problems such as zeros and fixed points of mappings and so on. A remarkable advantage of this method is that the algorithm generated by it is globally convergent under certain weak assumptions [10]. Recently, it has also been extended to dealing with inverse problems. Keller and Perozzi [11] may be the first researchers to have used the homotopy method to solve the inverse problems. Then, Vasco [12] combined a homotopy path-following formalism with the singularity theory to track multiple solutions of a geophysical inverse problem. Shidfar et al. [13] used a weighted homotopy analysis method to solve the inverse problem of identifying an unknown source term in a parabolic equation. Zhao et al. [14] developed an adaptive homotopy method for the parameter identification inverse problem of a nonlinear diffusion equation. Hu et al. [15] constructed a homotopy approach to improve the PEM identification of ARMAX models. Cao and Han [16] presented the convergence analysis of the homotopy perturbation method for solving nonlinear ill-posed operator equations.
Generally, a parameter inversion identifies parameters using only measurement data recorded on the surface of the object to be measured. However, these measurement data may have a low signal-to-noise ratio. To suppress the noise and improve the quality of inversion, constraints have been widely used in the inversion fields, such as remote sensing of the environment [17,18], atmospheric research [19,20], geological exploration [21,22], petrophysics [23,24] and so on. The main reason for this is that the constraint data, which are recorded from the interior of the object to be measured, may be less noisy than the surface data.
In this paper, we study the parameter inversion problem of partial differential equations considering not only the measurement data but also the constraint data. Firstly, the constrained parameter inversion problem is formulated as the constrained minimization problem; then, this constrained minimization problem is transformed into an unconstrained minimization problem by the use of the penalty function method. Secondly, the homotopy method is introduced and combined with Tikhonov regularization, a successive approximation method, to give a widely convergent inversion strategy. In general terms, the successive approximation method aims to reduce uncertainty, which could be measured as decreasing the entropy of the distribution of uncertain system parameters, at a computational cost that is analogous to energy supply. We call this method the regularization homotopy strategy. Finally, as a practical application, a distributed parameter inversion problem of the saturation equation in the two-phase porous media is solved.
The contents of this paper are outlined as follows. The constrained parameter inversion problem of partial differential equations is described in Section 2. A brief introduction of homotopy theory is given in Section 3. The regularization homotopy strategy is presented in Section 4. In order to test the method's effectiveness, numerical simulations are carried out for the parameter inversion of the saturation equation in the fractional flow formulation of the two-phase porous media flow equations in Section 5. Finally, in Section 6, some conclusions are given.

Inversion Model
Consider the following partial differential equation where x = (x 1 , x 2 , . . . , x n ) , Ω ⊂ R n is a bounded domain, u(x, t) is a sufficiently smooth function defined on Ω × (0, ∞), and L is the differential operator. The initial and boundary conditions are where ∂Ω is the boundary of Ω, and I and B are the initial condition operator and boundary condition operator, respectively. If p(x), ϕ(x), φ(x, t) are known, Equations (1)-(3) form the direct problem to determine u(x, t).
If p(x) is unknown, given the additional condition with Γ being a part of Ω and A the additional condition operator, Equations (1)-(4) form the general parameter inversion of partial differential equations.
In the practical problems, the parameter p(x) is usually in the discrete form. Let P = (p 1 , p 2 , . . . , p Q ) denote the parameters to be determined in the discrete form. As the solution u(x, t) of Equations (1)-(3) nonlinearly depends on p(x), a nonlinear operator K(P) can be defined as then, the parameter inversion is reduced to the output least squares problem: Denote the admissible set where p i 1 , p i 2 , . . . , p i d are obtained from the constraint data, and Then, the constrained parameter inversion can be transformed into the solution of the problem min P∈Σ K(P) 2 .
Let P = ( p i 1 , p i 2 , . . . , p i d ), where E is the extraction operator. It is obvious that EP − P = 0 holds for P ∈ Σ. Thus, the least squares problem (5) can be written, without constraint, as where α is a constraint parameter to determine the strength of the constraint. The solutions of Equations (5) and (6) are close to each other when α is large enough. Therefore, in order for solution of Equation (6) to approximate the solution of Equation (5) well, it is necessary that α is specified as large enough in the specific inversion process.

Homotopy Theory
A brief introduction of the homotopy theory, which is applicable to nonlinear operator equations such as K(P) = 0, is given in this section. For more details, interested readers are referred to the tutorial provided by Watson [10].
By adding a homotopy parameter κ ∈ [0, 1] and an artificially constructed function C(P) to the "target function" K(P), a homotopy map F(P, κ) is constructed so that F(P, 0) = C(P), F(P, 1) = K(P), where C(P) is usually chosen as a simple function, and we assume that the solution of C(P) = 0 is known. For example, we may choose C(P) = P − P 0 , and the known solution of C(P) = 0 is P 0 . P 0 can be chosen as the initial value of the whole computation. It is assumed that there is a curve P = P(κ), κ ∈ [0, 1] satisfying the homotopy equation Obviously, P(0) = P 0 , P(1) = P * , and P * is the solution of K(P) = 0. Thus, P(κ) is a continuous curve connecting P 0 with P * . If some numerical method is applied to trace this curve, the solution P * will eventually be reached.
Generally, there are two ways to trace the curve P(κ). The first method is on the basis of the hypothesis that F(P, κ) is differentiable with respect to P and κ. By differentiating each side of Equation (7), Because P satisfies the initial condition P(0) = P 0 , the problem of determining P(κ) is transformed into an initial value problem of ordinary differential Equation (8). An approximation of P * = P(1) can be obtained by using numerical integration. The second method is first to divide the interval [0, 1] into 0 = κ 0 < κ 1 < · · · < κ S = 1 and then to solve the operator equations sequentially: By induction, the solution P s of the sth equation is known, so P s can be used as the initial approximation of the (s + 1)th equation of Equation (9). P s is such a good approximation of P s+1 that we only use a locally convergent method to obtain P s+1 , as long as κ s+1 − κ s is small enough. In this article, the latter method is used as our path tracing method.

Homotopy Strategy
It is well known that the inverse problem in Equation (6) is ill-posed, so the regularization method must be employed: where P 0 and β are the initial estimate and regularization parameter, respectively. Obviously, Equation (10) is equivalent to the corresponding normal equation where * denotes the adjoint operator. For the sake of avoiding the impact of the second derivative, a successive linearization method is used to construct a basic iterative method. It is also interesting to use other methods, such as gradient [25], gradientless [26], minimization with constraints [27] and stochastic approaches [28]. Assume that the kth approximation P k of Equation (11) has been found. To compute the next approximation P k+1 , the linear function L k (P) = K (P k )(P − P k ) + K(P k ) is used instead of K(P), and the cost function in Equation (10) is replaced by

Its corresponding normal equation is
the solution of which is exactly the next approximation P k+1 ; that is, This iterative method has a fast convergence speed and good stability; however, it only has local convergence. As a matter of fact, Equation (12) is a variant of the iteratively regularized Gauss-Newton method [29] and has the same computational cost as the latter.
In an attempt to overcome the local convergence problem, the homotopy method is introduced to solve Equation (11). Consider the following fixed-point homotopy equation where κ ∈ [0, 1] is the homotopy parameter. In addition, Equation (13) can be rearranged as With the purpose of finding the solution of Equation (11), the interval [0, 1] is divided into 0 = κ 0 < κ 1 < · · · < κ S = 1, and for κ = κ s , the above basic iterative method is used to solve Equation (14) sequentially. The known solution P 0 of F(P, κ 0 ) can be regarded as the initial estimate of the next equation F(P, κ 1 ). Suppose that the solution P s of F(P, κ s ) has been found; the successive linearization method, which is performed with Equation (12), can be applied to F(P, κ s+1 ). Therefore, we can construct the iterative formula as follows: In reality, the final iterative result P S of Equation (15), which is the solution of F(P, 1), is none other than the regularized solution of Equation (11). This regularization homotopy algorithm is globally convergent to solve the constrained parameter inversion of partial differential equations. For the global convergence analysis of the homotopy method, see [16,30].
Let κ s = s S and s T ≡ 0; then, we have a simplified version of Equation (15): The iterative result P S of Equation (16) can be used as the initial estimate of Equation (12). Then, iteration using Equation (12) is repeated until the regularized solution of Equation (11) is found. Consequently, Equations (12) and (16) can form a stable, globally convergent and fast method, as described in Algorithm 1. The flowchart of Algorithm 1 is shown in Figure 1. Step 2, else if s ≥ S, then return P s . 4: Set k = 0, and P k = P s . 5: Compute P k+1 by Equation (12). 6: Let k = k + 1 and check if P k − P k−1 > ε, then go to Step 5, else if P k − P k−1 ≤ ε, then return P k .
If α = 0, Equation (10) is simply the ordinary parameter inversion of partial differential equations, and Equations (12) and (16) can respectively be rewritten as follows: and After the comparison of the proposed method and the other well-known iterative methods (such as the regularized Gauss-Newton method), it is found that their calculation amounts and storage requirements are the same at each step. The application of the firstorder derivative, the evaluation of the adjoint operator and the forward-modeling run are of primary importance; the most important factor is, however, that the convergence domain of the proposed method is much wider than the other iterative methods, such as the regularized Gauss-Newton method.

Numerical Simulations
To show the feasibility of the proposed method and illustrate its numerical behavior, the numerical simulation procedure is carried out for two concrete models in MATLAB R2018a environment. We choose 2 , s(x, y, t) = 0, N(u) = u 2 − u + 1, ϕ(x, y) = sin(πx) sin(πy), φ(x, y, t) = 0, ∆x = ∆y = 1/24, ∆t = 0.002, x 0 = 12/24, α = 10 4 , β = 10 −5 , S = 10, P 0 ≡ 5. The first true model is displayed in Figure 2. When there are 5%, 10%, 15% and 20% Gaussian noises in the measurement data, the inversion results of the homotopy method with constraints and the homotopy method without constraints are shown in Figures 3 and 4, respectively. When the noise level reaches 15% or 20%, the inversion result is still satisfactory for the homotopy method with constraints (see Figure 3c,d); however, the homotopy method without constraints is not convergent.  The second true model is the model of two anomalous bodies in a homogeneous medium, as shown in Figure 5. Figures 6 and 7 show the inversion results of the homotopy method with constraints and the homotopy method without constraints, respectively, with the noisy data at 5%, 10%, 15% and 20%. Figure 6c,d show the inversion results of the homotopy method with constraints when 15% and 20% Gaussian noises are added to the measurement data, respectively. If the same noise level is used, the homotopy method without constraints is divergent.
For comparison, we investigate three methods for the above two models: the homotopy method with constraints (Equations (12) and (16)), homotopy method without constraints (Equations (17) and (18)), and iterative method with constraints (Equation (12)). The relative errors of inversion results are tabulated in Table 1. Table 1 shows that with 5% and 10% Gaussian noises added, the inversion results by the homotopy method with constraints and the homotopy method without constraints both have very high accuracy; with 15% and 20% Gaussian noises added, the inversion results for the homotopy method with constraints are still satisfactory, but the inversion results for the homotopy method without constraints are not acceptable, which shows the necessity of the constraints; moreover, starting from the same initial estimate, the iteration of the iterative method with constraints is not convergent, which shows the necessity of introducing the homotopy strategy.  To compare the proposed method with other methods (wavelet multiscale method, nonlinear multigrid method, adaptive multigrid conjugate gradient method, homotopy perturbation method), we use all these methods for the first model with 5% Gaussian noise added. When the initial estimate P 0 ≡ 5, the relative errors and CPU run times of these methods are listed in Table 2; When the initial estimate P 0 ≡ 1, the relative errors and CPU run times of these methods are listed in Table 3. From Tables 1-3, it is easy to see the basic information about the advantages and disadvantages of the homotopy method with constraints compared to other methods, as presented in Table 4.
In order to study the sensibility of the inversion results on the parameters α, β, P 0 , we use the homotopy method with constraints for the first model with 15% Gaussian noise added, test several values of α, β, P 0 and list the relative errors in Table 5. From Table 5, it can be seen that the performance of the proposed method is highly sensitive to the values of α and β and insensitive to the value of P 0 . Usually, the optimal regularization parameter is unknown and must be determined by methods such as trial-and-error, cross-validation, L-curve criterion, the discrepancy principle and others. In this paper, we determine the parameters α and β by trial-and-error, for the sake of simplicity. From Table 5, we find that the constraint parameter α = 10 4 performs best, that the regularization parameter β = 10 −5 performs best and that the value of initial estimate P 0 does not affect the inversion result; thus, we let α = 10 4 , β = 10 −5 . It is because of the global convergence property of the homotopy method that the value of initial estimate P 0 can arbitrarily be chosen. In fact, the value of P 0 is chosen from range of 0 to 10, since the permeability values are from range of 0 to 10 in reality. The value of P 0 is usually chosen to be the mean value 5 of 0 to 10, when no prior knowledge is available about the permeability parameter.  Table 4. Advantages and disadvantages of the homotopy method with constraints compared to other methods.

Compared Method Advantage Disadvantage
Homotopy method without constraints better anti-noise ability × Iterative method with constraints larger convergence region × Wavelet multiscale method better anti-noise ability and larger convergence region more computation Nonlinear multigrid method better anti-noise ability and larger convergence region more computation Adaptive multigrid conjugate gradient method better anti-noise ability and larger convergence region more computation Homotopy perturbation method better anti-noise ability × In order to pictorially demonstrate how the homotopy method avoids the local minimum, we use the homotopy method with constraints and iterative method with constraints for the first model with no noise. Three locations ( 4 24 , 4 24 ), ( 4 24 , 20 24 ), ( 20 24 , 4 24 ) of the parameter model are selected as the (x, y, z) coordinates of three-dimensional plotting. The true parameter values of these locations are 1, 2 and 3, respectively, so the global minimum is (1, 2, 3). P 0 ≡ 5, so the initial guess is (5,5,5). Figure 8 pictorially shows how the homotopy method with constraints converges to the global minimum from the initial guess (5,5,5) and how the iterative method with constraints converges to the local minimum from the initial guess (5, 5, 5).  Submarine oil exploration is one of practical applications of the proposed solution. When reservoir workers value a block and plan to extract oil from the ocean floor, they first need to use methods such as seismic wave inversion and geophysical exploration to show where the oil and gas reservoirs may be. After choosing a location for oil extraction, they starts drilling wells, as shown in Figure 9. Well log data from the drilled well can be used as the constraint data. Based on the complementary characteristics of seismic data and logging data, the joint inversion technology of seismic data and logging data is studied. The regularization homotopy method proposed in this paper is combined with logging constraints to jointly invert the permeability parameters of the nonlinear convection-diffusion equation to obtain the complete distribution information of permeability parameters and further improve the inversion accuracy and anti-noise ability.

Conclusions
A regularization homotopy strategy is designed for the constrained parameter inversion of partial differential equations. As a practical application, this method has been used successfully to solve the constrained coefficient inversion of the saturation equation in the fractional flow formulation of the two-phase porous media flow equations. It turns out that the homotopy strategy can effectively widen the convergence region of the ordinary numerical methods, and the constraints can effectively improve the noise suppression ability of the inversion methods. Future work will extend this new strategy to the ill-posed image processing problems and the optimal control problems with nonlinear dynamical systems as constraint conditions. Author Contributions: Conceptualization, T.L.; methodology, T.L. and R.X.; software, T.L.; validation, T.L., R.X., C.L. and Y.Q.; investigation, T.L. and R.X.; resources, Y.Q.; writing-original draft preparation, R.X.; writing-review and editing, T.L. and R.X.; supervision, C.L.; funding acquisition, T.L. and Y.Q. All authors have read and agreed to the published version of the manuscript.