Djalil Chafa

Didier Concordet

A new method for the estimation of variance matrix with prescribed zeros in nonlinear mixed effects models

Keywords: nonlinear mixed effects models, maximum likelihood, expected maximisation algorithm, logitudinal data analysis, repeated measurements, iterated proportional fitting algorithm, Gaussian graphical models, stochastic inverse problems, pharmacokinetic/pharmacodynamics analysis

We propose a new method for the Maximum Likelihood Estimator (MLE) of nonlinear mixed effects models when the variance matrix of Gaussian random effects has a prescribed pattern of zeros (PPZ). The method consists of coupling the recently developed Iterative Conditional Fitting (ICF) algorithm with the Expectation Maximization (EM) algorithm. It provides positive definite estimates for any sample size, and does not rely on any structural assumption concerning the PPZ. It can be easily adapted to many versions of EM.

Introduction

Nonlinear mixed effects models are widely used in population Pharmacology for Pharmacokinetics & Pharmacodynamics (PK/PD) modelling. Such models can be seen as special cases of repeated measurement data, for which the asymptotics concern the number of individuals rather than the number of measures per individual, see for instance [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF] and references therein. In this article, we show that for nonlinear mixed effects models with Gaussian random effects, Expectation Maximization (EM) like algorithms for the computation of the Maximum Likelihood Estimator (MLE) can be coupled with Iterative Conditional Fitting (ICF) like algorithms in order to take into account a prescribed pattern of zeros (PPZ) in the variance matrix of the random effect. The ICF algorithm has been developed very recently [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF] in the context of directly observed Gaussian graphical models. Finding an adequate approach for generic PPZ in the context of nonlinear mixed effects models is a long standing problem. Our approach provides a true solution for the M step of EM in this context, for any PPZ. It is thus far more satisfactory than the standard approaches used in the existing software packages such as NLMIXED (SAS), NONMEM, nlme (S-Plus and GNU-R), or Monolix.

For instance, a traditional model used in population PK/PD is of the form

Y i = F (X i ) + g(X i , θ)ε i , 1 ≤ i ≤ N, (1)
where Y i is the vector of concentrations/effects observed on the i th individual for the drug of interest. Here F is a known function, often nonlinear. The q-vectors X i represent the unobserved individual parameters assumed independent and identically distributed N (m, Σ), the ε i are N (0, I) , unobserved, independent of the X i and independent. The matrix g(X i , θ) is the Cholesky transform of a positive definite matrix that depends on a parameter θ ∈ Θ ⊂ R p . For instance, Figure 1 represents the maximum concentrations of cortisol (Y i ) obtained after giving fixed doses of ACTH to N = 30 horses. Each individual has its own curve described by the parameter X i . The g(X i , θ) matrix defines the variance heterogeneity of concentrations obtained with different doses. This matrix is often assumed to be equal to θdiag(|F (X i |). The reader will later find a study of a model like [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF] for this cortisol data set.

One of the main goals in population PK/PD is to describe the distribution of the X i 's by observation of the Y i 's. This amounts to the estimation of (m, Σ, θ) from the Y i 's. Recall that the X i 's are not observed. A natural approach is to compute the MLE by maximizing

L(m, Σ, θ) = N i=1 1 |g(x i , θ)| φ g -1 (x i , θ)(Y i -F (x i )) φ m,Σ (x i )dx i (2)
where φ m,Σ (.) is the probability density function of the N (m, Σ) distribution and |g(x i , θ)| denotes the determinant of the matrix |g(x i , θ)|. Except for specific models, such as Gaussian linear mixed effects models, maximum likelihood estimators have no closed form. Several methods have been proposed for estimation of the parameter (m, Σ, θ) in these models. The methods suggested by Beal and Sheiner [START_REF] Beal | NONMEM user's guide. Nonlinear Mixed Effects Models for Repeated Measures Data[END_REF] or Lindstrom and Bates [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF] are based on a linearization of the conditional model ( 1) with respect to the vector X i about 0 or about a posterior mode. Pinheiro and Bates [START_REF] Pinheiro | Approximation of the log-likelihood function in the nonlinear mixed effects models[END_REF], Vonesh and Carter [START_REF] Vonesh | Mixed-effects nonlinear regression for unbalanced repeated measures[END_REF], and Wolfinger [START_REF] Wolfinger | Laplace's approximation for nonlinear mixed models[END_REF] proposed Laplacian approximations of the likelihood. Importance sampling approximations [START_REF] Pinheiro | Approximation of the log-likelihood function in the nonlinear mixed effects models[END_REF], Gaussian quadratures [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF], and pseudo-likelihood methods [START_REF] Concordet | A simulated pseudo-maximum likelihood estimator for nonlinear mixed models[END_REF] have also been investigated. The reader will find a detailed analysis of these methods in the book by Davidian and Giltinan [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF]. More recently, stochastic versions of the EM algorithm have been proposed, see for instance [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] and [START_REF] Wang | EM algorithms for nonlinear mixed effects models[END_REF]. These EM like algorithms converge to the MLE under some regularity and identifiability conditions. In many real situations, the kineticist's knowledge of the drug mechanism imposes a specific independence pattern on some components of X i . This means that the variance matrix Σ contains a PPZ. The estimation of (m, Σ, θ) in the presence of a PPZ in Σ is problematic due to the positive definiteness constraint in the optimization. Pinheiro and Bates [START_REF] Pinheiro | Approximation of the log-likelihood function in the nonlinear mixed effects models[END_REF] studied different parameterizations of Σ that ensure the definite positiveness of the estimate. In particular they suggested the usage of a Cholesky like parameterization. Unfortunately, except for the case where Σ is a block diagonal matrix up to coordinates permutation, Cholesky like parametrizations do not preserve the structure of the PPZ and are thus useless. Kuhn and Lavielle proposed estimating Σ in two steps in the implementation of their EM like algorithm. First, Σ is estimated without any constraint, then zeros are plugged according to the PPZ into the estimate provided by the first step. This method is widely used in practice. Unfortunately, by "forcing the zeros" in this way, nothing guarantees that the obtained estimate is still a positive matrix, and even when it is positive definite, it is not the maximum likelihood in general.

For Gaussian linear mixed effects models, the algorithm of Anderson [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF] deals with any linear hypothesis on the variance matrix of the random effect (a PPZ for instance).

Unfortunately, the estimate is not necessarily positive definite, see for instance [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF]. To our knowledge, no method is available for MLE of nonlinear models such as (1) when the variance matrix Σ of the random effect has a PPZ.

The aim of this article is to propose a general method for the estimation of (m, Σ, θ) in the presence of a PPZ in Σ. The method uses the ICF algorithm to perform the Maximization step of the EM algorithm. In other words, we couple EM and ICF in order to compute the MLE (or at least a stationary point of the likelihood) of (m, Σ, θ) when Σ has a PPZ.

The ICF algorithm was developed recently by Chaudhuri et al. in [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF] to estimate a variance matrix with PPZ of observed Gaussian random variables. In contrast, the random effects X i 's in (1) are Gaussian but not observed, and that is why we couple ICF with EM. The ICF converges towards positive definite saddle-points or local maxima of the likelihood function irrespective of the PPZ. The idea behind ICF is not new in the framework of graphical models, and is inspired by the famous Iterative Proportional Fitting (IPF) algorithm. We refer to [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF] for a review. Some alternative algorithms to ICF are available for specific PPZ, such as chain graph models [START_REF] Cox | Decompositions and estimation of a chain of covariances[END_REF] or non-chordal graph models [START_REF] Dahl | Covariance selection for non-chordal graphs via chordal embedding[END_REF]. The ICF algorithm is attractive because it does not rely on a specific structure of the PPZ.

The rest of the article is organized as follows. In section 2, we give some of the properties and drawbacks of the popular "zero forced" estimator, that consists of plugging zeros according to the PPZ into a full variance matrix. In section 3, we recall the main properties of the EM algorithm for models such as [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF]. Section 4 is devoted to the ICF algorithm, and to the coupling of ICF with EM. Section 5 contains the step by step analysis of a model like (1) for the cortisol data set depicted in Figure 1. In the last section, we perform a simulation study that quantifies the benefit of our EM+ICF approach on the model used for the cortisol data set.

The Zero forced estimator

Assume for example that for some specific model we get the following MLE for Σ:

Σ uc =   4 -3 3 -3 4 -3 3 -3 4   ,
without taking into account the PPZ in Σ. We will refer to this as the unconstrained estimation. If the PPZ consists of Σ 13 = Σ 31 = 0, the "zero forced" estimation of Σ is simply given by

Σ zf =   4 -3 0 -3 4 -3 0 -3 4   .
The unconstrained estimate Σ uc is a positive definite matrix but the "zero forced" estimate Σ zf is not. However we know that for a regular model, the unconstrained MLE is consistent. Therefore Σ uc converges componentwise towards the true matrix Σ with PPZ.

Consequently, there exists a random sample size from which the "zero forced" estimator is a positive definite matrix but this sample size is somewhat difficult to obtain. A possible ploy allowing to build a positive definite consistent estimator of Σ could be as follows. Compute the unconstrained estimator and denote it by Σ zf , the corresponding "zero forced" estimator. Nothing guarantees that its lower eigenvalue λ min is positive but since Σ zf is a consistent estimator of Σ, the quantity (λ min ) -max{-λ min , 0} is a random sequence of positive numbers that converges almost-surely to zero. Now, consider some auxiliary sequence of positive real numbers (u N ) that goes to zero with the sample size N (e.g. U N = 1/N 2 ), then, for any sample size N, the matrix

Σ zf + ((λ min ) -+ u N ) I
is a positive definite consistent estimator of Σ, and features the same PPZ. Its main drawback is that its diagonal terms are biased and that the choice of the (u N ) sequence is arbitrary. A better way to proceed is to directly consider the MLE of Σ with PPZ, which is precisely our aim in the next sections.

The EM algorithm

The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is a popular method to estimate parameters of a model with nonobserved or incomplete data. Let us briefly recall how its general form works as introduced by Dempster et al. The EM algorithm consists of iterations of an Expectation and a Maximization step. At the k th iteration, the E step computes the conditional expectation of the log-likelihood of the complete data (Y, X) with respect to the distribution of the missing, or non-observed, data X given the observed data Y at the current estimated parameter value ψ (k) :

Q ψ, ψ (k) = E log P (Y, X)|Y, ψ (k) .
The M step finds ψ (k+1) so that for all ψ in the parameter space Ψ

ψ (k+1) = arg sup ψ∈Ψ Q ψ, ψ (k) .
These two-step iterations are repeated until convergence. The essential property of the EM algorithm is that the likelihood increases monotonically along the iterations. Under some identifiability and regularity conditions, this algorithm converges to a stationary point of the likelihood, see for instance [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF].

Let us now describe more precisely this algorithm for model [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF]. We need first to define the parameter space on which the M step is to be performed. In this model, the parameter to be estimated is ψ = (m, Σ, θ). The variance matrix Σ lives in a subset of the set S + q of q × q symmetric positive definite matrices. More precisely, let Π be the set of subsets of {(i, j); 1 ≤ i < j ≤ q}. For any π ∈ Π, the set

S + q (π) {A ∈ S + q ; ∀(i, j) ∈ π, A ij = 0}
is formed by the symmetric positive definite matrices that have zeros located in π. The PPZ in Σ is represented by an element π of Π. We thus assume that for some π ∈ Π, ψ = (m, Σ, θ) ∈ Ψ M × S + q (π) × Θ where M and Θ are open subsets of R q and R p respectively.

At the k th iteration the Expectation step consists of the computation of

Q (m, Σ, θ|m k , Σ k , θ k ) = i E log 1 |g(X i , θ)| φ g -1 (X i , θ)(Y i -F (X i )) φ m,Σ (X i )|Y i , Σ k , m k , θ k where E (f (X, Y )|Y, Σ, m, θ) f (x, Y ) φ (g -1 (x, θ)(Y -F (x))) φ m,Σ (x) |g(x, θ)| φ (g -1 (u, θ)(Y -F (u))) /|g(u, θ)|φ m,Σ (u)du dx.
The maximization step computes

(m k+1 , Σ k+1 , θ k+1 ) = arg sup M ×S + q (π)×Θ Q (m, Σ, θ|m k , Σ k , θ k ) .
For model [START_REF] Anderson | Asymptotically efficient estimation of covariance matrices with linear structure[END_REF], the integral that appears in the E step can be split into two parts . The E step reduces to calculate

Q (m, Σ, θ|m k , Σ k , θ k ) = Q 1 (m, Σ|m k , Σ k , θ k ) + Q 2 (θ|m k , Σ k , θ k )
where

Q 1 (m, Σ|Y i , m k , Σ k , θ k ) i E (log φ m,Σ (X i )|Y i , m k , Σ k ) = - 1 2 i E (X i -m) ′ Σ -1 (X i -m) |Y i , m k , Σ k - N 2 log |Σ| = - N 2 tr 1 N i E (X i -m) (X i -m) ′ |Y i , m k , Σ k Σ -1 - N 2 log |Σ|, and
Q 2 (θ|m k , Σ k , θ k ) i E log φ g -1 (X i , θ)(Y i -F (X i )) -log(|g(X i , θ)|)|Y i , m k , Σ k , θ k .
It follows that the M step can also be decomposed into two parts sup

M ×S + q (π)×Θ Q (m, Σ, θ|m k , Σ k , θ k ) = sup M ×S + q (π) Q 1 (m, Σ|m k , Σ k , θ k ) + sup Θ Q 2 (θ|m k , Σ k , θ k ) .
Note that in the M step the maximization with respect to (m, Σ) is separated from that of θ. The function Q 2 depends only on θ via g.

Remark 1. In most applications, h is the probability density function of a standard Gaussian distribution, and θ is a variance matrix that can possibly contain a PPZ. In that case, its maximization can be performed using the ICF method, as for Σ, as described hereafter.

Maximization of Q 1 leads to

m k+1 = 1 N i E (X i |Y i , Σ k , m k , θ k ) (3)
and

Σ k+1 = arg inf Σ∈S + q (π) tr XΣ -1 + log |Σ| (4) where X = 1 N i E (X i -m k+1 ) (X i -m k+1 ) ′ |Y i , Σ k , m k , θ k . (5)
The matrix X is an empirical conditional variance matrix. The random vectors X i are not observed, however, the matrix X can be evaluated at each iteration of the algorithm. When π = ∅, that is when Σ has no PPZ, Σ k+1 reduces to X. When π = ∅ the maximum of Q 1 (Σ) must be sought in S + q (π). The next section deals with this problem.

Estimating the variance matrix with ICF

We have seen in the previous section that the M step of the EM algorithm involves a maximization problem such as

Σ k+1 = arg inf Σ∈S + q (π) tr XΣ -1 + log |Σ|. (6)
The difficulty here is that the optimization is not performed on the entire cone of symmetric definite matrices but only on a sub-cone that contains the matrices with PPZ. It is clear that a standard gradient-like algorithm does not fit these constraints. The usual method to get rid of the definite-positiveness constraint is to use a Cholesky like decomposition. Unfortunately, these decompositions do not preserve the PPZ when Σ is not a block diagonal matrix up to a permutation of the coordinates. The algorithm described hereafter allows one to move within S + q (π) whatever π may be. First, note that in the absence of PPZ in Σ, i.e. when π = ∅, the solution of ( 6) is Σ k+1 = X. We assume from now on that π = ∅. The case q = 2 is trivial since the only possible zero is Σ 12 = 0. In this case, the "zero forced" estimator is always positive definite and is the solution of [START_REF] Cox | Decompositions and estimation of a chain of covariances[END_REF]. We assume in the sequel that q > 2 and that π = ∅.

The method we propose is the core of the ICF algorithm presented in [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF]. Even if Chaudhuri & al. did not express it at such, it is mainly based on the specific properties of the Schur complement of a matrix. Let us recall the following classical result, which can be found in [START_REF] Horn | Matrix analysis[END_REF] or [START_REF] Zhang | The Schur Complement and its Applications, volume 4 of Numerical Methods and Algorithms[END_REF] for instance.

Theorem 1 (Schur). Let l be an integer in {1, . . . , q}. Consider two vectors U and V that respectively belong to R q-l and R l , the q-vector Z ′ = (U ′ , V ′ ) and a matrix Σ ∈ S + q that admits the following block decomposition

Σ = A B B ′ C
where A ∈ S + q-l , B is a (ql) × l matrix and C ∈ S + l . The Schur complement of A is the matrix S C -A ′-1 B. It belongs to S + l and moreover i)

