Identiﬁcation Problem for Nonlinear Gao Beam

: This paper deals with the identiﬁcation of coefﬁcients in the nonlinear beam model which was ﬁrst introduced by D. Y. Gao in 1996. For the identiﬁcation of coefﬁcients, an optimal control approach is used. The unknown coefﬁcients are material parameters of the beam and play the role of the control variables. The existence of at least one solution of the optimal control problem is proved. For the studied problem the ﬁnite element approximation is provided. Finally some illustrative examples are introduced.


Introduction
Beams are important components in many engineering applications. The Euler-Bernoulli beam is the most popular model but it is linear and its validity is limited only to small deflections. For moderately large deflections it is better to use the nonlinear Gao beam model which was firstly introduced in [1] as a static problem and a few years latter as a dynamic model in [2]. This nonlinear beam model respects the original Euler-Bernoulli hypothesis-i.e., the straight lines orthogonal to the midsurface remain straight and orthogonal to it even after the deformation.
The paper is devoted to the identification of coefficients representing material properties for a static Gao beam model where w = w(x) is a deflection of the beam, E stands for Young's modulus, I is the constant area moment of inertia while functions α = 3tb(1 − ν 2 ), µ = (1 + ν)(1 − ν 2 ) and f = (1 − ν 2 )q depend on the Poisson ratio ν, t stands for half-thickness, b is a width of the beam and q is the distributed transverse load. P is the axial force acting at the point x = L. The aim of the identification is to determine piecewise constant parameters E and ν to the beam with a priori known material distribution. The identification problem is defined by minimizing a cost functional J. A similar problem regarding the Euler-Bernoulli beam model was studied in [3], where the extension to higher-order differential equations was used. Inverse problems-i.e., problems to find simultaneously the solution and the coefficient of the Euler-Bernoulli beam equation-was studied in [4,5]. Here, the original problem was transformed into a higher-order well-posed problem following the idea of the method of variational imbedding. An inverse problem for a dynamic Euler-Bernoulli beam by using spectral data was described in [6]. As far as we know there is only one paper [7] dealing with the inverse problem for a dynamic Gao beam using a collage-based method, where data from [8] were considered for a numerical example.
In this paper, we focus on an identification of material parameters for the static Gao beam model. To the authors' knowledge, such a problem has not been studied yet. The material parameters are in fact coefficients of a nonlinear differential equation and their calculation can be based on the solution of the corresponding state problem for given data. For this purpose, an optimal control approach can be used. The similar idea for the identification of coefficients in scalar elliptical differential equations was used in [9], where a steady state groundwater flow problem was considered as a numerical example and a coefficient of hydraulic conductivity was identified.
The paper is organized as follows. In Section 2, the nonlinear Gao beam model is introduced and known results about the existence of a unique solution are reviewed. The problem for the identification of material parameters for Gao's beam is described as an optimal control problem in Section 3. The existence analysis is provided. In particular, we aim to prove the existence of a solution to this problem and show that the solutions of the state problem depend continuously on material parameters. Section 4 deals with a discretization of the studied problem and convergence analysis. Section 5 studies sensitivity analysis. Section 6 presents numerical examples. Concluding remarks are introduced in Section 7.
The paper deals with nonlinear Gao beam model, which is suitable for moderately large deflections. In many engineering applications it can be very useful to identify material parameters of this beam model with respect to given input data. The contribution of the presented paper can be seen in proof of the existence of at least one solution of the studied problem, the convergence analysis of its discretization and finally the sensitivity analysis that provides the gradient for the efficient numerical realization.

