ADJOINT-BASED PARAMETER AND STATE ESTIMATION IN 1-D MAGNETOHYDRODYNAMIC (MHD) FLOW SYSTEM

. In this paper, an adjoint-based optimization method is employed to estimate the unknown coeﬃcients and states arising in an one-dimensional (1-D) magnetohydrodynamic (MHD) ﬂow, whose dynamics can be modeled by a coupled partial diﬀerential equations (PDEs) under some suitable as- sumptions. In this model, the coeﬃcients of the Reynolds number and initial conditions as well as state variables are supposed to be unknown and need to be estimated. We ﬁrst employ the Lagrange multiplier method to connect the dynamics of the 1-D MHD system and the cost functional deﬁned as the least square errors between the measurements in the experiment and the numerical simulation values. Then, we use the adjoint-based method to the augmented Lagrangian cost functional to get an adjoint coupled PDEs system, and the exact gradients of the deﬁned cost functional with respect to these unknown parameters and initial states are further derived. The existed gradient-based optimization technique such as sequential quadratic programming (SQP) is employed for minimizing the cost functional in the optimization process. Fi-nally, we illustrate the numerical examples to verify the eﬀectiveness of our adjoint-based estimation approach.


1.
Introduction. Magnetohydrodynamic (MHD) is the study of the interaction of electrically conducting fluids (such as a liquid metal, a plasma) and electromagnetic forces. The description of MHD flows involves both the equations of fluid dynamics, the Navier-Stokes equations, and Maxwell's equations, which are coupled through the Lorentz force and Ohm's law for moving electrical conductors [5]. MHD flows have considerable theoretical and practical importance because of the wide range of applications in science and engineering problems, i.e., nuclear fusion technology with liquid metal [13], magnetic fusion plasmas [8], mixing enhancement for cooling systems [20], microfluidic devices [18,17] and channel fluid flow [14,24].
For MHD flow systems, a great amount of problems concerning in the areas of the flow control and the MHD stability have been investigated with different purposes. For example, In [20], a mixing enhancement by boundary feedback control law is proposed for a two dimensional MHD channel flow. In [24], the stabilization of a linearized MHD channel flow is achieved using the recently developed backsteppingbased boundary feedback control technique and extended to three dimensions for stabilizing the velocity and the pressure as well as the electromagnetic fields in an MHD flow [23]. Other recent developments involving nonlinear model reduction and optimal control have also been investigated [1,6]. These related control and stability problems for the MHD flows are all model-based approaches, in which the model parameters or the initial conditions are supposed to be well predefined and they are essential parameters to these spatial-temporal models. However, once such model parameters (i.e., empirical parameters and initial conditions) are unknown exactly or uncertain, it will usually cause the simulation results large biases and inconsistency between the output of established control system and the real physical one. Thus, it is fundamentally important to identify these various unknown or uncertain parameters in the direct system model so that the mathematical model can best fit the real physical system under consideration. Such problems are commonly referred to as parameter estimation problems or inverse problems in the areas of data assimilation. Under most circumstances, the uncertain parameters in the system models are estimated based on the prior knowledge or available observed data. How to effectively estimate these model parameters from existing data has always been a critical issue.
In this paper, we consider an estimation problem for a simplified MHD Hartman flow with characteristics of incompressible and Newtonian (constant viscosity) in 1-D space. The mathematical model of the 1-D MHD Hartman flow can be successfully built under some suitable assumptions and tightly coupled by the flow velocity and electromagnetic fields [9,22]. However, the coefficients of the Reynolds number and the initial conditions in the model are supposed to be unknown and need to be estimated. Only a limited number of output data at different locations in the system are measureable or available observed. Our purpose is to develop an effective computational method based on optimal control theory to estimate these coefficients and initial conditions inversely by utilizing information from the limited number of observed data in such 1-D MHD flow system. We define a cost functional as the least square errors between the observed data in the experiments and the numerical simulation values so that it can penalize the model-data misfit in time and space, and is also subject to our prior knowledge. We formulate the unknown parameter estimation problem as an approximation optimal parameter selection problem with constraints of the dynamic coupled 1-D MHD flow model, and then solve this typical optimal parameter estimation problem using adjoint-based method. The motivation of performing adjoint-based method is that, in general, it is a well-known variational approach for data assimilation (or inverse problem) and also belongs to the field of optimal control theory. The method aims at optimally adjusting a series of unknown decision parameters such as model parameters or any initial conditions as well as any control inputs based on the given data (observed data). The key role of adjoint method is that it can calculate the gradients of the cost functional efficiently with respect to the to-be-optimized parameters using the chain rule and perturbation techniques, and then an gradient-based optimization algorithm such as gradient descent method can be used to estimate the decreasing of the cost functional. Many parameter estimation problems involving lumped parameter systems (e.g., [2,3,12,10]) and distributed parameter systems (e.g., [16,7,19]) have been successfully solved by the adjoint-based method.
In present work, we first employ the Lagrange multiplier method to connect the dynamics of the 1-D MHD system model and the cost functional defined as the least square errors between the measurements in the experiment and the numerical simulation values. Then, we use the adjoint method based on the variational principle to the augmented Lagrangian cost functional, and get an adjoint coupled PDEs system, which can be solved numerically backward in time. As a result, the exact gradients of the defined cost functional with respect to these decision parameters and initial states are derived via the analytical equations which are related with the costate equations. Based on the gradients of the cost functional, the existed gradient-based optimization techniques such as SQP in the optimization process is used to minimize the cost functional and the optimal results are finally obtained. In this paper, we provide a clear demonstration of the effectiveness of the adjoint method to the estimation problem of 1-D MHD flow model based on the limited observed data at different locations. To the best of our knowledge, there are no other researchers have realized such an estimation problem in this 1-D MHD flow system.
The rest of this paper is organized as follows. In Section 2, we describe the dynamic model of the 1-D MHD flow system and formulate it as an optimal parameter and state estimation problem. In Section 3, we apply the adjoint-based method to the estimation problem and obtain an adjoint coupled PDEs. The exact gradients of the cost functional with respect to these decision parameters and states are also derived in Section 3. Then, in Section 4, we present a numerical scheme for solving the theoretical results obtained in Section 3. In Section 5, the numerical simulation results are presented. Finally, Section 6 concludes the paper by summarizing our results and proposing some further research topics.
2. Statement of the estimation problem.
2.1. System dynamics model. We consider a simplified MHD Hartman flow with characteristics of incompressible and Newtonian (constant viscosity) in 1-D space, and the case in which the fluid flows between two parallel solid plates and the velocity is perpendicular to the magnetic vector. Moreover, the pressure in the channel is supposed to be constant and the unit vectors of the velocity, the magnetic fields external electromotive force are generated right-handed. Then, the mathematical model of describing this 1-D MHD flow system can be derived from the viscous incompressible MHD equations [9,22], which can be formulated as the following dimensionless form: where the spatial variable x and time variable t belong to the set (x, t) ∈ [0, L] × [0, T ]; u(x, t) denotes the flow velocity, B(x, t) is the magnetic induction; f (t) is the pressure difference per unit of the channel length; the function β(t) is a given function of the induction of the external magnetic field and can be considered as the control input for the MHD flow; ν and ν m are two unknown coefficients which are related with the Reynolds number of the flow. We assume the boundary conditions on the solid surfaces are zero velocity of the viscous fluid and zero boundary conditions for the magnetic field corresponding to the continuity of the magnetic field strength. Thus, the boundary conditions for system (1) are given as follows: Furthermore, the initial conditions for the dynamic system (1) are given by However, the initial conditions u in 0 (x), B in 0 (x) are supposed to be unknown to us and need to be well estimated. Now, let us consider the estimation problem of the time and space evolution of states {u, B} when the coefficients ν and ν m are unknown, together with the unknown initial conditions u in 0 (x) and B in 0 (x). The only information available to us is a limited number of measurements or observed values of states {u, B} at different locations in the 1-D MHD flow system. Once the coefficients and initial conditions of system model can be successfully recovered by the limited observation data at different locations, then all the state profiles {u, B} of the 1-D MHD system will be fully rebuilt by the numerical simulation.

