The problem of determining parameters of the basic hydrophysical characteristics is studied. These parameters are defined by the type of the soil. To determine these parameters, a model of unsaturated water flow in porous media is considered. The modeled values of soil moisture at various depths are obtained as a result of solution of the initial boundary values problem for Richards equation. The parameters identification problem is stated as an optimal control problem. The objective function is mean-square deviation of simulated values of soil moisture at various points from some prescribed values. Discretized problem is proposed to solve by Marquardt-Levenberg method.
parameters. In [Takeshita, 1999], [Vrugt et al., 2001] the parameters were obtained by genetic algorithm. Over the past decades, to determine the parameters, many authors have applied the methods imitating the behavior of biological populations in conditions of lack of vital resources and migrating in order to find a place with favorable living conditions, algorithms imitating social behavior, see, for example, [Abbaspourt et al., 2001], [Yang & You, 2013].
In the present paper the parameter identification problem is stated as an optimal control problem in which the control is unknown parameters, and the objective function is mean-square deviation of calculated values of soil moisture at various depths from some prescribed values. Calculation of soil moisture is performed in according to the model of water transfer in soil. As a result of finite difference approximation, the problem is reduced to nonlinear programming problem. Numerical solution is obtained by Marquardt-Levenberg method. Jacobian of the moisture function is calculated by formulas of fast automatic differentiation [Aida-Zade & Evtushenko, 1989], [Griewank & Corliss, 1999], [Evtushenko, 1991], [Evtushenko, 1998]. 2
Consider an one-dimensional model of vertical water transfer in soil. Suppose that soil is homogeneous isothermal non-deformable porous media. Under these assumptions vertical water transfer in soil is well described by onedimensional nonlinear second order parabolic partial differential equation. Consider following initial boundary value problem: ∂θ ∂ ( ∂θ ) ∂K(θ) ∂t = ∂z D(θ) ∂z − ∂z θ(z, 0) = φ(z), z ∈ (0, L), θ(L, t) = ψ(t), t ∈ (0, T ),
( ∂θ ) − D(θ) ∂z − K(θ)
z=0 θmin ≤ θ(0, t) ≤ θmax, t ∈ (0, T ), ,
(z, t) ∈ Q, = R(t) − E(t), t ∈ (0, T ), where z is space variable; t is time; θ(z, t) is soil moisture at the point (z, t); Q = (0, L) × (0, T ); φ(z) and ψ(t) are given functions; D(θ) and K(θ) are diffusion coefficient and hydraulic conductivity – the hydrophysical characteristics of the soil; θmin = θr + ε and θmax = θs − ε, where θr and θs are, respectively, the residual moisture and the saturation moisture depending on the soil type, and ε is a constant such that 0 < ε ≪ θr; R(t) is precipitation; E(t) is evaporation, 0 ≤ E(t) ≤ M , t ∈ (0, T ), M is a constant such that M > 0.
The diffusion coefficient D(θ) and the hydraulic conductivity K(θ) appearing in this equation are found by the widely used van Genuchten formulas [van Genuchten, 1980]
K(θ) = K0S0:5[1 − (1 − S1=m)m]2, D(θ) = K0 1 − m αm(θs − θr)
S0:5 1=m × [(1 − S1=m) m + (1 − S1=m)m − 2], (1) (2) θ − θr θs − θr where S =
; and K0, α, m are some parameters. Described problem (1)-(2) will be called the direct problem.
Formulate the parameters identification problem. Let a function θˆ(z, t) be defined on some set Q0 ⊆ Q. Call this function θˆ(z, t) ”experimental data”. Introduce a set U = {u : u ∈ R3; 0 ≤ a[i] ≤ u[i] ≤ b[i], i = 1, 3}. Denote [K0, α, m]T by u. The problem is to pick up the parameters K0, α and m in such a way that corresponding solution of the direct problem (1)-(2) is as close as possible to the function θˆ(z, t) on the set Q0. More precisely, the problem is to find uopt, uopt ∈ U , and corresponding solution θopt(z, t) of the direct problem (1)-(2) which minimize functional
J = θin+1 − θin = τ 1 (Din++11=2 θin++11 − θin+1 − Kin++11=2 − Dn+1 θin+1 − θin+11 + Kn+1 h h i 1=2 h i 1=2 ) , 1 ≤ i < I;
0 ≤ n < N, θi0 = φi, 0 ≤ i ≤ I,
θIn = ψn, 1 ≤ n ≤ N.
Here θin, Din+1=2, Kin 1=2 are values of the functions θ(z, t), D(θ(z, t)), K(θ(z, t)) at the points (ih, nτ ), ((i + 1/2)h, nτ ), ((i − 1/2)h, nτ ), correspondingly.
Approximate the left boundary condition in the form θ0n+1 − θ0n = 2 (D1n=+21 θ1n+1 − θ0n+1 τ h h − K1n=+21 + Rn+1 − En+1) , 0 ≤ n < N, where Rn+1, En+1 are values of functions R(t) and E(t) at the points t = (n + 1)τ .
Thus, the discrete analog of the direct problem (1)-(2 has a form
( 1 2 ) 2 Φ0n = − τ + h D1n=2 θ0n + h
θmin ≤ θ0n ≤ θmax, Φin = h12 Din 1=2θin 1 + h12 Din+1=2θin+1 − 1 ≤ n ≤ N, { 1
τ
D1n=2θ1n + τ1 θ0n 1 + h2 ( + { θin 1 τ + h1 (Kin 1=2 − Kin+1=2 )} + h12 (Din+1=2 + Din 1=2 )}
θin+ = 0, 1 ≤ i ≤ I − 1,
1 ≤ n ≤ N, −K1n=2 + Rn − En) = 0, ΦIn = θIn − ψn = 0,
1 ≤ n ≤ N, θi0 = φi, 0 ≤ i ≤ I. (3) (4) (5)
The diffusion coefficient and the hydraulic conductivity at the intermediate points appearing in formulas (3) are calculated by the formulas
Din+1=2 =
Dn 1 + Din+11 , i 2
Kin+1=2 =
Kn 1 + Kin+11 , i 2 1 ≤ n ≤ N, 0 ≤ i < I. 4
Introduce a set Q0 = {(z, t) : z = ih, t = lτ, (i, l) ∈ A}, where A = {(i, l) : i = 0, 1, . . . , I, l = 1, . . . , d}, where 0 < d ≤ N , d is some natural number. Denote vector of desirable parameters by u, u ∈ U , where U = {u : u ∈ R3, 0 < a[i] ≤ u[i] ≤ b[i], i = 1, 3}. Define the objective in the form
Due to the form of the objective function, Levenberg–Marquardt algorithm [Levenberg, 1944], [Marquardt, 1963] of numerical optimization can be applied to the solution of the considered optimal control problem. This method is a combination of the Gauss-Newton algorithm with gradient descent method. And in this case, the exact values of the Jacobian of the soil moisture function [θ01, θ11, . . . , θI1, . . . , θ0N , θ1N , . . . , θIN ]T is proposed also to be calculated using FAD method. 4.1
Optimization by Levenberg-Marquardt Algorithm We rewrite the objective function in the form (θjn − θˆjn)2 , where A = {(i, l) : i = 0, 1, . . . , I, l = 1, . . . , d}, d is some natural number, d ≤ N . Denote [θ01, θ11, . . . , θI1, . . . , θ0N , θ1N , . . . , θIN ]T and [θˆ01, θˆ11, . . . , θˆI1, . . . , θˆ0N , θˆ1N , . . . , θˆIN ]T by Θ and Θˆ respectively. According to the Levenberg–Marquardt optimization algorithm at each iteration step k, the displacement vector ∆(uk) is determined from following system of equations: (J (uk)T J (uk) + λdiag(J T (uk)J (uk))∆uk = −J T (uk)(Θ − Θˆ ),
dF (z(u), u)/du = Fu(z(u), u) + GTu (z(u), u)p.
The vector p ∈ Rn from this formula is Lagrange multiplier which is determined as a result of the solution of following linear system of equations :
Fz(z(u), u) + GzT (z(u), u)p = 0n.
The system (10) is linear with respect to p and adjoint to the initial system (8).
Thus, in accordance with the formulas (8)-(10), we obtain the relations for computing gradient of V = θin, i = 0, I, n = 1, N dV (θ(u), u)/du = Vu(θ(u), u) + ΦTu (θ(u), u)p,
V (θ(u), u) + ΦT (θ(u), u)p = 0L, where ΦT = [Φ10, Φ11, . . . , ΦI1, Φ20, Φ21, . . . , ΦI2, . . . , Φ0N , Φ1N , . . . , ΦIN ], p ∈ RL is Lagrange multiplier, u ∈ U ⊂ R3, θT = [θ01, θ11, . . . , θI1, θ02, θ12, . . . , θI2, . . . , θ0N , θ1N , . . . , θIN ], L = (I + 1)N . where J (uk) is the Jacobian of the function Θ(uk): The parameter λ is positive and may be adjusted at the each iteration. The Jacobian of Θ(u) is proposed to be calculated using FAD formulas. 5
Fast automatic differentiation method allows to compute derivatives of complex functions whose variables are related by functional relationships. Briefly describe the essence of this method.
Let for vectors z ∈ Rn and u ∈ Rr continuously differentiable functions F (z, u) and G(z, u) define mappings F : Rn × Rr −→ R1 and G : Rn × Rr −→ Rn. Let z and u satisfy the system of n scalar algebraic equations
G(z, u) = 0n, where 0n is zero n-dimensional vector. Suppose that the matrix GzT (z, u) is not singular. We denote the matrix transposed to the matrix Gz(z, u) by GzT (z, u). Then according to the implicit function theorem, the relations (8) define continuously differentiable function z = z(u). And according to FAD method, the gradient of the function F (z(u), u) is calculated by following formula:
Described approach was applied to finding the numerical solution of the discrete optimal control problem. The problem was solved with following values of input parameters:
L = 100(cm), φ(z) = 0.3,
T = 1(d), z ∈ (0, L), θmin = 0.05(cm3/cm3), ψ(t) = 0.3, t ∈ (0, T ), θmax = 0.5(cm3/cm3), a = [0, 0.005, 0.01]T , b = [300, 0.1, 0.5]T .
The grid with I = 100 and N = 96 was used. The calculations were curried out in three stages. 6.1
At this stage, the direct problem (3)-(4) with the parameters K0true = 100(cm/d), αtrue = 0.01 and mtrue = 0.2 was solved. It is clear from the form of the system (3) and formulas for diffusion coefficient and hydraulic conductivity at intermediate points (4) that the system (3) can be split into N subsystems which of them corresponds to certain time layer. Each subsystem can be solved independently from others subsystems. For each such subsystem, the basic matrix is tridiagonal. Therefore, each subsystem was solved by tridiagonal matrix algorithm. Obtained solution was taken as a prescribed function θˆ(z, t). 6.2
At the second stage the numerical solution of the optimal control problem was searched by steepest descent method. Exact gradient of the objective function (5) was calculated by formulas of FAD. Step value along the chosen direction was determined as a result of procedure of one-dimensional optimization along this direction of the function interpolating the objective function by means of splines constructed on 40 points. We considered the optimal control problems with the objective function (5), where d varied from 1 to 10. For all problems numerical optimization was curried out with various initial approximations. Each initial approximation differed one from another by value of the first component, namely initial approximations for α and m were equal to 0.03 and 0.13 correspondingly and K0init was equal to 102, 105, 110, 120, 150, 180, 200.
Numerical calculations showed that for each initial approximation, the results improve slightly with increasing d from 1 to 10. At the same time, the deviation of the obtained values of the parameters αopt and mopt from the true values αtrue and mtrue depends on the initial approximation. So, with the change of K0init from 102 to 200, this deviation varies from 1.1% to 44.7% for α and – from 0.25% to 8.2% for m.
As to the parameter K0, its value does not practically change during the optimization process and remains very close to the initial value. And, the further the initial approximation of the parameter K0 from its true value, the smaller the difference between the initial value and the obtained value of the parameter K0. This difference does not exceed 1.44·10 5. Presumably, this inability of K0 to be optimized is due to the fact that the corresponding component of the gradient of the objective function differs from other components by 3-4 orders of magnitude.
The values of the parameters K0opt, αopt and mopt obtained for various initial approximations are presented in Figure 1. Under various initial approximations, we mean that the value of K0init changes, while the values of αinit and minit remain unchanged. The graphs, following from left to right in Figure 1, refer to K0opt, αopt and mopt correspondingly. These graphs are designated by solid line. Dashed line shows true values of the parameters. 200 180 160 140 120 100 0.014 0.0135 0.013 0.0125 0.012 0.0115 0.011 0.0105
Thus, these numerical calculations showed that the application of the steepest descent method to solving K0, α, m identification problem does not lead to satisfactory results. As the numerical experiments at the second stage showed, to identify parameters K0, α and m with good accuracy, another (not the steepest descent method) algorithm should be applied. Therefore, in order to solve the discrete optimal control problem, the Levenberg-Marquardt algorithm was applied. The objective function was defined by formula (6), where A = {(j, n) : j = 0, I, n = 1, . . . , d}. The cases of d = 2, 3, 4, 5, 6, 7, 8, 9, 10 were considered. At each iteration in determining the direction of the search, the exact derivatives of θjn, n = 1, N , j = 0, I, were calculated using FAD formulas (11)-(12). The process of the numerical optimization started with initial approximation K0init = 200, αinit = 0.03 and minit = 0.13. The parameter λ changed during the process of numerical optimization. At the beginning it was equal, as a rule, to 10 4 and then it was adjusted at each iteration. The process of the numerical optimization continued until the Chebyshev norm of the gradient of the objective function became less than 1.1·10 18 and value of the objective function became less than 2·10 24. The results of the calculations are presented in following table:
It can be seen from Table 1 that optimal values of the parameters K0, α and m are getting closer to their true values as d increases from 2 to 10. So, while d increases from 2 to 10, the deviation of K0opt from K0true decreases from 0.21% to 0.002%, the deviation of αopt from αtrue decreases from 0.11% to 0.0008%, and the deviation of mopt from mtrue decreases from 0.03% to 0.0002%. The time required to find the solution turned out to be approximately the same for all problems considered.
The value of d defines the set where measured and calculated values of soil moisture are compared. Therefore, the choice of d defines initial data required for determining the parameters. Thus, analyzing results of the numerical calculation, we can estimate how choice of one or another set of initial data will influence on the accuracy of the solution obtained.
Analysis of the results of the numerical calculations leads us to following conclusion.
• The application of the Levenberg-Marquardt method to solving parameters identification problem allows to determine these parameters with good accuracy. So, we can determine the parameter K0 with accuracy up to 0.002%, the parameter α – up to 0.0008% and the parameter mopt – up to 0.0002%. • The gradient method turned out to be ineffective in determining three parameters K0, α and m. In particular, the difficulties in solving this problem are due to the fact that one component of the gradient of the objective function differs from the other components by 3-4 orders of magnitude.
It should be noted the disadvantage of the proposed approach: the Levenberg-Maquardt algorithm used in the process of numerical optimization is one of the local methods. And therefore, there is the question of the possibility of applying in this situation a method of global optimization, for example, of the well-known uneven coating method [Evtushenko & Posypkin, 2013]. This work was supported by Russian Fund of Fundamental Researches, Project 15-07-08952.