det(Σ) = det(A) det(S), ii) Z ′ Σ -1 Z = U ′ A -1 U + (V -B ′ A -1 U) ′ S -1 (V -B ′ A -1 U) .
The Schur complement appears naturally when the random vector Z described in the previous Theorem is distributed according to N (0, Σ). More precisely we have:

L (V |U) = N (B ′ A -1 U, S).
Note that properties i) and ii) remain true after permutation of the rows of Σ. In the particular case where l = 1, we can easily derive the following property.

Corollary 1. We keep the same notations as in the previous proposition. For any j ∈ {1 . . . q}, let A = Σ -j,-j be the submatrix of Σ obtained by removing its j th row and column , B = Σ -j,j the j th column vector of Σ in which the j th row has been removed, C = Σ j,j . The column vector U and the positive real number V are respectively obtained by removing the j th row of Z and as Z j . Then, using these notations, the Schur complement of A = Σ -j,-j is the real positive number S given by the previous proposition and properties i) and ii) hold.

We are now able to solve the optimization problem (4) by running iteratively the decomposition of the Corollary over the columns of Σ. Set T i = X im k+1 and note that from [START_REF] Cox | Decompositions and estimation of a chain of covariances[END_REF] the function that has to be minimized can be rewritten as

K(Σ, T 1 , . . . , T N ) = 1 N i E T ′ i Σ -1 T i |Y i , m k , Σ k + log |Σ|. (7)
For the j th column of Σ we set A = Σ -j,-j , B = Σ -j,j , C = Σ j,j and S = C -B ′ A -1 B so that we can now write [START_REF] Dahl | Covariance selection for non-chordal graphs via chordal embedding[END_REF] as

K(Σ, T 1 , . . . , T N ) = K(A, U 1 , . . . , U N ) + K(S, V 1 -B ′ A -1 U 1 , . . . , V N -B ′ A -1 U N )
where V i is the j th component of T i and U i is obtained by removing the j th component of T i . Therefore, if A is fixed, the partial optimization of K(Σ, T 1 , . . . , T N ) with respect to (B, S) can be reduced to the global optimization of

K(S, V 1 -B ′ A -1 U 1 , . . . , V N -B ′ A -1 U N ) = 1 NS i E V i -B ′ A -1 U i 2 |Y i , m k , Σ k + log(S)
which is a standard least-squares problem. The optimization with respect to B and S leads to

B opt = i E (A -1 U i )(A -1 U i ) ′ |Y i , m k , Σ k -1 i E V i A -1 U i |Y i , m k , Σ k (8) = A i E (U i U ′ i |Y i , m k , Σ k ) -1 i E (V i U i |Y i , m k , Σ k ) . (9)
and

S opt = 1 N i E V i -B ′ opt A -1 U i 2 |Y i , m k , Σ k . ( 10
)
The vector B = Σ -j,j may contain some PPZ. These components are not optimized and are thus left at zero. This only decreases the dimension of the optimization problem. We deduce that after this step on the j th column of Σ, C = Σ j,j and B = Σ -j,j must be respectively updated with Σ new j,j = S opt + B ′ opt (Σ -j,-j ) -1 B opt and Σ new -j,j = B opt .

The striking property of this step is that it allows us to move within S + q without affecting the prescribed null components of Σ: A = Σ -j,-j and the null components of B = Σ -j,j are left unchanged. Since Σ is assumed positive definite C = Σ j,j cannot be zero.

As already mentioned, iterations of these steps converge to a local maximum of K(Σ, T 1 , . . . , T N ), see for instance [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF].

The cortisol data set

In this section, we use a practical example to illustrate the implementation of an EM algorithm coupled with the ICF algorithm. In order to explore the endocrine function of horse, a sample of horses (N = 30) was given eight doses of ACTH by intravenous route. The ACTH stimulates the adrenal gland that produces cortisol. The concentration profiles of cortisol in plasma were summarized by the maximal concentration reached after the ACTH administration (see Figure 1).

The seven doses of ACTH given to each animal were (in mg/kg) 0.005, 0.01, 0.1, 0.5, 1, 2, 10.

