Structural Identification Problem

The identification problem of the existing structures though the Quasi-Newton and its modification, Trust region algorithms is discussed. For the structural problems, which could be represented by means of the mathematical modelling of the finite element code discussed method is extremely useful. The nonlinear minimization problem of the L norm for the structures with linear elastic behaviour is solved by using of the Optimization Toolbox of Matlab. The direct and inverse procedures for the composition of the desired function to minimize are illustrated for the spatial 3D truss structure as well as for the problem of plane finite elements. The truss identification problem is solved with 2 and 3 unknown parameters in order to compare the computational efforts and for the graphical purposes. The particular commands of the Matlab codes are present in this paper.


Introduction
The problem of definition of the current condition of the existing structure is always one of the most important problems, which arises during the life cycle of the building.An estimation of the material properties is an everyday work of the surveying firms and companies, which deal with the renovation and repair of the buildings.
There exist two approaches for the evaluation of the material properties.The first is a destructive way in which the consistency is interrupted.For instance, taking samples and consecutive testing in the lab or the well-known in-situ pulling test belong to the first category of the surveys.In the other hand the non-destructive techniques are also widely used in such services and moreover for the buildings of cultural heritage.
The method, which is described here can be applied in both cases as for the entire structure as well for the part of the interest.The field of application has a wide range: from the simple truss structures to the planar surfaces, and here we would describe the cases for the 3D truss and CST finite element code.
Depending on the scope of the survey it can be utilized for the evaluation of the condition of the damaged structures or the current load-carrying state.