Gao Beam Model
This section is focused on a nonlinear beam model which was introduced by D. Y. Gao. This model is suitable for moderately large deflections and uses the Euler-Bernoulli hypothesis. We take into account that the material of the beam is isotropic and the beam has a uniform cross-section of a rectangular shape. Moreover both a transverse and axial loads are considered in this paper. The governing equation for the Gao beam is given by the fourth order nonlinear differential equation where and w is an unknown deflection, E is Young's modulus, ν is the Poisson ratio. The area moment of inertia I = 2 3 t 3 b is constant with 2t as a thickness and b as a width of the beam. The symbol L stands for the length of the beam. The distributed transverse load is denoted by q and P represents the constant axial force acting at the end point x = L. We distinguish two types of axial force cases, the case with an axial force causing compression P > 0 and the case with axial force causing tension P < 0, see Figure 1.
The beam model needs to be completed by one type of the following stable and unstable boundary conditions (see [10]): (B1) Simply supported beam: w(0) = w(L) = 0; w (0) = w (L) = 0; (B2) Fixed beam: w(0) = w (0) = w(L) = w (L) = 0; (B3) Propped cantilever beam: w(0) = w (0) = w(L) = 0; w (L) = 0; (B4) Cantilever beam: w(0) = w (0) = 0; w (L) = E I w (L) − Eα 3 (w (L)) 3 + P µ w (L) = 0. The unstable boundary conditions are separated from the stable ones by a semicolon. In order to give the weak formulation of the problem (2), we define the space V of the admissible displacements that contains the stable boundary conditions. The subsequent analysis is restricted to the simply supported Gao beam with the boundary conditions (B1), see Figure 1. The study for other types of boundary conditions (B2)-(B4) would follow similar steps. The space of admissible displacements is defined as where H 2 ((0, L)) is the Sobolev space-see [11]. The corresponding variational formulation of the Gao beam Equation (2) reads as follows: where Existence of a unique solution of the considered problem (4) has been studied in [12]. In particular it was shown that if (A1) f ∈ L 2 ((0, L)); (A2) E, I, α, µ are positive constants; then the problem (4) has a unique solution. The symbol P E cr stands for Euler's critical force given by , also shown in [13]. For the simply supported beam Euler's critical force P E cr can be computed by the formula also shown in [14]. Therefore, the critical value P has the following form In the recent paper [15] the derivation of the Gao beam model was analyzed and small corrections of the governing equation were presented. This correction consists in a different definition of the constant µ, instead of µ = (1 − ν 2 )(1 + ν) the modified form µ = (1 − ν 2 ) is proposed. This correction in µ respects the fact that the Gao beam is tougher than the Euler--Bernoulli beam. According to this our paper uses a modified Gao beam model in the form Due to this fact we consider the modified critical value instead of P defined by (7).