The production of cortisol is modelled as

Y ij = X 1i + X 2i d X 3i j X X 3i 4i + d X 3i j (1 + σε ij ) , 1 ≤ j ≤ 7, 1 ≤ i ≤ 30, ( 11
)
where Y ij is the maximal cortisol concentration observed in the i th horse after administration of a dose d j of ACTH, X ′ i = (X i1 , . . . , X i4 ) is a random vector that contains the individual parameters for the i th animal. We assume that the random vectors X i are independent and identically distributed N (m, Σ) and that the residual terms ε ′ i = (ε i1 , . . . , ε i7 ) are independent and identically distributed N (0, I 7 ). Moreover, the X i 's and ε i 's are assumed to be mutually independent. In this example, p = 4, θ = σ 2 and

g(X, θ) = σdiag X 1 + X 2 d X 3 j X X 3 4i + d X 3 j j=1...7
.

According to the kineticist, the correlations between X i1 and X i4 and X i3 and X i4 should be zero and thus Σ has the following structure

Σ =     . . . 0 . . . . . . . 0 0 . 0 .     .
In some problems, no a priori information is available for the possible zero correlation between the components of X i . The method of multiple testing of correlation, as described in Drton and Perlman [START_REF] Drton | Model selection for gaussian concentration graphs[END_REF], may be used in such cases to reveal the structure of Σ.

The estimation of model parameters requires evaluation of the conditional expectations of functions such as

E (f (X i , Y i )|Y i , Σ k , m k , σ 2 k )
. A standard approach is to use a stochastic version of EM that consists of the simulation of a Markov Chain (X

(l) i ) l with P (.|Y i , Σ k , m k , σ 2
k ) as unique stationary distribution by using a Metropolis-Hastings algorithm and the approximation of the conditional expectation by

E f (X i , Y i )|Y i , Σ k , m k , σ 2 k ≈ 1 L l l=1 f (X (l) i , Y i ).
For the analysis of the cortisol data we chose L = 500. We simulated the Markov chain with the Metropolis-Hastings algorithm with N (m k , Σ k ) as the proposal distribution. In this case, the acceptance probability of the Metropolis-Hastings algorithm reduces to

min φ (g -1 (X, σ 2 k )(Y -F (X))) |g(x (l) , σ 2 k )| |g(X, σ 2 k )|φ (g -1 (x (l) , σ 2 k )(Y -F (x (l) )))
,

1
which only depends on the conditional distribution of the observation.

The algorithm to estimate the model parameters can be summarized in the following scheme:

1. Start from some initial guess Σ 0 , m 0 , σ 2 0 and set k = 0 ; 2. Compute m k+1 from (3) and

σ 2 k+1 = 1 N i E (Y i -F (X i )) ′ C -1 (X i )(Y i -F (X i ))|Y i , m k , Σ k , σ 2 k , where C -1 (X i ) is a diagonal matrix whose j th term is (1/(F (X i )) j ) 2 .

Set Σ

(0) k = Σ k and and set l = 0 4. for j:=1 to q = 4 increment l compute Σ (l+1) j,j

= S opt + B ′ opt Σ (l) -j,-j -1
B opt and Σ

-j,j = B opt where S opt and B opt are respectively defined by ( 10) and (8) ;

5. if Σ (l-4) k and Σ (l) k are not close enough go to step 4. Otherwise set Σ k+1 = Σ (l) k ; 6. Stop if (Σ k , m k , σ 2
k ) and Σ k+1 , m k+1 , σ 2 k+1 are close enough. Otherwise increment k and go to step 2. For standard EM algorithms, i.e. when no constraint is imposed on the structure of Σ, steps 3), 4) and 5) of the previous algorithm should be replaced by the update of Σ k according to equation [START_REF] Concordet | A simulated pseudo-maximum likelihood estimator for nonlinear mixed models[END_REF].

It is well known that standard EM algorithms go quickly to a stationary point of the likelihood during the first iterations and then take time to converge. Since for stochastic EM algorithms the criterion being optimized changes randomly at each iteration, it is somewhat difficult to achieve and check convergence even when the length L of the simulated Markov Chain is large. Improvements have been proposed to overcome these problems. In particular Kuhn and Lavielle [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] suggested updates of the following form :

   Σ new = (1 -γ k )Σ k + γ k Σ k+1 m new = (1 -γ k )m k + γ k m k+1 σ 2 new = (1 -γ k ) k σ 2 k + γ k σ 2 k+1 ,
where (γ k ) k is a decreasing sequence of positive numbers, Σ k+1 , m k+1 , σ 2 k+1 is defined as in the previous algorithm and (Σ new , m new , σ 2 new ) is the update of (Σ k , m k , σ 2 k ). Note that this update scheme forces the algorithm to converge and preserves the PPZ as well as the positive definiteness of Σ.

The sequence (γ k ) k should satisfy k γ k = +∞ and k γ 2 k < +∞. These two conditions are fulfilled when γ k = a/k b with a > 0 and b ∈ (0, 1). Choosing γ k = a/k speeds-up convergence of the algorithm but the choice of a has to be made sample by sample. Choosing the same a for all samples can lead to poor estimations. As practical advice, we suggest choosing γ k = 1/k 0.8 . The algorithm takes more time to converge but a fine tuning of a is unnecessary.

A well known drawback of the EM algorithm is that it does not produce standard errors as a by-product. We implemented the method proposed by Jamshidian and Jennrich [START_REF] Jamshidian | Standard errors for em estimation[END_REF]. This method relies on numerical derivation and seems well-suited to the method we propose. Even if standard errors are helpful for comparing the results obtained with these two models, the Fisher information matrix gives pertinent quantitative information only when the sample size is large enough. However, N = 30 is probably not a large sample size. In the next section we use simulations to weight the performance of the estimation proposed for the cortisol data.

The estimation of the model parameters for the cortisol data requires some initial estimates to be provided. Thanks to the model parametrization, we can directly read reasonable values for m 0 on Figure 1. Since the four components of X respectively represent the basal value of cortisol, the maximal increase, the "slope" of the sigmoid and the ACTH dose for which half the maximal increase is obtained we roughly get m 0 = (50, 70, 1, 0.1). We initialize Σ with the following diagonal matrix: Σ 0 = Diag (0.01m 2 0 ). Finally, for this heteroscedastic model, σ can be interpreted as the coefficient of variation of the cortisol for a given dose. We set it at 20% that is σ 2 0 = 0.2 2 . We estimated Σ with EM alone (no constraint was imposed) and with EM+ICF that preserves the PPZ. In this example, the algorithm seems to converge in less than 400 iterations. We implemented this algorithm in C++ with a matrix library. Estimates of the parameters obtained with EM alone were: [START_REF] Ye | Modelling of covariance structures in generalised estimating equations for longitudinal data[END_REF]. We can see that these likelihoods are about the same and a likelihood ratio test would not reject the PPZ proposed by the kineticist. The residual variance estimates are also very close. Surprisingly, there are quite large differences between the estimates of the third component of m and the non null components of Σ.

Σ =     21.
Remark 2 (Modelling). The general problem of mean and variance modelling for longitudinal data is delicate, and several choices are possible, see for instance [START_REF] Davidian | Nonlinear Models for Repeated Measurement Data[END_REF], [START_REF] Wolfinger | Heterogeneous variance-covariance structures for repeated measures[END_REF], [START_REF] Pourahmadi | Joint mean-covariance models with applications to longitudinal data: unconstrained parameterisation[END_REF][START_REF] Pourahmadi | Maximum likelihood estimation of generalised linear models for multivariate normal covariance matrix[END_REF], [START_REF] Zimmerman | Parametric modelling of growth curve data: an overview[END_REF], and [START_REF] Cepeda | Bayesian modeling of joint regressions for the mean and covariance matrix[END_REF]. Our model [START_REF] Horn | Matrix analysis[END_REF] belongs to a standard family of models in PK/PD and was chosen with the kineticist. This relatively simple model is heteroscedastic with a constant coefficient of variation. An examination of the "individual" residuals shows that they are centered, which is quite satisfactory.

Simulations

The aim of this section is to quantify the potential benefit of directly estimating a variance matrix with PPZ. We simulated 100 data sets using model [START_REF] Horn | Matrix analysis[END_REF] with parameters close to the estimate found in the cortisol data analysis : N = 30,m ′ = (50, 70, 1.5, 0.08),

σ 2 = 0.015 and Σ =     20 -4.5 2.5 -0.3 -0.1 0.05 0 -2 × 10 -3 0 10 -5     .
Both EM and EM+ICF estimates were calculated. Results are given in Table 1.

As expected, the standard errors given in the example are smaller than those of Table 1. For such sample sizes, which are often encountered in practice, asymptotic statistics should be interpreted with care.

The mean parameter m seems to be well estimated. At least on these simulations, the Σ structure influences the estimation of m. However, we notice that the estimates obtained with EM+ICF have a smaller standard error and mean quadratic error (M.Q.E.) than those obtained without any constraint. On the whole EM+ICF also gives estimates with lower bias. This suggests that the mean and variance estimations are heavily dependent. This sheds light on approaches, such as the 'zero forced" method, that rely on estimating the full variance matrix first and modify it by forcing the PPZ: since all the non zero entries are estimated with the assumption that the variance matrix does not have prescribed zeros, they could be poorly estimated. This is consistent with the results obtained by Ye and Pan [START_REF] Ye | Modelling of covariance structures in generalised estimating equations for longitudinal data[END_REF] who concluded, in a different context, that misspecification of the working variance structure could lead to a large loss in efficiency of the estimators of the mean parameters.

Likelihood ratio tests were performed to test

H 0 : Σ ∈ S + 4 (π) H 1 : Σ ∈ S + 4 for π = {(1, 4
); [START_REF] Cepeda | Bayesian modeling of joint regressions for the mean and covariance matrix[END_REF][START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF]}. Note that whatever the value of π, H 0 is not on the boundary of S + 4 . Consequently, the likelihood ratio statistics follow asymptotically a Chi-square distribution under H 0 . Since the data have been simulated under H 0 , the P-values distribution should be close to a uniform law on (0, 1) at least for large N. The Q-Qplot of the P-Values is represented in Figure 2. This figure shows that the P-Values are not distributed according to a uniform distribution and thus the distribution of the likelihood ratio statistics is not close to a χ 2 distribution. Consequently, N = 30 is probably not large enough to trust asymptotic statistics.