The identification problem
In the general framework the identification problem can be extremely useful in the cases when the points of interest are actual stress state or the material properties.The basic concept of the identification problem for the 3D truss structure consists in the procedure of application of the static load with the known value and the following estimation of the nodal displacements {u}.Due to the absence of the possibility to test the real structure, we simulate the measurement by the direct procedure of the Finite Element code.In this case we are able to obtain the displacements which are provoked by the certain load condition.
After the construction of the stiffness matrix and definition of the solution which is now and later called as known.We need to repeat the entire procedure again with inserting instead the known quantities the symbolic variables.
As the consequence, we can construct the vector of unknowns {uUnknown}, which in each string contains the expression by means of symbolic variables, as it is written in the Equation (1).
Finally, we can obtain the equality ^`^`^` 0 uUnknown uKnown , which describes the compatibility of the nodal displacements in the certain directions.In order to obtain the values of unknown parameters we need to apply the minimization procedure [10,16].

Least squares method
The basic formulation of the minimization problem can be defined through the least squares method, which finds the parameters that minimize the sum of the squares of the errors between experimental data {uUnknown} and model predictions {uKnown}.In the other words in our problem we are identifying material properties by means of formula (2): Where, k is the number of load cycles, E is an objective function to minimize.It is should be noted that the previous expression is nonlinear and depends on the number of desired parameters to identify.

Advanced formulation of least squares method
In addition to the basic least squares method there exists more advanced formulation, which takes in account the variance-covariance matrix which takes in account the statistical information and denoting by x the desired parameters, so then the objective function E(x) becomes: The challenge of construction [C(x)] matrix arises the problems of definition of statistical data based on the observations or calculated using a model, which for the sake of simplicity are omitted in the further calculations [14].

Unconstrained nonlinear optimization
If to denote by x our desired values then the problem can be formulated as MIN E(x) which means that we are trying to identify a vector x that is a local minimum to a scalar function E(x), without any restrictions to the place of x.
[A] In our case we can consider the unconstrained minimization problem to minimize E(x) by means of trust region.Suppose you are at a point x in k-space and you want to improve, i.e., move to a point with a lower function value.
The basic idea is to approximate E(x) with a simpler function q, which reasonably reflects the behaviour of function f in a trust region (neighbourhood) N around the point x.A trial step s is computed by minimizing (or approximately minimizing) over N.This is the trust-region sub-problem, MIN {q(s), s∈N}.
The current point is updated to be x + s if E(x + s) < E(x); otherwise, the current point remains unchanged and N, the region of trust, is shrunk and the trial step computation is repeated.
The sequence of steps in the unconstrained minimization using trust-region can be shown as: x Formulate the trust-region sub-problem x Solve MIN ^1 ; 2

T T s Hs s g suchthat Ds d '
where g is the gradient of E at the current point x, H is the Hessian matrix (the symmetric matrix of second derivatives), D is a diagonal scaling matrix, Δ is a positive scalar, and . is the 2-norm.This is done in order to determine the trial step s x If E(x + s) < E(x), then x = x + s x Adjust Δ These four steps are repeated until convergence.The trust-region dimension Δ is adjusted according to standard rules.In particular, it is decreased if the trial step is not accepted, which can be stated in the Equation ( 4).

E x s E x t (4)
The basic requirement of this optimisation procedure is the presence of the gradient of the function, which we wish to minimize [2].
[B] Another way to define the minimum of the function is the possibility to use Newton or Quasi-Newton algorithm.They require higher computation effort, but the derivative of function is not needed.
These methods are described by the formulation of Equation ( 5) and build up curvature information at each iteration to formulate a quadratic model problem of this form.
where, H is the Hessian matrix (positive definite symmetric matrix), c is a constant vector, b is a constant, The optimal solution for this problem occurs when the partial derivatives of the Equation ( 6) with respect to x go to zero ) 0 ( f x Hx c (6) And so the optimal solution x* becomes of the form of the Equation ( 7) [15].

Bayesian approach
This method does not provide a single value, for the identified parameter, but probability distribution is of material properties is to be estimated.In Bayesian approach we are trying to identify the distribution of material properties X having an instance x.This vector is characterized by probability density function ( ) The vector of measurements u Known is assumed to be contaminated by uncertainty errors meas e , stemming from the random variable meas 6 .Then the true but unknown response is represented by the Equation ( 8).

meas
Unknown Known e u u (8) and we can denote by meas true U the random variable of true response, so that the Equation ( 8) can be translated by the Equation (9).
Recalling the definition of x -our seeking random value, we can say state that the Equation ( 9) can be converted to the Equation ( 10), using the definition [7].
where, model 6 is a random variable due to the approximations involved in the input parameters p at the modelling stage, u po (x) is the function of material properties and also other p input parameters [3].Bayesian identification involves calculating, for given material properties -x, the likelihood that the true response is the same whether derived from the measurement or from the model.

Code for the direct problem
In the problem solution we have used the code, which solves every kind of truss structures under the hypothesis of linear elastic and first order approach [13].
In order to simulate the measurements of the nodal displacements we need to run the direct procedure, which estimates the values by means of the equilibrium.
The models are drawn by means of a pre-processor software GMesh in which the user inserts the nodes and the rods of the truss structure.The drawn data are contained in a file with a .geoextension that can be read by Matlab like a text file.
Through the script read_input_analysis.m Matlab acquires from the file .geo the number of nodes, number of elements, Coordinates Vector and connectivity vector.Every node and every beam have an identification number.
Young Modulus and area of every beam, the number of restrained nodes and loaded nodes, the degree of restrained of every restrained node, the load on every loaded node are also required data.
Following the directions in Matlab Command Window, these data can be inserted manually and can be saved in a file ".mat" that can be loaded in order to run the analysis.
Once input analysis phase is concluded, the code starts to initialise vectors and matrixes that will be used in the computation.
The main part of the code is the assembly of stiffness matrix, obtained performing a loop over the elements.
It starts with the identification of the number of nodes and their coordinates then computes the length of the element and the transformation matrix.
The transformation matrix is a rotation matrix to transform a vector from local reference system to global reference system.
Let us consider a beam in the 3D description of the rotation matrix, represented by the Equation ( 12).
By C and S are denoted corresponding trigonometric functions Sine and Cosine.Knowing the rotation matrix, we can write the stiffness matrix of the element in the global reference system.After that identifying the degrees of freedom of the element, the code assembles the global stiffness matrix K.
In order to consider the restrained nodes, the code computes the load vector and modifies the stiffness matrix after the assembling.
Defining an α parameter such as a big number and summing it in K matrix terms, it is possible to take into account the infinite rigidity given by restrains.
Finally, it is possible to obtain the solution of the linear system in terms of the vector {u} representing the displacement field (Table 1).%Solution of the linear system u=K\F; Knowing {u} vector it is possible to calculate the reaction forces and the strain in every beam [4].

Code for the inverse identification problem
The identification procedure starts with the initialization of a vector of unknowns in which there are the Youngs modulus written in symbolic form (Table 2).

YoungUn=sym('E',[1 length(Young)]);%Initialization of vector of unknowns
After that, the same procedure of direct problem is performed until the solution of linear system, by doing this it is possible to obtain a vector of displacement that depends by the modulus of elasticity of the elements.
At this point, a vector of displacement as a function of Young modulus and the numerical results given by the direct procedure are available.
In order to find the value of Young modulus it is possible to state the Equation ( 13) should be admitted.

SPbWOSCE-2015
01008-p.5 ^`1 By estimating the norm in L 2 of the vector {u diff } it is possible to obtain the function to minimize [4,17,18].