Cost functional.
There are usually two approaches for solving the parameter estimation problems. One way is to use the the direct empirical formula with observation data to obtain the best parameter values. However, this method is only suitable for some simple parameter estimation problems. Another popular way is the minimization theory based on optimal control theory by minimizing a defined cost functional. This approach can solve many more complexity parameter estimation problems. In our current work, due to the complexity of the 1-D coupled MHD model, we use the second method to minimize the errors between simulations and some lumped observation data of variables u(x, t) and B(x, t). Thus, we propose the following cost functional to be minimized: where T denotes the optimization time horizon; L denotes the spacial length where the system takes place; N is the number of observation values of u(x, t) and M is the are first guessed values of initial conditions of u(x, t) and B(x, t) respectively; ν F and ν mF are first guessed values of parameters of undeterminate constants ν, ν m , respectively. γ 1 , γ 2 , γ 3 , and γ 4 are weighting factors. The term δ A (x − x j ) is an approximation of Delta-Dirac function described by a Gaussian function with a very small variance ). Now we state our optimal estimation problem formally as follows.
Problem P 0 . Given the 1-D MHD model (1) with the boundary conditions (2) and unknown initial conditions (3), our object is to find the optimal coefficients ν and ν m and initial conditions u in 0 (x), B in 0 (x) based on a limited set of observation data at different locations so as to recover all the state distributions of the 1-D MHD system.
3. Gradient computation via variational analysis. Actually, Problem P 0 is an optimal parameter selection problem with decision parameters ν, ν m and conditions u in 0 (x), B in 0 (x). In principle, such problems can be solved as nonlinear optimization problems using the SQP method or other nonlinear optimization-based methods. However, to do this, we need the gradients of the cost functional with respect to these decision parameters and condtions. Since the cost functional (4) is an implicit function of the decision variables depending on the state trajectories u(x, t) and B(x, t), thus computing the gradient of (4) is a non-trivial task. We now develop a computational method, analogous to the costate method in the optimal control of ordinary differential equations [21,11], for computing these gradients.
Theorem 3.1. The gradients of cost functional (4) with respect to the initial con- and the gradients of cost functional (4) with respect to the parameter ν, ν m are given as where and in which the variables λ 1 (x, t) and λ 2 (x, t) are satisfied the following costate systems: with boundary conditons and termial condtions Proof. We firstly define the following augmented Lagrangian cost functional as where λ 1 (x, t), λ 2 (x, t) are the Lagrangian variables. Then, we take the first variation of the Lagrange function (9) with respect to the system variables u(x, t), B(x, t), initial conditions u in 0 (x), B in 0 (x) and parameters ν, ν m . As a result, the first variation δJ can be given by For the second term, called A 1 of L, we first use the integration by part technique, and can obtain Then, the first variation of A 1 with respect to the system variables u(x, t), B(x, t), u in 0 (x) and B in 0 (x), are developed as follows: For the third term, called A 2 of L, we use the integration by part, and can obtain Then, the first variation of A 2 with respect to the system variable u(x, t), B(x, t), u in 0 (x) and B in 0 (x), is developed as follows: All the terms multiplied by δu in (10), (12) and (14) can be collected together and set to zero so as to cancel all the values of δu. This will give the first order of optimality condition or the adjoint variable λ 1 (x, t) described as Similarly, All the terms multiplied by δB in (10), (12) and (14) can be collected together and set to zero. This will give the adjoint variable λ 2 (x, t) described as follows As a result, the gradients of cost functional with respect to the unknown initial conditions u in 0 (x) and B in 0 (x) can be obtained as Similarly, the gradients of cost functional with respect to the parameters ν and ν m can be obtained as and After reforming the adjoint systems and gradients, all the remain terms in (10), (12) and (14) are gathered together and set to zero to satisfy the optimization conditions. As a result, we can obtain the following initial and boundary conditions for the adjoint systems (15)-(16) Recalling that the boundary conditions u(0, t), u(L, t), B(0, t), B(L, t) in system (1) are fixed beforehand. Thus, we can easily know that their variations with respect to t are all zeros. Thus, the last four terms of (21) are all zeros, and both of the third term and the forth term will be vanished when we make each of λ 1 (0, t), λ 1 (L, t), λ 2 (0, t), λ 2 (L, t) are equaled to zeros. Thus, we have the following boundary conditions for adjoint systems (15)-(16) We complete the proof. 4. Numerical scheme and optimization process.