Conclusion

We have proposed a method for the estimation of the variance matrix with PPZ in nonlinear mixed effects models. This method, which consists of coupling an ICF like algorithm with an EM like algorithm gives more efficient estimates than standard EM that ignore the PPZ. For the sake of simplicity, we have only presented the estimation algorithm for independent and identically distributed observations. Extension to different numbers of observations per individual is straightforward. We also restricted our study to models with Gaussian ε. More general models in which the distribution of ε is not Gaussian and depends on a parameter θ 2 can also be considered. This simply requires the Metropolis-Hastings chain to be chosen accordingly. We deliberately chose to show in section 2 a columnwise ICF implementation that can be extended using theorem 1 to blocks of Σ.

Of course, our approach can be adapted without much effort to many versions of EM and many alternatives to ICF. For pedagogical reasons, we presented our EM+ICF coupling on a low dimensional example. The method of course is particularly suited to large variance matrices with a high percentage of prescribed zero entries.

Figure 1 :

 1 Figure 1: Maximum cortisol concentrations observed after IV administrations of ACTH in 30 horses.

Figure 2 :

 2 Figure 2: The Q-Q plot of the P-values of the likelihood ratio test versus the uniform distribution on (0, 1).

Table 1 :

 1 Empirical mean, standard error and square root of mean-quadratic-error of the estimates (M.Q.E.) obtained with EM and EM+ICF. The Mean Quadratic-Error is defined as bias 2 + Variance.

		True value Mean	EM S.E.	√	M.Q.E. Mean	EM+ICF S.E.	√	M.Q.E.
	Σ 11	20	17.3	10.74		11.07	18.88	9.09	9.16
	Σ 12	-4.5	-3.64	3.16		3.28	-4.1	2.39	2.42
	Σ 22	2.5	2.27	1.48		1.5	2.74	1.25	1.28
	Σ 13 × 10 2 Σ 23 × 10 2 Σ 33 × 10 3 Σ 14 × 10 4 Σ 24 × 10 4 Σ 34 × 10 5 Σ 44 × 10 6 m 1	-30 -10 50 0 -20 0 10 50	-8.04 -7.57 74.36 -7 49.78 306.23 86.53 19.95 95.6 91.62 -43.69 63.29 17.84 14.15 51.18 1.38	89.27 20.1 98.66 91.89 314.08 76.9 16.18 1.82	-12.27 -10.66 73.67 0 104.01 0 18.57 48.96	39.85 14.65 65.71 0 204.02 0 13.68 1.23	43.62 14.66 69.85 0 238.76 0 16.15 1.61
	m 2	70	70.73	2.60		2.70	69.52	2.44	2.48
	m 3	1.5	1.54	0.22		0.22	1.49	0.21	0.22
	m 4 × 10 3 σ 2 × 10 3 log like.	80 15	91.44 15.91 -749.32 11.39 4.08 1.77		12.15 1.99	90.06 14.07 -751.54	3.87 1.57 11.68	10.78 1.82

Acknowledgements. We thank Dr. Alain Bousquet-Mélou who provided the cortisol data set and for valuable advice regarding the model and the specific structure of Σ. We also thank Dr. Mathias Drton for interesting discussions on the ICF algorithm and for providing the last version of [4] before publication. The final form of this article has greatly benefited from the questions and suggestions of two anonymous reviewers.


Zhibin Yang
email: yangzhibin168@163.com

Zhikai Qiu

Yong Zhou
email: zhouyong@nuaa.edu.cn

Zhiqiu Huang
email: zqhuang@nuaa.edu.cn

Jean-Paul Bodeveix
email: bodeveix@irit.fr

Mamoun Filali
email: filali@irit.fr

C2AADL Reverse: A Model-Driven Reverse Engineering Approach for Development and Verification of Safety-critical Software

Keywords: Safety-Critical Systems, Model-Driven Development, Model-Driven Reverse Engineering, AADL, Compositional Verification

The safety-critical system communities have been struggling to manage and maintain their legacy software systems because upgrading such systems has been a complex challenge. To overcome or reduce this problem, reverse engineering has been increasingly used in safety-critical systems. This paper proposes C2AADL_Reverse, a model-driven reverse engineering approach for safety-critical software development and verification. C2AADL_Reverse takes multi-task C source code as input, and generates AADL (Architecture Analysis and Design Language) model of the legacy software systems. Compared with the existing works, this paper considers more reversed construction including AADL component structure, behavior, and multi-threaded run-time information. Moreover, two types of activities are proposed to ensure the correctness of C2AADL_Reverse. First, it is necessary to validate the reverse engineering process. Second, the generated AADL models should conform to desired critical properties. We propose the verification of the reverse-engineered AADL model by using UPPAAL to establish component-level properties and the Assume Guarantee REasoning Environment (AGREE) to perform compositional

Introduction

Safety-critical systems (SCS) are the systems whose failure could result in loss of life, substantial economic loss, or damage to the environment [START_REF] Leveson | Engineering a safer world: Systems thinking applied to safety[END_REF]. There are many well-known examples in different domains such as aircraft flight control, space missions, and nuclear systems. The SCS communities have been struggling to manage and maintain their legacy software systems because upgrading such systems has been a complex challenge. As surveyed by FAA (Federal Aviation Administration), reverse engineering (RE) has been increasingly used in many industries, including aircraft applications [START_REF] George Romanski | Reverse engineering for software and digital systems[END_REF].

In contrast with forward engineering, reverse engineering can be defined as the process of examining an already implemented software system to create a higher abstraction level representation in a different form. Reverse engineers typically start with a low-level representation of a system (such as source code, or execution traces), and try to build more abstract representations from these (such as architectural models, or use cases, respectively) [START_REF] Van Deursen | Software reverse engineering[END_REF]. The main objective of RE is to provide a better understanding of the software systems current state, which can be used to correct (e.g. fix bugs), update (e.g. alignment with updated user requirements), upgrade (e.g. add new capabilities), or even completely re-engineer the system under study [START_REF] Rugaber | Model-driven reverse engineering[END_REF].

Generally, reverse engineering a software system is a time-consuming and error-prone process. Its difficult to predict how much time RE will require and there are no standards to evaluate the quality of the result of RE [START_REF] Rugaber | Model-driven reverse engineering[END_REF]. To overcome these difficulties, model driven reverse engineering (MDRE) [START_REF] Rugaber | Model-driven reverse engineering[END_REF][START_REF] Raibulet | Model-driven reverse engineering approaches: A systematic literature review[END_REF][START_REF] Bruneliere | Generic model-based approaches for software reverse engineering and comprehension[END_REF] has been proposed to enhance the traditional reverse engineering processes.

MDRE is the application of model driven engineering (MDE) principles and techniques to RE in order to generate relevant model-based views on legacy systems, thus facilitating their understanding and manipulation.

There have been several past works on MDRE which can be classified as two categories: specific and general solutions. This is determined depending on whether they aim to reverse engineer the system from a single technology and/or with a predefined scenario in mind (e.g., a concrete kind of analysis), or to be the basis for any other type of manipulation in later steps of the reverse engineering process [START_REF] Bruneliere | Modisco: A model driven reverse engineering framework[END_REF]. Manev et al. [START_REF] Manev | Facilitation of IoT software maintenance via code analysis and generation[END_REF] propose a tool, called ITACG (IoT software Analysis and Code-Generation tool), for performing reverse engineering and extraction. This is accomplished by scanning the source code of the target system and extracting architectural information from it, which is stored into a UML model. Umair Sabir et al. [START_REF] Sabir | A model driven reverse engineering framework for generating high level UML models from java source code[END_REF] present a MDRE framework named Src2MoF to generate UML structural and behavioral diagrams from the Java source code. In order to address several kinds of scenarios relying on different legacy technologies, Hugo Bruneliere et al. [START_REF] Bruneliere | Modisco: A model driven reverse engineering framework[END_REF] give an extensible and generic model driven reverse engineering: MoDisco. MoDisco has three layered architecture i.e. infrastructure, technologies and use case layers. It denes a basic meta-model approach for MDRE based on Knowledge Discovery Meta-model (KDM) specication to provide support for XML, JSP and Java. MoDisco only deals with structural aspects and does not support the MDRE for behavioral aspects from source artifacts.

Most of the existing works of MDRE mainly consider general domains such as desktop or business applications. In this paper, we consider MDRE in the domain of complex embedded systems, especially the safety-critical systems.

50

Complex embedded software systems are typically special-purpose systems developed for control of a physical process with the help of sensors and actuators.

They are often the systems requiring a deep combination of software, runtime operational system and hardware platform. Typical non-functional analysis of the requirements in this domain, such as safety, schedulability, and so on, 55 needs the modeling of architecture, functional behaviors and runtime. These characteristics already make it apparent that complex embedded systems d-iffer from desktop and business applications. Compared with the modeling languages used in the existing works of MDRE such as UML, AADL (Archi-tecture Analysis and Design Language) [START_REF]SAE, Architecture Analysis & Design Language[END_REF] is a powerful modeling language 60 for complex embedded system, which provides a unified formalism for the modeling of architecture, functional behaviors, and runtime. This paper proposes C2AADL_Reverse, a MDRE approach for safety-critical software development and verification. C2AADL_Reverse takes multi-task C source code as input, and generates AADL model of the legacy software systems. Moreover, when MDRE 65 exists in the domain of safety-critical systems, validation of the MDRE process and verification of the resulted models are highly desirable because such software systems have to undergo development regulations and certification restrictions.

Therefore, the reverse-engineered AADL components become the basis for applying MDD development in the same application domain, and should be 70 analyzed and verified.

Research Problems