Parameter Identification
The aim of the identification of parameters is to determine coefficients of a given differential equation using experimentally measured values. This section deals with the identification of the material parameters given by the Young modulus E and Poisson ration ν in the modified Gao beam Equation (8) by using the optimal control approach. According to the definition of α and f in (3), we consider the following equation in (0, L) as a state problem. Here, coefficients E and ν play the role of control variables. We suppose that the interval (0, L) is decomposed into mutually disjoint open intervals K i , called material where 0 < E min < E max < ∞ are given constants and P 0 (K i ) is the set of constant functions on the subinterval K i . Therefore, U ad is the closed, convex subset of couples of piecewise constant functions on the partition of (0, L). The variational formulation of the state problem (10) reads as where using the definitions of α, f . Taking into account the conditions (A1)-(A3), the modified critical value P given by (9) and the definition of U ad , we arrive at the following existence and uniqueness result for (12): let then for any (E, ν) ∈ U ad there exists a unique solution w := w(E, ν) to the problem (12). This can be proven using the results established in [12]. Note that For the subsequent analysis we assume that (A1)'-(A3)' hold. We will study the parameter identification problem as an optimal control problem-see [16,17]-defined as follows where w(E, ν) solves the state problem (12) and J : V −→ R is a cost functional. (13) If (E * , ν * ) is a solution to (13) and w * := w(E * , ν * ) solves (P (E * , ν * )), then the pair ((E * , ν * ), w * ) is called an optimal pair of problem (13).
For further consideration we need to introduce an appropriate continuity concept in the admissible set U ad . The pair (E, ν) ∈ U ad is a piecewise constant vector function over the partition of the interval (0, L), so that it can be identified with a vector (E, ν) = (( . . , r. By convergence in U ad , we mean the convergence of components in the following sense: As the admissible set U ad can be identified with a bounded and closed subset of R 2r , then from the Bolzano-Weierstrass theorem, it follows that U ad is a compact set. In order to prove the existence of a solution to problem (13), we need a continuous dependence of the solution w(E, ν) to the state problem (12) on material parameters.
and w n := w(E n , ν n ) ∈ V be the solution to (P (E n , ν n )). Then Proof. The proof consists in three steps. Firstly, we show that the sequence {w n } is bounded in V. Step We set v := w n and get From (11) it is clear that 1 − ν 2 n > 0, since 0 ≤ ν n ≤ 0.5. Therefore, In order to estimate the term with the axial force P, we apply modified Almansi's inequality-see [18]. Let a function y be such that y (a) = y (b) and b a y( The inequality (17) can be used with y = w n (x), owing to the boundary condition (B1) w n (0) = w n (L) = 0 and holds for any P ≥ 0. Let P ≥ 0 be given. Then by using (11), (16), (18) and Friedrich's inequality: where c 2 := E min I − PL 2 π 2 > 0 using the assumption (A3)'. If P < 0 then (19) trivially holds with c 2 = E min I. For the right hand side in (15) we get where Hölder's inequality and (11) were used. Finally, from (15) , (19) and (20) we see that {w n } is bounded in V. Therefore there exists its subsequence, for simplicity we denote it as {w n } again, such that w n n→∞ w (weakly) in V.
Step 2. Next we show that w solves (12). We prove that using that E n → E in L ∞ ((0, L)) and (21). Similarly (14), the compact impact embedding H 2 ((0, L)) into C 1 ([0, L]), (21) and the formula a The last limit passage (24) can be done in a similar way. Finally, for n → ∞, we have as w solves (12). Since there exists a unique solution to this problem, the whole sequence {w n } tends weakly to w ∈ V.
Step 3. To prove strong convergence, it is sufficient to show that is the norm with respect to which V is complete. Then From this, strong convergence of {w n } to w in V follows. Now we prove the existence at least one solution of the identification problem (13). To this end we suppose that the cost functional J is continuous in V,-i.e., Theorem 2. Let U ad be given by (11) and J satisfy (25). Then the identification problem (13) has a solution.
There exists a large class of functionals J, satisfying (25). For example, the standard least square form where w(E, ν) is the solution to the state problem for given (E, ν) ∈ U ad , z is a target deflection of the beam and · stands for an appropriate norm in V. It is easy to see that the cost functional J is continuous if, for example, z ∈ L 2 ((0, L)), and · is given by the L 2 ((0, L))−norm.
In practice, the function z can be given by discrete measured values which are provided from an experiment-i.e., we have at our disposal point-wise measurements of z at a finite number of points t 1 , . . . , t m , t i ∈ (0, L), ∀i. In this case, the cost functional J is defined by The function y ∈ V is represented by a continuous function, owing to the fact that H 2 ((0, L)) is compactly embedded into C 1 ((0, L)), see [19]. Thus the function J is continuous and well-defined.

Discretization and Convergence Analysis
Since the analytic solution of (13) is hard to find, we need to use an appropriate approximation of the considered problem. We start with a discretization of the state problem (12) by using a finite element approach, see [20]. Let be the nodal points defining the partition of 0, L into intervals I j = x j−1 , x j of length h j = x j − x j−1 , j = 1, . . . , p and h = max{h j } be the norm of D h destined to tend to zero. We shall suppose that the system {D h } is consistent with the (fixed) partition of 0, L into material elements {K i }, meaning that the end-points of K i , i = 1, . . . , r belong to D h for any h > 0. In addition, we suppose that there exists a constant β > 0 which does not depend on h such that h/h i ≤ β, i = 1, . . . , r. With any D h we associate the finite-dimensional space V h defined by where P 3 (I j ) is the space of cubic polynomials defined on the subintervals I j . The discretization of (13) reads as follows: where w h (E, ν) solves the state problem (29) defined by Using classical continuity and compactness arguments we obtain the following existence result. Theorem 3. Let U ad be given by (11). Then (28) has a solution for any h > 0.
Next, we shall study the relation between the optimal control problem (13) and its discretization (28) for h → 0 + . We use notation (E h , ν h ) to point out that (E h , ν h ) ∈ U ad is used in (29) with (E, ν) := (E h , ν h ). Let us recall that U ad itself does not depend on h.
Let v ∈ V be arbitrary. Then there exists a sequence due to the fact that the system of the discrete spaces {V h }, h > 0 is dense in V. Passing to the limit with n → ∞ in the weak formulation and using (30) and (31) we getā The existence of a unique solution of P (E, ν) is guaranteed, therefore (30) holds for the whole sequence {w h }. Strong convergence can be shown similarly as in Theorem 1.
The main result of this section is the following convergence theorem. For simplicity we denote a sequence and its subsequence by the same symbol.
as follows from Theorem 4. Let (Ē,ν) ∈ U ad be arbitrary but fixed. The classical convergence result says that In addition, the definition of (28) yields and from continuity of J in V we get -i.e., ((E * , ν * ), w(E * , ν * )) is an optimal pair of (13). From the proof it is clear that any other accumulation point of {(E * h , ν * h ), w * h } is also an optimal pair of (13).
Next, we derive the algebraic formulation of the discretized problem using a finite element approach. Let us recall that any (E, ν) ∈ U ad can be identified with the vector (E, ν) = ((E 1 , ν 1 ), . . . , (E r , ν r )) ∈ R 2r , where (E i , ν i ) = (E, ν)| K i , i = 1, . . . , r and the admissible set U ad with the compact set U ⊂ R 2r : Since V h is finite dimensional, one can express the solution w h ∈ V h of (29) as . . , N and dim(V h ) = N. Then the discrete state problem (29), (E, ν) ∈ U ad leads to the system of nonlinear algebraic equations where w ∈ R N is the vector whose components are the coefficients α i of the linear combination, (I K 1 (E) − P K 2 (ν)) w corresponds to the bilinear formā in (12), and K 3 ((E, ν), w)w to the nonlinear form π. Vector q(ν) is the force vector corresponding to the linear form L. The matrices K 1 (E), K 2 (ν) and the vector q(ν) are computed using the material constants (E, ν) ∈ U ad : where M ij = supp ϕ i ∩ supp ϕ j and Q ij = supp ϕ i ∩ supp ϕ j . To compute the nonlinear matrix we use Reddy's formula-see [20]: Using the nodal basis system {ϕ i }, the components of vector w are the values of w h and w h at the nodal points x i = ih, i = 1, . . . , p, Let J : U −→ R be defined by J (w(E, ν)) := J(w h (E, ν)), where w h (E, ν) ∈ V h is the solution to (29). Thus, the nonlinear programming problem is given by where w(E, ν) solves (33). Next we assume that J is generated by the least square function: where S ∈ R m×(2p+2) is the matrix representing the restriction mapping of R (2p+2) onto R m . Finally, the vector z = (z 1 , . . . , z m ) is given, where z j = z(t j ), j = 1, . . . , m, where the points t 1 , . . . , t m are such that t j ∈ (0, L), t j ∈ D h ∀j = 1, . . . , m. For numerical realization of (35) the nonlinear conjugate gradient method will be used-see [16], with a suitable step size. For this purpose we need the gradient of J , which can be computed by using the adjoint state technique. For more details, see [9]. The nonlinear Equation (33) will be solved by using Newton's method.

Sensitivity Analysis
The aim of this section is to compute the derivatives of the discretized cost functional with respect to the control variables. The sensitivity analysis will be performed using the algebraic approach. We focus on the following least square function J -see (36): The partial derivatives of J with respect to E i and ν i , i = 1 . . . , r will be computed using the classical chain rule of differentiation: where ∇ w denotes the gradient of J with respect to w, Recall that w(E, ν) solves the nonlinear system R((E, ν), w(E, ν)) = 0, where The partial derivatives w E i and w ν i at (E, ν) will be computed by using the implicit function theorem. From (39) it follows: , ν), w(E, ν)), , ν), w(E, ν) where and Since ∇ w K 3 ((E, ν), w(E, ν))w(E, ν) = 2K 3 ((E, ν), the expression (42) can be rewritten as , ν), w(E, ν)).
To eliminate w E i (E, ν) and w ν i (E, ν) in (40) we use the adjoint state problem defined by: The vector (m(E), m(ν)) is the adjoint state. Multiplying the first equation in (43) by w E i (E, ν) and the second equation by w ν i (E, ν) we get From (38) and (40) we obtain , ν), w(E, ν)), This, together with (41), defines the components of the gradient J with respect to E and ν. For detailed information we refer to [9,21].