4.1.
Numerical simulation of the state system and adjoint system. In order to solve the 1-D MHD system dynamics (1) and adjoint equations (7) numerically, we will developed the finite-difference method to discretize these PDEs. This method involves discretizing both the spatial and the temporal domains into a finite number of subintervals, i.e., where M and N are positive integers and h = L/M and τ = T /N . Using the Taylor expansion, the numerical approximations for u(x, t) at each point (x i , t j ) can be expressed as follows: where O(τ ), O(h) and O(h 2 ) denote, respectively, omitted first-and second-order terms such that O(τ ) → 0 as τ → 0 and h −1 O(h 2 ) → 0 as h → 0. Similarly, we can also obtain the numerical approximations for B(x, t) at each point (x i , t j ), and numerical approximates for adjoint states λ 1 (x, t) and λ 2 (x, t). As a result, substituting these approximate formulae (e.g., (24)) into 1-D MHD model (1) and adjoint systems (7) and iterate sequentially, we can obtain all the values of u(x, t), B(x, t), λ 1 (x, t) and λ 2 (x, t). More details about finite-difference methods available in [4].

Numerical integration.
Recall the cost functional (4) and the cost functional's gradients from (6). Clearly, both the cost functional (4) and its gradients (6) involve evaluating double integrals of the form To evaluate these integrals, we partition the space and temporal domains using the same equally-spaced mesh points x 0 , x 1 , . . . , x m and t 0 , t 1 , . . . , t n . These subintervals define step sizes h = L/M and τ = T /N . The integral in (25) can be written as the following iterated integral: Applying the composite Simpson's rule [4] twice, we obtain the following approximation: We also note that both the objective (4) and the costate system (7) also involve evaluating integral of the form such as 1 0 g(x)dx. By the same method of composite Simpson's rule [4], we can also obtain the following approximation: More details on numerical integration algorithms are available in [4]. Moreover, based on the the spatial discretizing, the gradients of cost functional with respect to the initial conditions from (5) at point x j , ∀x j = j∆x where 0 < j ≤ M − 1 with ∆x = L/M , M being the number of points, can be further simplified as follows 4.3. Gradient-based optimization process. By incorporating these gradient formulae in (5) and (6) into a nonlinear programming algorithm such as SQP, we can solve the estimation Problem P 0 numerically. Now we give the Algorithm 1 for solving Problem P 0 . Note that Steps 4 − 7 of Algorithm 1 can be performed automatically by some standard nonlinear optimization solvers such as fmincon function in MATLAB.
Algorithm 1. Gradient-based optimization procedure for solving Problem P 0 .
Step 1: Choose the initial guess ν, ν m and u in 0 (x), B in 0 (x).
Step 2: Solve the 1-D MHD model (1) forward in time corresponding to initial guess values from t = 0 to t = T , and solve the adjoint coupled PDE systems (7) backward in time corresponding to initial guess values from t = T to t = 0 to obtain u(x, t), B(x, t), λ 1 (x, t), λ 2 (x, t).
Step 4: Use the gradients information obtained in Step 3 to perform an optimality test. If ν, ν m and u in 0 (x), B in 0 (x) are optimal, then exit; Otherwise, go to Step 5.
Step 5: Use the gradient information obtained in Step 3 to calculate a searching direction.
Step 6: Perform a line search to determine the optimal step length.
Step 7: Compute the new points ν, ν m and u in 0 (x), B in 0 (x) and return to Step 2.