Currently, there are several researches on AADL automatic code generation (i.e. forward engineering). For instance, OCARINA [START_REF] Hugues | From the prototype to the final embedded system using the ocarina aadl tool suite[END_REF] and RAMSES [START_REF] Rahmoun | Multi-objectives refinement of AADL models for the synthesis embedded systems (mu-RAMSES)[END_REF] support automatic code generation from AADL to C, Ada and Java. Regarding the reverse generation of AADL models, Wang et al. [START_REF] Wang | Studying on AADL-based architecture abstraction of embedded software[END_REF] propose an approach 75 for extracting AADL models from existing embedded software in order to reduce maintenance costs. In an effort to bridge the semantic and syntactic gaps between the two languages, they have defined a set of mapping rules from C to AADL models. For Integrated Modular Avionics (IMA) systems, Lesovoy et al. [START_REF] Lesovoy | Extracting architectural information from source code of AR-INC 653-compatible application software using CEGAR-based approach[END_REF] present an approach to extract the AADL models from source code 80 of ARINC 653-compatible application software. They apply the ideas of counterexample and path feasibility check to the task of extracting the architectural information from source code. As mentioned before, safety-critical software of-ten run on various embedded platforms, reverse engineering needs to deal with 85 the information such as static structure, dynamic run-time, and functional behavior. However, the existing approaches mainly deal with structural aspects instead of behavioral and run-time aspects of source artifacts. Safety-critical software systems are large and intricate, often constituting hundreds of components. Thus, the challenge is to be able to derive the information about 90 the functional behaviors and the runtime dynamics of a system. In particular, as multi-core processors are widely used in safety-critical software [START_REF] Salman | A systematic methodology to migrate complex real-time software systems to multi-core platforms[END_REF], the reverse engineering of multi-task synchronization, mutex, communication, and task scheduling has become an important problem.

Moreover, how to evaluate or measure a MDRE effort? On the one hand, problem. An approach to deal with the state-explosion problem is the use of compositional verification [START_REF] Posse | Contract-based specification and analysis of aadl models[END_REF][START_REF] Bensalem | Compositional verification for component-based systems and application[END_REF][START_REF] Cofer | Compositional verification of architectural models[END_REF] which leverages the structure of the system. The basic idea is to apply divide-and-conquer approaches to infer global properties of complex systems from properties of their components.

To overcome the above-mentioned research problems, we have implemented a complete framework for the proposed approach, that is C2AADL_Reverse, erties and the Assume Guarantee REasoning Environment (AGREE) [START_REF] Ghassabani | Proof-based coverage metrics for formal verification[END_REF][START_REF] Gacek | The JKind model checker[END_REF] to perform the compositional verification of the architecture.

.c

Main Contributions

The main contributions of the paper can be summarized as follows:

• A new MDRE approach named C2AADL Reverse: The transformation from multi-task C source code to AADL is divided into three parts:

-Structural aspect: the transformation from global variables, local variables, data types, function definitions and multi-task structures to AADL components;

-Behavioral aspect: the transformation from function and task execution behavior to AADL behavior annex [START_REF]SAE, Architecture Analysis and Design Language (AADL) Annex D: Behavior Model Annex[END_REF], which involve various types of branch statements, assignment statements, and function call statements;

-Run-time aspect: the transformation from multi-task communication, multi-task synchronization and mutex, and task scheduling to AADL execution-model properties.

• Validation and verification approach of C2AADL Reverse: Two types of activities are proposed to ensure the correctness of C2AADL_Reverse.

First, it is necessary to validate the reverse engineering process. Second, the generated AADL models should conform to desired critical properties.

We propose the verification of the generated AADL model by using the model checker UPPAAL to establish component-level properties and the AGREE environment to perform the compositional verification of the architecture. This combination of verification tools allows us to iteratively explore design and verification of detailed behavioral models, and to scale formal analysis to large models.

• The prototype tool : The C2AADL_Reverse prototype tool adopts a modular architecture, which is implemented based on the AADL open source environment OSATE [START_REF] Osate | Plug-ins for front-end processing of AADL models[END_REF], in which an intermediate model is proposed to facilitate the transformation from C source code to AADL.

• Case study: A real-world aerospace industrial case, the rocket launch control subsystem, is used to show the feasibility of the method presented in the paper.

Outline

The rest of this paper is organized as follows. Section 2 introduces the AADL language, the principle of compositional verification, and the AADL compositional verification tool AGREE. Section 3 presents the details of the transformation rules of the C2AADL_Reverse approach. Section 4 introduces the formal verification method of the generated AADL models. In section 5, we give the prototype tool. Section 6 illustrates a real-world aerospace industrial case study to show the eectiveness of the C2AADL_Reverse approach. Section 7 discusses related work and Section 8 provides concluding remarks and plans for future work.

Preliminaries

In this section, we rst provide an overview of the AADL language, and then introduce the principle of compositional verification and the AADL compositional verification tool AGREE.

AADL

The SAE Architecture Analysis and Design Language (AADL) is a textual and graphical language used to design and analyze the software and hardware architecture of embedded real-time systems. AADL is used to describe the structure of such systems as an assembly of software components mapped onto an execution platform, as shown in Fig. often faces the so-called state-explosion problem. An approach to deal with the state-explosion problem is the use of compositional verification which leverages the structure of the system. In these techniques, the verification of a composite system is reduced to the verification of its parts.

195

A well-known compositional approach is based on assume/guarantee contracts [START_REF] Posse | Contract-based specification and analysis of aadl models[END_REF][START_REF] Bensalem | Compositional verification for component-based systems and application[END_REF][START_REF] Cofer | Compositional verification of architectural models[END_REF][START_REF] Uriagereka | Guiding assurance of architectural design patterns for critical applications[END_REF] where each component is annotated with a contract consisting of an assumption specifying how the component expects its environment to behave, and a guarantee specifying the behaviour guaranteed by the com-ponent if the assumptions hold. The component implementations can be abstracted with contracts that specify the behavioral aspects that are relevant for the system-level properties. environment OSATE [START_REF] Osate | Plug-ins for front-end processing of AADL models[END_REF]. AGREE makes use of the AADL annex mechanism, named AGREE annex1 , to annotate AADL models with contracts. A contract contains a set of assumptions about the component's inputs and a set of guarantees about the component's outputs. The assumptions and guarantees may also contain predicates that reason about how the state of a component evolves over time. When AGREE is proving a system correct, it takes the specied behavior of the system along with its assumptions and the guarantees of any sub-components to verify the guarantees of the system being proven. In addition, AGREE can perform a standard iterative unrolling of the transition relation to find counter-examples, and if one is found, this counter example will be provided to the user. Please note that, AGREE currently mainly handles synchronous architectural models in which execution proceeds in a deterministic discrete sequence of steps [START_REF] Backes | Requirements analysis of a quad-redundant flight control system[END_REF]. Support for modeling components that execute asynchronously (or quasi-synchronously) will be added to AGREE.

C2AADL Reverse: The MDRE Approach for Generating AADL from C

C2AADL_Reverse takes multi-task C code as input, and generates an AADL model. It mainly contains two steps: code analysis and code-to-model transformation.

Code analysis

AADL is used to describe the model at the architecture level which is higher than C code. In order to fill in the syntactic and semantics gaps between AADL and C, we first build a code structure model from the source code, based on the parsing of the program. The simplified meta-model of the multi-task C code structure is given in Fig. 4. The top layer is Project. Each project can have several tasks and their communications, i.e., the task layer.

Code-to-Model transformation 260

The transformation from C to AADL consists of three parts: structure, behavior, and run-time information.

Structure transformation

As shown in Table 1, the structure transformation rules include the trans-

Behavior transformation

The behavior of a function or a task is defined by the statements inside the body of the function or task. The statements always include assignment, if, switch-case, for, while, and function call.

To be more intuitive, we present the description of the AADL behavior annex The transformations will be more complex in presence of inner statements.

We recursively apply transformations for each kind of statements until complex statements have been eliminated.

Run-time information transformation

It needs to consider the use of the platform's API in the transformation of the run-time information. Without loss of generality, in this paper we consider TI SYS/BIOS Real-time Operating System (SYS/BIOS) [START_REF] Instrument | TI SYS/BIOS v6. 33 real-time operating system users guide[END_REF] which is broadly used in the aerospace domain. The transformation rules mainly include task communication, task synchronization, and task scheduling. As shown in Fig. 7,

to

Validation and Verification Approach of C2AADL Reverse

365

Two types of activities are required to ensure the correctness of the reverse engineering approach proposed in this paper. First, it is necessary to validate the reverse engineering process. Second, the generated AADL models should conform to desired critical properties.

The validation and verification approach of C2AADL_Reverse is shown in

Model checking of leaf components

As the first step of compositional verification of a component-based system, it is necessary to prove that the behaviors of the leaf components satisfy their component-level contracts. Thus, in this paper, UPPAAL is used to verify AADL Behavior annex models. UPPAAL is a tool supporting the simulation and the verification of models defined by a set of timed automata communicating through synchronous channels and shared variables. The state transition of the Behavior annex is similar to the automata in the UPPAAL model. Therefore, the transformation from AADL behavior annex into UPPAAL model is straightforward.

First, the formal definitions of behavior annex and time automaton are stated as follows:

Definition 1 Behavior annex BA = S, S 0 , V, G, A, T , where:

• S is a set of states, in which a state may be qualified as initial state, final state, or complete state, or combinations thereof. A state without qualification will be referred to as execution state.

• S 0 is the initial state, S 0 ∈ S .

• V is a set of local variables.

• G (Guards) is a set of state transition conditions.

• A (Actions) represents the actions that need to be performed during state transitions.

• T is a set of state transitions, T ⊆ S × (G × A) × S.
Definition 2 A timed automaton [START_REF] Behrmann | A tutorial on uppaal, Formal methods for the design of real-time systems[END_REF] T A = L, l 0 , V, C, A, I, E , where:

• L is a set of locations;

• l 0 is the initial state, l 0 ∈ L; • V is a set of variables; • C is a set of clocks;
• A is a set of actions; It would be highly desirable if the formalism of the contracts expressed with the AGREE annex is consistent with the formalism of the properties used in UP-PAAL. Now, we manually generate the required component-level properties for verification in UPPAAL from the contracts. In the future, the translation between the AGREE contracts and the UPPAAL properties will be automated.

• I = L → B(C)
• E ⊆ L × A × B(C) × 2 C × L is

Compositional verification of AADL architecture model

To formally argue that the system satisfies its requirements, assume-guarantee contracts provide an appropriate mechanism for capturing the information from requirements or source code to reason about system-level properties. A contract specifies precisely the information that is needed to reason about the component's interaction with other parts of the system. Furthermore, the contract mechanism supports a hierarchical decomposition of the verification process that follows the natural hierarchy in the system model.