Numerical Examples
In this section numerical realization of the studied identification problem is presented. We consider the Gao beam model with the boundary conditions (B1) and with the following input data: half-thickness t = 0.1 m, width b = 0.1 m and length L = 1 m. Further, q = −4 × 10 6 Nm −1 is the constant vertical load applied on the entire beam. We use two the axial forces P = 10 6 N and P = −10 6 N in the numerical example. The beam is split into two intervals, K 1 = 0, 2L 3 and K 2 = 2L 3 , L -i.e., r = 2, such that each interval is characterized by the different material constants (E i , ν i ), i = 1, 2. The vector z representing the measured data is given for positive axial at the points 0, L 6 , L 3 , L 2 , 2L 3 , 5L 6 , L. These vectors z were determined in the following way. For material parameters (E 1 , ν 1 ) = (21 × 10 10 , 0.2), (E 2 , ν 2 ) = (25 × 10 9 , 0.25) and corresponding axial load, we computed the values of deflection of the beam at the given points and to these values we added the white noise from standard normal distribution.
Therefore, the axial forces P = 10 6 N and P = −10 6 N guarantee the existence and uniqueness of the solution to the state problem. For the construction of V h we used the equidistant mesh 0 = x 0 < x 1 < · · · < x 36 = L.
The nonlinear mathematical programming problem has been solved by our implementation of the nonlinear conjugate gradient method-see [16]. The stopping criterion is taken as the combination of the relative error of gradient and the relative error of the cost functional. The initial approximation is (E 0 1 , ν 0 1 ), (E 0 2 , ν 0 2 ) -i.e., on the first interval K 1 is (E 0 1 , ν 0 1 ) = (19 × 10 10 , 0.22) and on the second K 2 is (E 0 2 , ν 0 2 ) = (27 × 10 9 , 0.22). The numerical results are summarized in Table 1 where it denotes the total number of iterations of the nonlinear conjugate gradient method. The deflection for computed material parameters (E * , ν * ) is shown in Figures 2 and 3.  Numerical computations were realized by Matlab.

Conclusions
In the theoretical part of the paper, we analyzed the identification problem for the Gao beam model. Firstly, the Gao beam equation was introduced. Further, the identification Young's modulus and Poisson ratio, was formulated as an optimal control problem. Next, the existence theorem for simply supported beam was presented together with continuous dependence of the solution to the state problem on material parameters. Moreover, the discretization and convergence analysis of the studied identification problem was discussed. The theoretical results were completed by the numerical example.
Future research will be extended to more general cases, namely to an identification problem for the nonlinear Gao beam model supported by an elastic deformable foundation.

Conflicts of Interest:
The authors declare no conflict of interest.