Numerical simulations.
In this section, we will illustrate a numerical example to verify the effectiveness of our computational algorithm proposed in Sections 3-4. Our numerical simulations are conducted within the MATLAB programming environment running on a desktop computer with the following configuration: Intel Core i5-6500 3.20GHz CPU, 8.00GB RAM, 64-bit Windows 7 Operating System. Our codes implement the gradient-based optimization procedure in Algorithm 1 by combining the optimization tool called fmincon function in MATLAB. The gradients (5)-(6) computed in Section 3 are provided to the fmincon function. In our numerical validation procedures, all the measurements or observed values of system variables u(x, t) and B(x, t) are generated by the exact numerical simulation of the 1-D MHD flow model with respect to given parameters and initial conditions in MATLAB programming environment using numerical simulation technique such as finite difference method. The purpose is that it will sufficiently verify the effectiveness of our algorithm once our estimated values can converge to these given real values. Moreover, considering that the observation data usually has measurement noise in practical applications, we add the Gaussian normal distribution noise to the observation data in our numerical simulation process. These observation data with noise are also used as the input of the optimization process. In present work, we also assume the positions of the measurements are selected manually. For all the following numerical simulations, we set space length L = 1, observation time horizon T = 2, space step x = 0.04, time step t = 0.001, and weighting factors are set γ 1 = γ 2 = γ 3 = γ 4 = 10 −4 . The input function f (t) and β(t) in model (1) are assumed to be known and given by f (t) = 0.01, β(t) = 2t 2 − 2t.
First, we need to validate the correctness of the gradients obtained by the adjoint equations. Using the method described in [15], we define the following function of α where Θ denotes the decision vector in the cost functional, α is a very small scalar but not too close to zero and h is a vector of unit length (i.e., h = ∇J(Θ)||∇J(Θ)|| −1 ).
The correctness of the gradients can be validated once we can obtain a value for φ(α) that is close to 1. Through our numerical simulation in MATLAB, we note that when the values  of α are chosen as 10 −4 ≤ α ≤ 10 −12 , the values of φ(α) are 1.0080 ≤ φ(α) ≤ 1.0168, which are close to 1. This indicates that a unit value of φ(α) can be obtained and thus demonstrates the gradients obtained by the adjoined equations are correct. Now we consider the first scenario. The real values of coefficients ν and νm are given by 0.12 and 0.02, respectively. The real values of initial conditions u(x, 0) and B(x, 0) are also given, as shown in Figure 1 and Figure 2, respectively. The measurements or observed values u meas j (xj, t), B meas k (x k , t) are firstly generated by these given coefficients and initial conditions. We choose M = 20, M is the number of point for the observation data. Based on our optimal computational method proposed in Sections 3-4 and combined the optimization function fmincon, our program takes 3.3789 seconds and converges after 99 iterations. The estimated values of coefficients ν and νm are obtained 0.1218 and 0.0195, and the relative errors between the real values and estimated values are 1.5%, 2.5%, respectively. The estimated initial conditions u in 0 (x), B in 0 (x) are also obtained, as shown in Figure 1 and Figure 2, respectively. The maximum relative errors are 2.61% and 3.79%. In Figure 3, we also give the detailed iteration process of the cost functional (4) in the numerical optimization process. As we can see from the numerical results, all the estimated values converge to the expected real values indeed with some very small bias values at some points. The reasons for the estimation biases are possibly caused by the numerical errors on direct and adjoint system simulation by the finite-difference method or convergence limitation of the optimization function fmincon. Furthermore, at the end of our optimization process, the minimal value of cost functional (4) is 4.4388 × 10 −3 , which is sightly larger than zero, the theoretical value, because of small estimation biases. The numerical results validate the effectiveness of our estimation method and also show that our gradient-based method has strong robustness to find the optimal values. As a result, based on our optimal estimated coefficients ν, νm and initial conditions u in 0 (x), B in 0 (x), all the state information of the MHD flow system can be successfully rebuilt by simulation, as shown in Figure 4 and Figure 5. That means we can recover all the unknown states of the system just only by using the limited observed data at different locations. This point is very important and essential to the estimation problems with unknown dynamic system model, especially in the fields of data assimilation.
Then, we consider another scenario. In this case, the real values of coefficients ν and νm are set as 0.08 and 0.05, respectively. The real values of initial conditions are also given, as shown in Figure 6 and Figure 7, respectively. Similarly, the measurements or observed values u meas j (xj, t), B meas k (x k , t) are also generated by these given coefficients and initial conditions. We also choose M = 20. Based on our optimal computational   method and combined the optimization function fmincon, our program takes 12.8597 seconds and converges after 616 iterations. The estimated values of coefficients ν and νm are obtained 0.0798 and 0.0499, respectively, and the relative errors between the real values and estimated values are 0.25%, 0.02%, respectively. The estimated initial conditions u in 0 (x), B in 0 (x) are also obtained, as shown in Figure 6 and Figure 7, respectively, and the maximum relative errors are 0.48% and 0.99%. Similarly, we give the detailed iteration process of the cost functional (4) in the numerical optimization process, as shown in Figure  8. The numerical results show that all the estimated values also converge to the expected real values indeed with some very small bias values at some points. Furthermore, at the end of our optimization process, the minimal value of cost functional (4) is 9.1620 × 10 −4 , which is sightly larger than zero, the theoretical value. The results also valid the effectiveness of our estimation method. Similarly, based on our optimal estimated coefficients ν, νm and initial conditions u in 0 (x), B in 0 (x), all the state information of the MHD flow system can be successfully recovered by simulation, as shown in Figure 9 and Figure 10, respectively. 6. Conclusions. In this paper, an inverse problem in 1-D MHD Hartman flow system has been studied. We formulated the problem as an optimal parameter estimation problem, in which the unknown parameters and states are successfully estimated under our proposed numerical optimization framework. The adjoint equations coupled with the original system and the gradients of cost functional with respect to these decision variables are derived based on adjoint method. Both direct and adjoint system are solved numerically by means of finite-difference scheme. The numerical examples are illustrated and verified the effectiveness of our proposed computational method in this paper. The optimization process indicates our methodology is very effective to estimate such unknown variables of the dynamic system just only using some limited measurements or observed data. In fact, the estimation scheme proposed in this paper can also be extended to other state and parameter identification problems, i.e., 1-D hydrological system, 2-D transport model in plasmas and switching PDE systems. More general estimation problems in MHD flow control can also be considered further and solved by our approach, and these will be our future research topics.