In the AGREE framework, it uses the AADL AGREE annex to specify the contracts of the component of each layer and the underlying formalism of the AGREE annex is the past-time operator subset of PLTL. We thus establish that the properties of the top-level system are proved given that the properties of the lowest layer i.e. leaf-level components are true (by model checking with UPPAAL shown in Section 4.1). As shown in Fig. 10, we take the LCU CT -ContOrd function of the rocket launch control system case study (Section 6.2) as example. It is reversed into a component containing two sub-components, and the sub-component contract is used to prove whether the upper-level component meets the requirements.

Prototype Tool Support

As mentioned above, the C2AADL_Reverse prototype tool adopts a modular architecture, which is implemented with the Eclipse plug-in technology.

Moreover, an intermediate model is proposed to facilitate C code information 460 extraction and AADL model generation. The tool is implemented as a plug-in of OSATE, which is shown in Fig. 13.

Intermediate model

495

It mainly provides three functionalities, such as project management, reverse engineering, code/model viewing and modification.

Project management

Code/model viewing and modification

Reverse engineering

Figure 13: C2AADL Reverse tool 6. Evaluation

Industrial Case Studies

The Rocket Launch Control System (RLCS) is a critical subsystem of the rocket launcher system (RLS). The function of the RLCS is to control the rocket to perform various operations and automatically execute rocket launch operation from launch command receipt to the the time when the rocket leaves the launcher. RLCS running on the launch control unit (LCU) computer and the LCU computer interacts with other modules of the RLS through a bus and/or a network. The RLCS can ensure the normal execution of the launch of the rocket through a series of hardwaresoftware interactions. A simplified architecture of RLCS is shown in Fig. 14, which mainly includes a management chassis, a launcher module, a power supply chassis, a power-on control module, a launch control unit detection board and a rocket control module. The power-on control module includes three sub-modules: power-on 1, power-on 2 and power-on 3.

In addition, the launcher module controls several devices such as exhaust cover and hatch, etc. The structure of RLCS can be divided into four layers, namely driver layer, driver management layer, process layer, and main control layer, among which:

• Driver layer: Describes the lowest-level module, which is responsible for the interaction of software and hardware, such as network communication, bus communication, and sending and receiving operations. The driver layer processes the received data and sends them to the tasks of the driver management layer.

• Driver management layer: The driver management layer analyzes and assembles the data corresponding to the different communication methods of the driver layer, and sends the assembled data to the process layer, such as the receiving and sending of power-on command;

• Process layer: The process layer parses the assembled data sent by the driver management layer into various commands. The functions in the main control layer are realized by the command processing in the process layer, such as power control, control of hatch covers, cartridges, thermal batteries, etc;

• Main control layer: The main control layer starts the task of the driver layer according to the command to perform corresponding function execution, such as sequence maintenance processing, electric blast tube in-spection, data receiving processing, sending control parameter processing, command execution processing, and so on.

Results and Analysis

C2AADL transformation

We use the C2AADL_reverse tool to reverse a part of source code of the RLCS to AADL model. This section takes power-on control as an example to analyze the generated AADL model from three perspectives: structure, run-time properties, and behaviors. 2. The reason why the coverage rate of the generated model does not reach 100% is that some codes cannot be expressed in the behavior annex, such as bit operation and type mandatory conversion, etc. However, AADL supports to describe the behaviors of thread/subprogram with source code directly, through properties such as Source Language, Source Name, Source Text, etc. Thus, it is used as a complementary of the results of C2AADL_Reverse.

Validation and verification

1) Validation of C2AADL reverse

We refined the generated AADL model by adding processors, virtual processors and scheduling mechanisms, and used the C code generator of OCARINA [START_REF] Mkaouar | An ocarina extension for AADL formal semantics generation[END_REF] to produce an executable version from our model of RLCS. We can compare the generated codes with the original ones by using code reviewing and executing. Fig. 16(a) shows that the execution sequence of the tasks in the generated code is consistent with the original code. Fig. 16(b) shows the comparison between a part of the original codes and the generated codes. OCARINA mainly considers the POK operational system [START_REF] Delange | an ARINC653-compliant operating system released under the BSD license[END_REF]. The pok buffer send and pokbuffer receive used for inter-task communication in the POK OS are consistent with the Mailbox post and Mailbox pend in SYS/BIOS. Thus, the C2AADL_-Reverse method can generate adequate models.

(LCS), Exhaust cover control (ECC), Rocket hatch control (RHC), and Rocket power-on control (RPC), etc. These requirements were developed hierarchically following the system architecture. A simplified requirements hierarchy is shown in Fig. 17.

Launch control system (LCS)

Power-on control (SR4) Exhaust cover control (SR1)

Hatch control (SR2) SR4C1 SR4C2 SR4C3 SR4C1H1 SR4C1H2 SR4C1H3 WD1 PF1 CPO1 SSPO1 SR1C1 SR1C2 WDSR2 SSR2 SR2C1 SR2C2
…… Figure 17: A simplified LCS requirements hierarchy 580

In order to illustrate the relationship between the requirement decomposition and the system structure decomposition of the RLCS, this paper takes the rocket power-on control (SR4) as an example to illustrate the further detailed decomposition of the requirement. Rocket power-on control consists of three sub-components: rocket power-on 1 control (SR4C1), rocket power-on 2 control 585 (SR4C2) and rocket power-on 3 control (SR4C3), and each power-on control is used to control different battery components. For instance, power-on 1 can be decomposed into rocket self-check result query function (SRC1H1), rocket power-on 1 function (SR4C1H2), and rocket type identification code query func-tion (SRC1H3). Furthermore, SR4C1H2 can be divided into watchdog function 590 1 (WD1), power-off 1 (PF1), continuous power-on 1 (CPO1), and single-step power-on 1 (SSPO1). AGREE currently mainly handles synchronous architec-tural models in which execution proceeds in a deterministic discrete sequence of steps, so we consider the model which execute synchronously. Based on the requirement decomposition, the requirements of each component are formalized 595 as contracts. Several contracts are shown in the following paragraphs.

34

of the exhaust cover control subsystem. The system-level requirement of the exhaust cover control subsystem is SR1, which can be decomposed into exhaust cover switch cover processing (SR1C1) and exhaust cover switch cover control (SR1C2). SR1C2 can be divided into exhaust cover self-inspection watchdog (WDSR2) and exhaust cover self-inspection function (SSR2). The compositional reasoning result is shown in Fig. 19. According to the results of compositional verification, it can be known that the compositional verification of the exhaust cover control subsystem involves 13 contracts, among which the system-level contract SR1G is satisfied. The major objective of this subsection is to evaluate the eectiveness of our approach. This objective is decomposed into the following research questions:

645

RQ1: Can we use our approach to obtain AADL models for dierent case studies from the aerospace domain?

Rationale: As a model generation approach, C2AADL_Reverse should be able to generate various kinds of legacy systems that implemented in C language.

To answer the question, we apply our approach to several case studies from the aerospace domain. Table 5 provides the name, the short description, and the number of AADL models from three industry cases. For the description of Case-A, see Section 6.1. Case-B is is a core system supporting the orbiting operations of spacecraft, which undertakes the tasks of determining and controlling spacecraft attitude and orbit. For Case-C, the radar information processing subsystem refers to a system that uses modulated wave forms and directional antennas to transmit electromagnetic waves to a specic airspace in space, and extract guidance information from the received echoes. Moreover, since C and AADL are broadly used, the approach proposed in this paper can be also applied into other safety-critical domains such as automotive system.

RQ2: Can our approach generates high quality models when compared with other approaches?

This RQ checks how much information our approach can model. To answer 665 the question, we compare the C2AADL_reverse tool with several existing tool-s:

MoDisco [START_REF] Bruneliere | Modisco: A model driven reverse engineering framework[END_REF], fREX [START_REF] Bergmayr | fREX: fUML-based reverse engineering of executable behavior for software dynamic analysis[END_REF], RE-CMS [START_REF] Trias | RE-CMS: a reverse engineering toolkit for the migration to CMS-based web applications[END_REF], Src2MoF [START_REF] Sabir | A model driven reverse engineering framework for generating high level UML models from java source code[END_REF], srcYUML [START_REF] Decker | A tool for efficiently reverse engineering accurate uml class diagrams[END_REF], and Wang [START_REF] Wang | Studying on AADL-based architecture abstraction of embedded software[END_REF]. In

Y N N N Src2MoF Java UML Y Y N N srcYUML C++ UML Y N N N Wang [13] C AADL Y N N N C2AADL re- verse C AADL Y Y Y Y

Related Work

Model-driven reverse engineering

We distinguish the exiting work on MDRE into two main families: specic vs. general purpose solutions.

General Frameworks of MDRE

Fleurey et al. [START_REF] Fleurey | Model-driven engineering for software migration in a large industrial context[END_REF] propose a semi-automatic round-trip model-driven migration and modernization process for the migration of large industrial software.

This migration solution includes the automatic analysis of the source code, the generation of abstract models into target platform models, and generation of code for the target system.

The approach proposed by Favre et al. [START_REF] Favre | Formalizing mda-based reverse engineering processes[END_REF] is placed within the context of obtaining models from object-oriented code and formal techniques at the metamodel level to maintain consistency in reverse engineering process, according to the MDA standard. This approach exploits static and dynamic analysis to generate PSMs and PIMs from code and to analyze the consistence of the performed transformations from code to models and between models.

A semi-automatic approach for MDRE to extract business rules from Java applications is presented in [START_REF] Normantas | Extracting business rules from existing enterprise software system[END_REF]. The authors address issues concerning the extraction of knowledge from software artifacts and the representation of the extracted knowledge into the KDM metamodel with the objective to abstract the business logic implemented in a system.

MoDisco [START_REF] Bruneliere | Modisco: a generic and extensible framework for model driven reverse engineering[END_REF] and extraction. This is accomplished by scanning the source code of the target system and extracting architectural information from it, which is stored into a UML model.

Regarding the reverse generation of AADL models, Wang et al. [START_REF] Wang | Studying on AADL-based architecture abstraction of embedded software[END_REF] proposes an approach for AADL models extraction from existing embedded software and 725 reduce maintenance costs. In an effort to bridge the semantic and syntactic gaps between the two languages, they have defined a set of mapping rules from C to AADL models. However, this method does not consider behavior and runtime information, and does not validate the reverse process. In Integrated Modular Avionics (IMA) systems, Lesovoy et al. [START_REF] Lesovoy | Extracting architectural information from source code of AR-INC 653-compatible application software using CEGAR-based approach[END_REF] proposes an approach to extract the 730 AADL models from source code of ARINC 653-compatible application software.