Code for the direct analysis
The code solves linear elastic analyses for plane problems in plane strain or plane stress formulations.The models are drawn by means of a pre-processor software GMesh in which the user inserts the nodes and the contour of the section.The drawn data are contained in a file with a .geoextension that can be read by Matlab like a text file.The finite element mesh is obtained by using GMesh and is stored in a .mshfile [20].
From a text file prepared by the user, and through the script read_input_analysis.m, Matlab acquires information about the material, such as Young modulus and Poisson coefficient, and obtains boundary conditions [21].
Once vectors and matrix are initialized, the main goal of the code is the assembly of stiffness matrix [19], obtained performing a loop over the elements.
Then, the element stiffness matrix is computed by means of a script function stiff_t3.mwhich takes in account the choice of plane strain or plane stress [4].
The last step before the solving system is to consider the boundary conditions and the nodal forces.Finally, it is possible to obtain the solution of the linear system in terms of the vector u representing the displacement field.Stresses are computed as post-processing [4].

Code for the inverse identification problem
The identification procedure consists in a repetition of the direct procedure with an initialization of an unknown elastic modulus written in symbolic form and inserted into the element stiffness matrix script function, which becomes stiff_t3un.m.
After that, the same procedure of direct problem is performed until the solution of linear system, by doing this it is possible to obtain a vector of displacement that depends by the modulus of elasticity of the elements (Table 3).At this point, a vector of displacement as a function of Young's modulus and the numerical results given by the direct procedure are available.
In order to find the value of Young modulus it is possible to state the Equation ( 13) should be admitted, that in explicit form leads to the Equation (14).
6 Numerical results

Truss 3D
Considering a structure composed by three beams and four nodes at the Figure 1 (a), it is possible to illustrate the way of computations.
The coordinates of nodes are P_1=(333,333,333) P_2=(0,0,0) P_3=(1000,0,0) P_4=(0,1000,0).The value of nodal load applied to the top is V=10000000 N The elastic parameters are E1=2.0e+05E2=3.0e+05E3=2.5e+05All quantities are expressed in the SI units (N, mm).Performing a structural linear elastic analysis of this model with the commercial finite elements software SAP2000, the following result in terms of displacements are obtained (Figure 2).At this point the displacement field is related the following deformed shape is displayed at Figure 1 (a).The results obtained through the Matlab direct code are collected to the Table 6 and the deformed shape is displayed at Figure 3.The vector {UKnown} = {2.4852,1.3301, -6.7586} mm is obtained.Then, using the inverse procedure it is possible to get a function, which depends on three elastic modules E1 E2 E3.This is exactly the function, which is needed for the further computations and estimation of the unknowns in the next step.
In order to use the Quasi-Newton algorithm of the Matlab Optimization Toolbox it is needed to set the initial guess.For the physical problems the range of the output values can be predicted, that is why starting guess [1.5e+05 1.5e+05 1.5e+05] was chosen.The results of the solution of our nonlinear minimization problem by means of Matlab code are listed at the vector ELA = [2e+05 3e+05 2e+05] are obtained.
Finally, it can be easily noticed that vector ELA of the desired parameters perfectly matches to the vector of the inserted elastic moduli during the direct simulation procedure.