This approach applies the ideas of counterexample and path feasibility check to the task of extracting the architectural information from source code.

However, since safety-critical software often run on various embedded plat-forms, reverse engineering needs to deal with the information such as static

Verification of AADL models

A number of works have been proposed respect to formal analysis of AADL 745 models. However, they always consider a subset of AADL because that the AADL language is very rich. In [START_REF] Johnsen | Automated verification of AADL-specifications using UPPAAL[END_REF], a methodology for translating AADL to UPPAAL has been proposed. This work is primarily focused on flow and deadlock analysis rather than behaviorial specification. Chkouri et al.,[START_REF] Chkouri | Translating AADL into BIP -application to the verification of real-time systems[END_REF],

proposed a method for translating AADL models into BIP, which makes it possi-750 ble to make use of the BIP toolset for verification. Bodeveix et al. [START_REF] Bodeveix | Towards a verified transformation from AADL to the formal component-based language FIACRE[END_REF] proposed a verification tool chain for AADL models through its transformation to the Fiacre language, and prove the correctness of the transformation from AADL into Fiacre. Yang et al. [START_REF] Yang | From AADL to timed abstract state machines: A verified model transformation[END_REF][START_REF] Hu | Exploring AADL verification tool through model transformation[END_REF] proposed a methodology for translating AADL to TASM, through which the properties of AADL can be analyzed by the toolset de-755 veloped for TASM. Yu et al. [START_REF] Yu | Polychronous modeling, analysis, verification and simulation for timed software architectures[END_REF] proposed to apply the POLYCHRONY toolset, based on the synchronous language SIGNAL, for timing modeling, analysis and validation of AADL. Hugues et al. [START_REF] Mkaouar | Towards a formal specification for an AADL behavioural subset using the LNT language[END_REF][START_REF] Mkaouar | A formal approach to AADL model-based software engineering[END_REF] proposed the formal verification of real-time systems modelled with the AADL language and its behaviour annex, and defines a formal semantics of an AADL behavioural subset using the LNT 760 (Lotos NT) language.

Compositional verification has attracted significant research attention because of its viability as a scalable technique for reasoning about complex sys-tems.

Backes et al. [START_REF] Backes | Requirements analysis of a quad-redundant flight control system[END_REF] applied a compositional verification approach to a level properties by using the Simulink Design Verifier to establish componentlevel properties and AGREE to perform the compositional verification of the architecture, and verifies a realistic medical cyber-physical system. Cofer et al.

[18] described the compositional reasoning framework for proving the correctness of a system design, and provide a proof of the soundness of their compositional reasoning approach. An aircraft flight control system is provided to illustrate the method. Posse [START_REF] Posse | Contract-based specification and analysis of aadl models[END_REF] proposed an AADL annex sub-language for annotating components with assume/guarantee contracts and a prototype verier that performs compositional analysis. They also provide a prototype plug-in for OSATE supporting an annex language which is called AGCL.

Compared with existing works, this paper focuses on the verification of the resulted models of the reverse process from C to AADL. Because of the broadly use of UPPAAL in the verification of functional behaviors of components, we assemble verification of system-level properties by using UPPAAL to establish component-level properties and the AGREE environment to perform the compositional verification of the architecture. In [START_REF] Whalen | Hierarchical multi-formalism proofs of cyberphysical systems[END_REF], the authors consider a hierarchical multi-formalism proofs of cyber-physical systems by using AGREE and UPPAAL. However, they focus on multiple abstraction layers. AGREE is used in the high-level analysis, while UPPAAL is applied in the low-level one.

Therefore they propose a translation from AGREE contracts to timed automata.

In this paper, we combine AGREE with UPPAAL in the same abstraction layer, in which the model checker UPPAAL establishes the behaviors of leaf components and the AGREE environment performs the compositional verification of the architecture.

Conclusion and Future Work

This paper has presented a model-driven reverse engineering approach for safety-critical software development and verification, namely C2AADL_Reverse. We will further carry out the following future work:

• The source code cannot explicitly express non-functional properties of the software system (such as period, execution time, resource consumption, and so on). At present, we apply third-party dynamic tools (such as WCET analysis tools) to measure timing properties and add them to the corresponding AADL model.

• Inspired by the restricted natural language approach proposed in our previous work [START_REF] Wang | An approach to generate the traceability between restricted natural language requirements and AADL models[END_REF], the automatic transformation from natural language requirements into AGREE contracts is currently being developed. As well, the translation from AGREE contracts to TCTL properties in UPPAAL will be also automated.

• We are considering the extension of AGREE to support for modeling components that execute asynchronously (or quasi-synchronously), and formalizing the reasoning rules in the theorem prover Coq [START_REF] Chlipala | Certified Programming with Dependent Types: A Pragmatic Introduction to the Coq Proof Assistant[END_REF].

  we can use the generated model of MDRE to produce another version of the 95 original software and make the comparison between the two versions to validate the MDRE process. On the other hand, automatic formal verification techniques such as model-checking can be used to analyze the behaviours of the generated model. Since the increasingly size of the source code, formal verification of reverse-engineered AADL models often faces the so-called state-explosion

as shown in Fig. 1 .

 1 It includes five phases, (1) analysis of the original source code, (2) extraction of an intermediate model, (3) generation of an AADL model, (4) validation of the C2AADL process, and (5) formal verification of the generated AADL model. Compared with the existing AADL RE method, this paper considers more reverse constructions including AADL component structure, behavior, and multi-threaded run-time information. For the validation of the reverse process, we generate a second version of the original software and compare the two versions of code. Moreover, we propose the verification of the generated AADL model by using UPPAAL to establish component-level prop-

Figure 1 :

 1 Figure 1: The framework of C2AADL Reverse

Fig. 3 Figure 3 :

 33 Figure 3: An example of compositional verification

Figure 4 :

 4 Figure 4: Simplified meta-model of multi-task C code structure

Figure 5 :

 5 Figure 5: An example of the construction of code structure model

265

  formations from basic data types and composite data types to AADL data components, from function definitions to AADL subprogram components, and from task structures to AADL thread components. Since variables in C language and data components in AADL describe similar data types, we map variables to data components. In order to make the structure of AADL clearer, we encapsulate the data components into two types of packages: one package represents 270 the basic data type named Base Types, and the other package represents the extended data type (composite data type) named User Define. For basic data types, they can be mapped into AADL data types straightforward, while the pointer is represented as a data access feature of component. AADL can also express signed and unsigned integers. For composite data types with internal 275 implementations, the elements of the composite data are mapped to subcomponents declared in the implementation of data component. In addition, the function definitions and task structures are mapped to AADL subprogram and thread components respectively.

285 as a graphicFigure 6 :

 6 Figure 6: The transformation rules of behavior

Figure 7 :

 7 Figure 7: The transformation rules of run-time information

Fig. 8 .Figure 8 :

 88 Fig.8. At first, the validation of the reverse engineering is performed by code

  assigns invariants to locations, B(C) is the set of conjunctions over simple conditions of the form x c or xy c, where x, y ∈ C, c ∈ N and ∈ {>, ≤, =, ≥, <};

  a set of edges between locations with an action, a guard and a set of clocks to be reset. On the one hand, BA is transformed into TA. Globally, the state set of BA is transformed to the states of TA (BA.S → T A.L); the initial state of BA corresponds to the initial state of TA (BA.S 0 → T A.l 0 ); the local variable of BA is transformed to the variable in TA (BA.V → T A.V ); the state transition of BA is transformed to the state transition of TA (BA.T → T A.E); the operation on time of BA such as computation() is converted into the invariant(I) of the state in TA, for instance computation(5ms) is represented by cl <= 5; there are two clocks in the C of TA: cl and global cl , cl records the time of the operation performed by each state transition and global cl records the time when the automaton is executed.On the other hand, the Assume of the contract ( Assume, Guarantee ) of the component is transformed into the initialization operation when the time automaton starts to execute, and the Guarantee of the contract is transformed into the property written in TCTL language that needs to be verified by UP-PAAL. As mentioned in Section 2.2, AGREE makes use of the AGREE annex to annotate AADL models with contracts. The underlying formalism of the A-GREE annex language is a subset of Past-Time Linear Temporal Logic (PLTL).

Fig. 9

 9 Fig.9shows an example of the transformation from BA to UPPAAL. It verifies that the component meets the contract without deadlock and timeout.

Figure 9 :

 9 Figure 9: The transformation of BA to UPPAAL

Figure 10 :

 10 Figure 10: The framework of compositional verification for architecture

As shown in Section 3 . 1 ,

 31 the meta-model of the multi-task C code structure represents the abstract concepts of source code, which is helpful for the description of the transformation rules. Here, we propose an intermediate model named 465 CAInterM to facilitate the implementation of the C2AADL_Reverse tool. The intermediate model contains the elements of the meta-model of the multi-task C code structure, the file structures, and the connection information among tasks or functions. The intermediate model CAInterM is shown in Fig.11.

Figure 11 : 5 . 2 .Figure 12 :

 115212 Figure 11: Intermediate model: CAInterM

Figure 14 :

 14 Figure 14: The rocket launch control system

Fig. 15 (

 15 Fig.15(a) shows the top-level system structure of the generated AADL model. The AADL model is consistent with the RLCS system introduced in Section 6.1, several thread components (such as LCU DM NetRecv,FRAME LCU DT -NetDataAna, etc) implement the power-on function of the system through port communication. Fig.15(b) shows the communication between multiple components, thread LCU DM NetRecv receives data from the network and sends it to the thread FRAME LCU DT NetDataAna via the NetRecv. Fig.15(c) shows an individual component, including input/output parameters, data access and behavior annex.

Figure 15 :

 15 Figure 15: Generated AADL models of RLCS

Requirement 1 :Requirement 2 :

 12 When the task FRAME LCU DT NetDataAna of the driver management layer receives the data sent by the task LCU DM NetRecv of the driver layer, it sends the encapsulated data to the main control task (contType=0 means single-step power-on, MisNnum=1 means power-on 1). It is formalized in AGREE annex (R1G) shown as follows: 600 thread FRAME_LCU_DT_NetDataAna features NetRecv : in data port Base_Types::Integer; MainCont : out data port Globle_Define::UNION_Cont_Order.impl; properties Priority=>5; annex agree {** property judge_MainCont= MainCont.Info=1 and MainCont.STRUCT_Bits.MisNnum=1 and MainCont.STRUCT_Bits.ContType=0; assume "A:FRAME_LCU_DT_NetDataAna receive data from NetRecv" : NetRecv = 0; guarantee "G:FRAME_LCU_DT_NetDataAna send data to MainCont" : judge_MainCont; **}; end FRAME_LCU_DT_NetDataAna; The task LCU CT MainCont of the main control layer, after receiving the data of single-step power-on 1, calls related functions to process the data, sends an order to the task FRAME LCU DT 289ASend of the driver management layer, and sends a data with the same MisnNum to the task FRAME LCU DT AnaOrder of the driver management layer. It is 605 formalized as the contract R2G:

Figure 19 :

 19 Figure 19: The compositional verification results of the exhaust cover control subsystem model

  transitions T_0: s0-[msg!="NULL"]->s1; T_1: s1-[msg.STRUCT_Bits.ContType=0]->s2; T_2: s2-[msg.STRUCT_Bits.ContOrder=1]->s3{ FRAME_PC_MisReady1!(msg.STRUCT_Bits.MisNnum,Car_State); LCU_WD_PlatInfo!(msg.STRUCT_Bits.MisNnum)}; T_3: s3-[msg.STRUCT_Bits.ResChk=1]->s4{ LCU_WD_PowerOn1Chk!(msg.STRUCT_Bits.MisNnum)}; T_4: s3-[not(msg.STRUCT_Bits.ResChk=1)]->s4; T_5: s2-[not(msg.STRUCT_Bits.ContOrder=1)]->s4; T_6: s1-[not(msg.STRUCT_Bits.ContType=0)]->s4; T_7: s0-[not(msg!="NULL")]->s4; **}; end LCU_CT_ContOrd.impl;

Figure 20 :

 20 Figure 20: Verification of leaf component based on UPPAAL

  720

  735 structure, dynamic run-time, and functional behavior. In particular, as multi-core processors are widely used in safety-critical software, the reverse engineering 42 of multi-task synchronization, mutex, communication, and task scheduling has become an important problem. Moreover, when MDRE exists in the domain of safety-critical systems, validation of the MDRE process and verification of the resulted models are highly desirable because such software systems have to undergo development regulations and certification restrictions. To the best of our knowledge, this paper presents a first effort on the validation and verification of the reverse process from C to AADL.

  Compared with the existing works, C2AADL_Reverse considers more reversed construction including AADL component structure, behavior, and multi-threaded run-time information. Moreover, when MDRE exists in the domain of safetycritical systems, validation of the MDRE process and verification of the resulted models are highly desirable because such software systems have to undergo development regulations and certification restrictions. We use reverse reverse engineering to validate the reverse engineering process, and verify the generated AADL models by using the model checker UPPAAL to establish componentlevel properties and the AGREE environment to perform the compositional verication of the architecture. To the best of our knowledge, this paper presents a first effort on the validation and verification of the reverse process from C to AADL. Finally, the effectiveness of C2AADL_Reverse is demonstrated using a real-world aerospace case study.

Behavior annex Hardware components Execution model Software component Software component functional behavior

		Initial	Condition Action	Final/Return
						C	Ada
		Thread group	Thread		Thread	Thread
			Subprogram	Data		Subprogram group
		Process			Process
			Scheduling, Communication, Dispatch, Mode Change ,etc.
		Virtual	Virtual	Virtual	Virtual	Virtual	Virtual
		Processor 1	Processor 2	Processor n	Processor 1	Processor 2	Processor n
		Multi-core processor		Multi-core processor
				AADL Bus
		Memory	Device 1	Device n
				System
		(such as Processor, Virtual Processor, Memory, Device, Bus, and Virtual Bus)
		Figure 2: AADL basic modeling concepts
		and their connections are used to describe the hardware architecture of the sys-
		tem. In addition, AADL describes the run-time state of the system through 2.2. Compositional Verification and AGREE
	185	properties such as Dispatch, Communication, Scheduling, Mode Change, etc., Complex embedded systems are always hierarchically organized by using
		that is, the Execution model. Finally, software components, hardware compo-component-based architecture. Automatic formal verification techniques such
		nents, and execution models are combined with System components to establish as model-checking can help to analyze the behaviours of such systems. For
		a hierarchical system architecture model. Furthermore, it contains an extension instance, thanks to a model-checking tool, one can create a model and analyze
		mechanism (called an annex, e.g. the Behavior annex [21][23]) that can be used all of the behaviors of the components in the architecture model. Actually,
	190	to extend the language to support additional features. The Behavior annex is most of the time, the architecture model is flattened. Nevertheless, doing so,
		defined for the refinement of thread/subprogram behaviors including functional
		behaviors of thread/subprogram and dispatch behaviors of thread ports. The
		behaviors of thread/subprogram can also be defined by the traditional program-
		ming language such as C, Ada, etc.

2. Software components and their connections are used for software architecture modeling, including Data, Thread, Thread Group, Process, and Subprogram components. Hardware components

Table 1 :

 1 The transformation rules of structure information

	C language	AADL
	int a; int *a
	char a; char *a
	bool a; bool *a
	float a; float *a;
	struct name a;
	struct name *a;

Base Types::Integer; requires data access Base Types::Integer Base Types::Character requires data access Base Types::Character Base Types::Boolean requires data access Base Types::Boolean Base Types::Float requires data access Base Types::Float User Defined::struct name.impl; requires data access User Defined::struct name.impl; data Integer 32 extends Integer properties Data Model::Number Representation ⇒ Signed; end Integer 32; data Unsigned 32 extends Integer properties Data Model::Number Representation ⇒ Unsigned; end Unsigned 32; struct dataname{ type spec var name; }; enum; union; ...... data dataname properties Data Model::Data Representation⇒ (Struct/Union/Enum); end dataname data implementation dataname.impl subcomponents var name: data package name::type spec

Table 2 :

 2 The statistics of the AADL model of the launch control

	system
		AADL mod-	Threads Subps Coverage
		el (line)
	Exhaust cover control	1800+	4	14	93%
	Rocket hatch control	1600+	3	12	92%
	Control of rocket launch prepa-	1600+	4	11	94%
	ration/cancellation
	Rocket power-on control	4000+	18	32	93%
	Control of rocket hatch unlock-	2000+	6	13	95%
	/lock

Table 3 :

 3 The statistics of verification results of launch control sys-to illustrate the verification of AADL leaf components based on UPPAAL. As shown in Fig.20, the left side is the AADL model, and the corresponding UPPAAL model is given on the right side. We create the initial function (start()) of the time automata according to the Assume in the AADL contract, and write the verification formula of UPPAAL according to the Guarantee in Fig.10. The verification results are shown in Table4.

	tem
		Contract Verification	Correctness
		time
	Exhaust cover control (SR1)	13	5s	13/13

635

Table 4 :

 4 The verification results of the LCU CT ContOrd function

		TCTL		Description	T/F
	A[]!deadlock			No deadlock		T
	E <> P rocess.T ime out		Reachability,	automata	F
				will eventually go to the
				Time out state
	A[]	Process.s4	imply	Functional correctness, af-	T
	(ProAnaContOrder==1	and	ter the LCU CT ContOrd
	DT289ASend.STRUCT Bit-		is executed, two data are
	s.MisNnum==1)		sent through the port
	6.2.3. Analysis
	• Research Questions

Table 5 :

 5 Application in different case studies

	Case	Description	Number of Require-
			ments
	Case-A	Rocket Launch Control Subsys-	17600+
		tem (RLCS)
	Case-B	Autonomous Guidance, Naviga-	13400+
		tion and Control (AGNC)
	Case-C	Radar Information Processing	11700+

Table 6

 6

		, we delimit six parameters to compare the most renowned re-
		searches. In 1) Source Artefact and 2) Target Model, as we seen, most
		of the existing works of MDRE mainly consider general domains, thus sever-
	670	al tools support the reverse engineering from object-oriented source code into
		UML. This paper focuses on the domain of complex embedded systems, espe-
		cially the safety-critical systems. 3) Structural indicates whether the MDRE
		methods can generate structural models. 4) Behavioral indicates whether the
		MDRE methods can generate behavioral models. fREX and Src2MoF mainly
		consider UML activity diagram, while we consider AADL Behavior annex. 5)

675

Runtime indicates whether the MDRE methods can generate execution models, which are important for the embedded software systems. 6) V&V indicates validation and verification.

Table 6 :

 6 Comparisons with a part of MDRE tools

	Parameters	Source	Target	Structural Behavioral Runtime V&V
	Tools	Artefact	Model
	MoDisco	Java,JSP UML	Y	N	N	N
	fREX	Java	UML	N	Y	N	N
	RE-CMS	PHP	AST of
			PHP

  is an extensible and generic Eclipse plug-in for model driven reverse engineering proposed and implemented by Bruneliere et al., based on model discovery and model understanding. MoDisco has three layered archi-

		7.1.2. Specific MDRE Solutions
		Several past works have already described how to migrate from a particular
	technology to another one using dedicated components.
	710	In the automotive industry, General Motors (GM) has been using a custom-
		built, domain specic modeling language, implemented as an internal proprietary
		metamodel, to meet the modeling needs in its control software development. Se-
		lim et al. [37] explores applying model transformations to address the challenges
		in migrating GM-specic, legacy models to AUTOSAR (AUTomotive Open Sys-
		tecture i.e. infrastructure, technologies and use case layers. It denes a basic
		metamodel approach for MDRE based on Knowledge Discovery Meta-model
		(KDM) specication to provide support for XML, JSP and Java.

715

tem ARchitecture) equivalents. They have built and validated a model transformation using the MDWorkbench tool, the Atlas Transformation Language, and the Metamodel Coverage Checker tool.

In the field of IoT, Manev et al.

[START_REF] Manev | Facilitation of IoT software maintenance via code analysis and generation[END_REF]

proposes a tool, called ITACG (IoT software Analysis and Code-Generation tool), for performing reverse engineering

Available at: http://github.com/smaccm/smaccm

realistic avionics system and demonstrated the effectiveness of the AGREE tool in performing this analysis. Murugesan et al.[START_REF] Murugesan | Compositional verification of a medical device system[END_REF] assembled proofs of system

Acknowledgements

This work is supported by National Natural Science Foundation of China (62072233) and Aviation Science Fund of China (201919052002).