Truss 3d with 2 unknown parameters
Using the geometric parameters of the problem discussed before we can solve the problem with two unknowns.
The vector of simulated measurements is presented at the Table 7.For the inverse procedure we set the symbolic vector of unknowns as [E1 , E2 , 2,5*10 5 ] and solve the problem with two unknowns.
Obtaining the function which we need to minimize it is possible to plot the 2D plot depicting the point of the minimum and 3D surface (Figure 5) from which is clearly seen the coordinates of extreme point E1=2, 0*10 5 and E2=3, 0*10 5 .By means of the inverse procedure of Matlab code we can observe in the Figure 9 the value of the function that goes to the minimum point and the number of iteration.
Matlab performs 32 iteration to obtain the values of the modules of elasticity which are slightly different from the exact ones E1=2.0005*10 5 and E2=3.0001*10 5 .
The variation between exact value of the input and the result of minimization can be provoked by the threshold for the iteration procedure.

CST Finite element
Given the high requirement of computational resources needed for inverting the system in symbolic manner, the CST finite element identification problem has been performed on simple structure and with a very reduced number of finite elements.Nevertheless, after several tests performed on different samples, a rather good result has been observed.
Let us consider a concrete block with unknown elastic modulus subjected to a vertical distributed load q=100 N/mm, shown at the Figure 6.
In the real case procedure, it should be supposed to measure directly the nodal displacements under the load condition.Although, assuming a Poisson ratio 0.3 and a square cross section 400mm x 400mm, it is possible to obtain the displacements field from the direct CST procedure by imposing a value to the Young modulus (in this example YOUNG=24500 N/mmq).Therefore, by means of a plane strain analysis the following displacement field is obtained in Table 9.Using the Quasi-Newton algorithm of the Matlab Optimization Toolbox with a starting guess = [20000] and after 28 iterations (Figure 8), the correct value of the desired parameter ELA = [2.4500e+04] is obtained through the minimisation of the function drown at Figure 8.

Conclusions
In the real cases of the truss structures, the parameters of the most components are defined from the very beginning.The solution is the subject for identification is the matter several parameters.For instance having a complex structure as the object for survey we can definitely name the values related to the substructures.At least the groups of elements, which have the same properties can be defined.What is important is to mention that the structures that have uniform or classified distribution of sectional and material parameters require lower computational effort.This can significantly decrease the number of unknowns and as the consequence the complexity of the problem.
During the testing procedures of the Matlab code, we tried to run the structure with great number of rods, where each parameter is different and unknown and we mentioned that it requires the high computational capacity in order to construct the matrix of unknowns.
The surveying procedures can vary depending on the possibility to achieve the important parts of the structure and the used equipment.Sometimes the measurements are easier to perform in the terms of strains, in this case, the slight modification in the procedure are required, but the philosophy of the solution remains the same.

Figure 1 .
Figure 1.a) 3D Truss model in SAP; b) Obtained displacement in SAP.

Figure 2 .
Figure 2. Results of calculation in SAP.

Figure 7 .
Figure 7. a) Deformed configuration CST.b) Displacement diagram in the direction u1 and u2At this point, by means of the inverse code it is possible to obtain a displacement vector as function of the unknown parameter E.

Figure 9 .
Figure 9.The function to minimize.

Table 6 .
Results of Matlab direct procedure.

Table 7 .
Results of Matlab direct procedure.

Table 8 .
Displacement field of CST model.