Gaëtan Desrues

Hervé Delingette

Maxime Sermesant

Towards Hyper-Reduction of Cardiac Models using Poly-Affine Deformation

Keywords: Model Order Reduction, Finite Elements Method, Affine Transformation, Meshless

This paper presents a method for frame-based finite element model in order to develop fast personalised cardiac electromechanical models. Its originality comes from the choice of the deformation model: it relies on a reduced number of degrees of freedom represented by affine transformations located at arbitrary control nodes over a tetrahedral mesh. This is motivated by the fact that cardiac motion can be well represented by such poly-affine transformations. The shape functions use then a geodesic distance over arbitrary Voronoï-like regions containing the control nodes. The high order integration of elastic energy density over the domain is performed at arbitrary integration points. This integration, which is associated to affine degrees of freedom, allows a lower computational cost while preserving a good accuracy for simple geometry. The method is validated on a cube under simple compression and preliminary results on simplified cardiac geometries are presented, reducing by a factor 100 the number of degrees of freedom.

Introduction

Patient specific cardiac modelling is important for understanding pathologies, planning a therapy or rehearsing a surgery. In order to be used in routine by clinicians, models have to be fast enough and provide an accurate solution to given boundary and initial conditions. The finite element method is a classical approach to solve physically based deformation problems where the computational domain depends on an underlying mesh, implying that the computation time and the accuracy will depend on the mesh discretization. Particle-based meshless frameworks have recently been proposed [START_REF] Gilles | Frame-based interactive simulation of complex deformable objects[END_REF][START_REF] Lluch | Breaking the state of the heart: meshless model for cardiac mechanics[END_REF] and have the advantage of being less dependent from the topology. Particles, called control nodes, can be placed arbitrarily. It has also been shown in the literature [START_REF] Mcleod | Spatio-Temporal Tensor Decomposition of a Polyaffine Motion Model for a Better Analysis of Pathological Left Ventricular Dynamics[END_REF][START_REF] Zhou | 3D harmonic phase tracking with anatomical regularization[END_REF] that cardiac motion can be well represented by mixing affine tranformations located at the center of sub-regions of the American Heart Association (AHA) segments, see Fig. 1. Furthermore, the decomposition into a mixture of affine transformations provides clinically interpretable parameters related to cardiac strain.

The key idea of this paper is to combine these approaches to propose an efficient reduced model. Compared to other reduced models using statistical learning of a reduced basis, this approach builds upon physiological knowledge of the cardiac function to propose a novel reduction.

We first define the degrees of freedom (DOFs), regions and integration points. The internal forces resulting from the deformation of the solid are computed based on the minimization of an elastic energy density, itself resulting from the relation between strain and stress. Fig. 2 presents the famework used in our simulation of elastic solids. The different elements and the link between them through mappings will be described in further sections. Fig. 2. Three levels continuum mechanics. The shape functions map the degrees of freedom, supported by the control nodes, to the deformation gradients located in the sample regions (grey square). Then, the material law computes the internal forces from the stress and strain. These forces are finally mapped back to the degrees of freedom. This process is repeated at each simulation step. Source: Flexible plugin documentation (Sofa-framework).

Frame-based deformation

Our scheme relies on a reduced number of affine frames as degrees of freedom acting on a tetrahedral mesh V. If N is the number of vertices of V, the classical approach is to define a displacement field based on the displacement of each vertex, which makes this field of dimension 3N . Yet the cardiac motion usually does not require such high number of DOFs since the motion of neighboring vertices is highly correlated. This is why we propose to discretize the displacement field as a combination of n affine frame motions which makes it a 12n dimensional space. Instead of the 3-coordinate vector displacement field classically used, we propose to use a poly-affine displacement field for which at each control node q i , 1 ≤ i ≤ n, the 12 DOFs are the 12 coefficients of the affine transformation T(q i ) = T i ∈ M 3×4 (3 coefficients for the 3D translation and 9 for the linear transformation combining rotation, scaling and shearing).

The whole computational domain Ω is partitioned into n non overlapping regions Ω i , centered on q i , each region being defined as a set of tetrahedra.

At a point x in the domain, the local affine transformation T(x) is defined as a weighted sum of DOFs T i :

∀x ∈ Ω, T(x) = n i=1 φ i (x)T i
The weights φ(x) are the shape functions and represent the influence of the control frames over the domain. Their definition depends on the mesh geodesic signed distance D i (x) of a point x with respect to region i:

∀x ∈ Ω, ∀i ∈ 1, n , D i (x)    = 0 on the border of Ω i < 0 inside Ω i > 0 outside Ω i
This distance is computed after solving an Eikonal problem with fast marching method [START_REF] Pernod | A multi-front eikonal model of cardiac electrophysiology for interactive simulation of radio-frequency ablation[END_REF]. It is normalized as Di (x) such that its minimum value is -1 at its center q i . If σ a is the sigmoid function, the shape function φ i is defined as:

φ i (x) = σ a Di (x) , σ a (x) = 1 1 + e ax
The parameter a controls the overlap of the kernel function over neighboring regions. We observed better results with a = 2.25, thus φ i (q i ) = 0.90. It has been noticed that a big overlap (a small) leads to a smooth solution while a small one leads to more independent regions, in both cases introducing unrealistic values in the solution. An example with a = 2.25 is presented in Fig. 3(c).

Numerical integration method

The elastic energy of a deformable solid is the work done by the elastic forces between the undeformed and deformed positions, integrated across the whole domain. The numerical integration of the elastic forces is classically performed at the tetrahedron level in regular finite element methods. In our approach, we partition the domain Ω into M integration regions V m consisting of a set of tetrahedra. In the spirit of our affine frame control nodes, we use a high order integration rule called elastons introduced in [START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF] and generalized in [START_REF] Gilles | Frame-based interactive simulation of complex deformable objects[END_REF]. The classical integration approach assumes a constant force within each integration region V m . The elaston framework relies on a first order Taylor expansion of the field to be integrated in order to reach a higher level of accuracy. More precisely, a field f (x) is locally approximated as f (x) ≈ F p where F is a vector containing the value of the field and its derivatives at the center of V m , and p is the polynomial basis of order α in dimension d, for example [1, x, y, z] at order 1, in 3D. Finally, the integration is performed as :

Ω f (x) dx = m V m f (x m ) dx = m F m V m p dx
In our approach, the strain and stress are approximated by a first order polynomial function in the vicinity of each integration sample which is estimated thanks to a Generalized Moving Least Square (GMLS) interpolation scheme [START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF].

Kinematic Mappings

The mappings are templated functions allowing to project forces, displacement and their derivatives between the affine frame DOFs and the tetrahedral mesh. We present the well known linear blend skinning and the GMLS mapping, very suitable with our first order framework.

Position and Force Mappings

The mesh vertices are displaced based on a linear blending of the affine frame DOFs using the shape functions defined in section 2.

Let xj be the rest position of the j th mesh vertex. Its deformed position x j is computed as a function of the affine transformations T i :

∀j ∈ 1, N , x j = n i=1 φ i (x j )T i xj 1
Conversely, to compute the resulting force F i applied to the i th frame DOF, we gather all the forces f j applied to the vertices with the weighted sum:

∀i ∈ 1, n , F i = N j=1 φ i (x j ) f j

GMLS mapping

The GMLS method was shown [START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF][START_REF] Fries | Classification and overview of meshfree methods[END_REF] to be an efficient method to approximate globally a function from sparse discrete values with a minimisation problem. This problem seeks new shape functions N i (x) representing accurately the function over the whole domain.

T(x) = n i=1 T i N i (x)
In our framework, this mapping is used to compute the deformation gradients from vertices position at integration points. See [START_REF] Martin | Unified simulation of elastic rods, shells, and solids[END_REF] for an analytical formulation of N i (x).

Results

The proposed method was tested on different geometries, increasing in complexity. The simulation is performed on an Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz laptop, using the SOFA framework. The method widely uses the Flexible plugin approach of Fig. 2.

Cube under simple compression

First, we simulate the compression of a cube [0, 1]

3 m discretized with 420 points and 1402 tetrahedra. Forces of 2000N are applied to the bottom and upper faces while other faces are free. The Young modulus is E = 15M P a and Poisson ratio ν = 0.49.

At each simulation, regions are randomly generated by picking seeds among the tetrahedra, then extended according to a classic Voronoï sampling followed by a Lloyd relaxation.

The integration samples are generated with the same process within the regions. Their number is proportional to the volume of the zone. To evaluate the simulation results, we use a full order finite element model referred as FEM in which the control nodes are placed at the mesh vertices and the DOFs are the 3 components of the displacement vector. Barycentric shape functions are used within the elements and the integration is performed within the tetrahedra.

Table 1 shows the computation time increasing accordingly to the number of affine control nodes and integration points (from 0.088 to 2.669 seconds). Choosing suitable numbers of DOFs and integration points will consist in a trade off between speed and accuracy, depending on the complexity of the deformation to approximate. By selecting as many integration points as the number of tetrahedra, the simulation time of the proposed method is similar to the FEM's.

The thickness of the cube under compression is studied for one, five and nine control nodes. The relative error compared to FEM is listed in Table 1. The compression, consisting in a simple scaling, is exactly represented by one affine transformation. In this case, one integration sample is sufficient to represent the deformation. The results for one region validate the simulation in this sense but are not physically interesting. The other experiments show that a denser sampling will lead to a more accurate solution. This is expected since the deformations are captured in more local areas. Since the regions are randomly generated, their size and position can vary, leading to areas less represented in integration points. We observed that the stiffness is higher in the well sampled areas. Such a heterogeneity causes the upper and bottom faces to tilt. This phenomena can be observed for a high number of regions and low number of integration points, corresponding to a high uncertainty in Table 1. For example, the top left figure in Fig. 4 shows the random sampling of integration points. At the bottom of the cube, the lack in integration points introduces errors in the deformations leading to the bending of the cube's edge.

Inflation of a truncated ellipsoid

The second experiment deals with Dirichlet boundary conditions on an axially symmetric shape. We study the deformations of a truncated ellipsoid discretized with 1801 points and 6181 tetrahedra. The regions are linearly generated. One fixed-size region is placed at the apex. A pressure force of 1000N is applied on the interior triangles along the normal surface. The following results are obtained with 4000 integration points, E = 1GP a and ν = 0.49.

The top mesh vertices are fixed with a projective constraint. To do so, we modify the control node location and shape function of all regions sharing at least one mesh vertex with the boundary. For these regions, the affine control node is arbitrary placed on the boundary. The shape function is equal to 1 on the boundary vertices and then classically decreases on the remaining vertices. We finally project the constraint on the previously selected DOFs, thus constraining the mesh vertices on the boundary.

The orange curve on Fig. 5 represents the converged FEM deformation. Our method tends to converge toward the FEM one while increasing the number of control nodes. These results are two to three times faster than the full order method.

Discussion

On the ellipsoid case, we can notice that the external surface position is directly dependent on the number of regions. Considering the circular shape of the ellipsoid, tangential forces in large regions cancel each other resulting in a lower pressure force mapped to the control node. This effect tends to minimize the inflation for small number of regions.

Furthermore, it is good to notice that with the projective constraint, the slope of the curve near the Dirichlet boundary condition is well approximated. Another solution is to apply very stiff forces on the top nodes as boundary condition. However, this leads to stiffer top regions. In this case, the forces are mapped to the control nodes, thus influencing other vertices contained in these regions. This problem does not appear in the classical FEM since the shape functions are defined at the element level. With the projective constraint applied to the control nodes on the boundary and smoothly decreasing shape functions, this problem is tackled.

Finally, we observed small bumps on the deformed mesh, more visible for a large number of regions and probably due to the shape functions definition. The inflation seems more pronounced at the center of the regions than the reference one as it is visible on the ellipsoid Fig. 5. This phenomena is reduced when a, the parameter of the sigmoid shape function, is high but leads to a smaller deformation.

Conclusion

The method presented in this paper shows promising advances in reducing the computation time and number of integration samples. Deformations on simple geometries are well approached with suitable parameters and validate the method. More complex geometries, such as the truncated ellipsoid, which is a simplified model of a ventricle, are also well approached by poly-affine deformations providing a sufficient number of DOFs. To be efficient, this hybrid method still requires a fine mesh for the offline computation of the geodesic distance as well as the regions volume. With circular shapes, the required number of control nodes to obtain an accurate solution is still high to efficiently use the large regions defined by the AHA segments. Moreover, further work needs to be done to understand the bumps issue. Finally, simple boundary conditions are well handled providing an extra care on the control nodes location and the mapping between the mesh vertices and the basis of affine degrees of freedom.

Fig. 1 .

 1 Fig. 1. 29 AHA regions (left) and oriented control nodes (right)

Fig. 3 .

 3 Fig. 3. Example on a cardiac topology

Fig. 4 .

 4 Fig. 4. Thickness of the cube under compression. Top left figure corresponds to the deformed cube with 9 regions and 200 integration points. The yellow crosses represent the deformation gradients at these points. Frames indicate the position of the control nodes.

Fig. 5 .

 5 Fig. 5. Left: ellipsoid at rest (blue with edges) and deformed with full (orange wireframe) and reduced (plain brown) models. Right: position of the points on the external surface in the plane x + , 0, z . nz corresponds to the number of zones on the vertical axis and nt along the circumference of the ellipsoid (plus one zone for the apex). The number of DOFs of each model is indicated in parenthesis.

Table 1 .

 1 Cube computation time and relative error comparison. The error uncertainty represents the tilting of the upper and bottom faces. FEM computation time is 2.41s.

	Gauss Points / Zones	1	5	9
	200	Time (s) Error (%)	0.088 0.0	0.575 3.796 ± 9.864	1.157 4.511 ± 4.469
	500	Time (s) Error (%)	0.116 0.0	0.891 1.142 ± 1.297	1.823 2.192 ± 3.987
	800	Time (s) Error (%)	0.144 0.0	1.125 0.312 ± 1.228	2.263 0.646 ± 1.608
	1100	Time (s) Error (%)	0.173 0.0	1.409 0.104 ± 1.203	2.669 0.115 ± 1.064
	1400	Time (s) Error (%)	0.204 0.0	1.293 0.081 ± 0.272	2.402 0.081 ± 0.380


domain: math sde

text:
Janelle K K Hammond
email: janelle.hammond@inria.fr

Ruiwei Chen

Vivien Mallet

Meta-modeling of a simulation chain for urban air quality

Keywords: Model order reduction, Surrogate model, meta-model, Reduced basis methods, Gaussian dispersion model, Air quality simulation, Modeling chain

come

Introduction

Air quality simulations at urban scale are a key tool for the evaluation of population exposure to particulate matter and gaseous air pollutants. The simulations are however subject to costly computational requirements and complicated implementation. Studies in exposure estimation or uncertainty quantification, for example, require many solutions to the model. Advanced models can be rendered feasible in this context if we can reduce the computational cost without significant loss of accuracy.

Let us consider a generic stationary model over a physical domain Ω ⊂ R and parameter domain

D ⊂ R Np M : Ω × D → R N p → c(p)
The model output for a given parameter vector p ∈ D, c(p) ∈ R N , will be a large-dimension vector representing the solution over a grid coverin Ω. M can represent various types of atmospheric pollution models, from highly complex formulations based on partial differential equations and fluid dynamics [START_REF] Milliez | Computational Fluid Dynamical Modelling of Concentration Fluctuations in an Idealized Urban Area[END_REF][START_REF] Tominaga | CFD simulation of near-field pollutant dispersion in the urban environment: A review of current modeling techniques[END_REF] to simpler, and more commonly operational, formulations such as Gaussian dispersion models. Even in the case of the (comparatively) simpler models, the computational time necessary for the solution of M in practical applications over large domains with many parameters (e.g. emissions sources) can be high. This would make numerous solutions to the model too costly in practice. Methods of Model Order Reduction (MOR) can reduce computational costs without introducing significantly increased model error, and for a range of varying parameters p ∈ D.

Various MOR techniques have been studied in the context of air quality models (AQMs). In [START_REF] Mallet | Meta-modeling of ADMS-Urban by dimension reduction and emulation[END_REF] the metamodeling technique described in section 2 was tested on pollutant concentration fields over Clermont-Ferrand approximated by the ADMS-Urban model [START_REF] Carruthers | Development of Adms-Urban and Comparison with Data for Urban Areas in the UK[END_REF] using daily profiles for traffic emissions. In [START_REF] Lee | Emulation of a complex global aerosol model to quantify sensitivity to uncertain parameters[END_REF], statistical emulation was used to evaluate the sensitivity of some input parameters on a global aerosol model. A Gaussian process emulation was used for the study of model uncertainty in [START_REF] Armand | Probabilistic safety analysis for urgent situations following the accidental release of a pollutant in the atmosphere[END_REF] for accidental release scenarios. Gaussian process emulation was also used in [START_REF] Girard | Emulation and Sobol' sensitivity analysis of an atmospheric dispersion model applied to the Fukushima nuclear accident[END_REF] for the Sobol' sensitivity analysis of a dispersion model representing the Fukushima event.

In this paper, we will consider a modeling chain for air quality modeling over the agglomeration of Clermont-Ferrand and surrounding area in France. Air quality models are known to commit significant errors [START_REF] Shorshani | Modelling chain for the effect of road traffic on air and water quality: Techniques, current status and future prospects[END_REF][START_REF] Russell | NARSTO critical review of photochemical models and modeling[END_REF][START_REF] Milliez | Computational Fluid Dynamical Modelling of Concentration Fluctuations in an Idealized Urban Area[END_REF][START_REF] Zhang | Real-time air quality forecasting, part II: State of the science, current research needs, and future prospects[END_REF], however these errors are strongly dependent on the calibration and inputs to the model. Providing more precise input data, such as data on pollutant emissions from road traffic, can greatly improve the accuracy of the modeled concentration field. The advantage of a modeling chain is the use of the best (most precise) information available on various inputs by using traffic and emissions models. In [START_REF] Shorshani | Modelling chain for the effect of road traffic on air and water quality: Techniques, current status and future prospects[END_REF], the authors provide a review of modeling chain techniques for traffic pollutant emissions, atmospheric dispersion, and effects on water quality.

The modeling chain studied here consists of the dynamic traffic assignment model LADTA [START_REF] Leurent | On network assignment and demand-supply equilibrium : an analysis framework and a simple dynamic model[END_REF][START_REF] Leurent | Large problems of dynamic network assignment and traffic equilibrium : Computational principles and application to Paris road network[END_REF], an emissions model Pollemission [START_REF] Chen | Pollemission software computing traffic emissions of atmospheric pollutants with copertiv formulations[END_REF] based on COPERT-IV emissions database [START_REF] Ntziachristos | COPERT: A European Road Transport Emission Inventory Model[END_REF], and a Gaussian AQM, Sirane [START_REF] Soulhac | The model SIRANE for atmospheric urban pollutant dispersion; part I, presentation of the model[END_REF]. The computation of a pollutant concentration field over the agglomeration for any given time requires the solution of each model in the chain, which proves costly for long time periods.This brings us back to MOR techniques. However in this case, we have a chain ofmultiple models to reduce, which leads us to questions on the implementation of MOR techniques : could a single reduction over the full chain be feasible, or will satisfactory results require a chain of meta-models ? How can we treat the large parameter dimension of the chain ?

We resort to projection-based MOR techniques based on Reduced Basis (RB) [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] to construct cheap and accurate meta-models. A projection-based meta-model for the dynamic traffic model was built in [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF]. Here we will complete the model chain with the conversion between traffic assignment and emissions model outputs and pollutant dispersion model inputs, then construct a meta-model for the AQM to build a low-cost meta-model chain for the entire system. The motivation for this choice will be explained in section 3.

In section 2, we will describe the meta-modeling technique based on RB methods. In section 3, we will describe the case study over Clermont-Ferrand : input and measurement data, computational domain, and selected models.

In section 4, we will summarize the results of the meta-model on the AQM chain, studying accuracy, precision, and computational savings.

meta-modeling Methods

Computation times for large problems are commonly on the order of hours, making many-query contexts, such as sensitivity analysis and optimization, hardly feasible. Model reduction methods are of great interest to applications of parametrized problems involving many-query or real-time study. We will rely on a projection-based method of model order reduction in which the output solution space to the model

X = {c(p)|p ∈ D} ⊂ R N ,
where the parameter dimension is N p , is represented by a reduced basis of small dimension. While the model output is of high dimension N , the reduced order solution will be of dimension N N . We will begin here by detailing the MOR method as applied to the AQM part of the chain, and we will discuss the details of the full meta-model chain in section 3.

Reduced Basis Method

Let us consider a model, or model chain, M which takes input parameter vector p ∈ D ⊂ R Np and computes an output vector c(p) over a grid of N points. Our objective is to construct a reduced basis {Ψ AQ n } 1≤n≤N of N basis functions approximating the concentration solution space X such that the projection of any simulated state, Π N c(p), onto the reduced basis is sufficiently precise. The basis representing atmospheric concentration fields will be denoted by AQ (air quality). To construct a RB we first need to sample a large number of solutions in X . This so-called training set should represent the variability in the solution states. We will sample the solution space by Latin Hypercube Sampling (LHS). Then we will construct the RB by Principal Component Analysis (PCA).

We

c(p i ) -c - N n=1 Ψ AQ n Ψ AQ n T (c(p i ) -c) 2 2 = N k=N +1 λ k , (1)
for eigenvalues λ arranged in decreasing order. N principle component basis functions ψ AQ n are selected to represent I N = 98% of the variability in the concentration state. This means that the error of projecting any member of the training ensemble onto the basis, Err N , will be bounded by the tolerance

Err N ≤ N = √ 1 -I N .
For any new parameter, we can thus represent the solution as

c(p) Π N c(p) = c + N n=1 α AQ n Ψ AQ n (2) with projection coefficients α AQ n = Ψ AQ n T (c(p) -c).

Statistical Emulation

Once we have constructed the reduced basis by PCA, we need a reduced order modeling scheme to approximated new solutions. Classical reduced basis methods which replace the approximation space with the reduced basis space are intrusive and require the modification of the computational code. We would like to use a non-intrusive method which can be applied to a black-box model or model chain, which is particularly pertinent in the context of operational models. We consider meta-modeling by the emulation of projection coefficients α AQ n , 1 ≤ n ≤ N .

First we select a linear trend, which will be a least squares regression R n (p) = Np k=1 β n,k p k , calculated from the training simulations {c(p i )} 1≤i≤Ntrain . To this we add an interpolation term on the residuals α n (p i ) -R n (p i ),

I N (p) = Ntrain i=1 ω n,i φ d θ (p, p i ) .
We chose to compute this interpolation using radial basis functions (RBF). We chose cubic RBFs φ and a weighted Euclidean distance d θ (•, •) to represent the varying ranges of each input parameter.

d θ (p 1 , p 2 ) = Np i=1 θ i (p i 1 -p i 2 ) 2 , ( 3
)
max p∈D p i -min p∈D p i 2 .
We then define the emulated projection coefficients as follows.

αAQ n (p) = Np k=1 β n,k p k Least squares regression + Ntrain i=1 ω n,i φ d θ (p, p i ) Residual interpolation . ( 4
)
The weights {ω n,i } 1≤n≤N ;1≤i≤Ntrain are chose such that the interpolation is exact for the sample points

{p i } 1≤i≤Ntrain , αAQ n (p j ) = α AQ n (p j ) = Np k=1 β n,k p j,k + Ntrain i=1 ω n,i φ d θ (p j , p i ) . ( 5
)
The emulated solution is finally

ĉN (p) = c + N n=1 αAQ n Ψ AQ n . ( 6
)
The regression represents the relation between the model parameters and the RB projection coefficients, and computed from the training set (p i , α(p i )) 1≤i≤Ntrain . This provides an initial trend to be corrected by the interpolation. In practice, the interpolation of the residual is the most important part of the emulation. In [START_REF] Mallet | Meta-modeling of ADMS-Urban by dimension reduction and emulation[END_REF] this method of approximating projection coefficients is compared to approximation by Kriging. The two meta-models showed similar results, and we chose RBF emulation for its simpler (and thus more accessible in operational applications) implementation and lower computational cost.

Case Study on Clermont-Ferrand

In this work we will apply the meta-modeling method described in section 2 to a modeling chain over the city of Clermont-Ferrand in France. We will build a meta-model chain representing road traffic emissions and the dispersion and reaction of pollutants over the urban agglomeration and surrounding area using data over a two-year period form 2013-2015. The model chain is represented in figure 1.

Traffic Emissions Modeling

Traffic emissions modeling is done using the dynamic traffic assignment model LADTA. A meta-model was constructed [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF] to represent the traffic flow and speed simulations over a road network of 19, 628 oriented links, where nearly 45, 000 traffic flow observations are available each day. Emissions of NOx and PM are computed using Pollemission code [START_REF] Chen | Pollemission software computing traffic emissions of atmospheric pollutants with COPERT-IV formulations[END_REF] based on the COPERT-IV emissions database [START_REF] Gkatzoflias | COPERT 4 : Computer programme to calculate emissions from road transport[END_REF][START_REF] Eea | EMEP/EEA air pollutant emission inventory guidebook -Part B.1.A.3.b.iiv Road transport[END_REF]. A detailed description of this section of the modeling chain and its input parameters can be found in [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF]. The varying input parameters consist of 23 traffic parameters and 6 emissions parameters. These parameters are time-dependent or considered sources of uncertainty. They include temporal traffic demand, computed using traffic observations, the capacity and speed limits of traffic network links, multiplicative coefficients on origin-destination matrices representing spatial trends of traffic, traffic direction (morning versus evening), engine size, bype, and emission standards of the vehicle fleet, and ratio of heavy-duty vehicles to personal cars.

The emissions model provides traffic emissions estimations for NOx and PM . In the deterministic case, we set the ratio N O2 N Ox = 0.15 [START_REF] Carslaw | Evidence of an increasing NO2/NOX emissions ratio from road traffic emissions[END_REF][START_REF] Beevers | Trends in NOx and NO2 emissions from road traffic in Great Britain[END_REF][START_REF] Kurtenbach | Primary NO2 emissions and their impact on air quality in traffic environments in Germany[END_REF], and the ratio P M2.5 P M10 = 0.75 [START_REF] Gillies | On-Road Particulate Matter (PM2.5 and PM10) Emissions in the Sepulveda Tunnel[END_REF][START_REF] Querol | Speciation and origin of PM10 and PM2.5 in selected European cities[END_REF]. In order to construct a meta-model which can account for varied or uncertain speciation ratios, we will draw LHS parameters for the training ensemble in the intervals

(p N O2 , p P M2.5 ) ∈ [0.1, 0.25] × [0.65, 0.8].
The output of the traffic-emissions coupling is the emissions on each link of the traffic network in g/15min.

Air Quality Modeling

Air quality modeling is done using the Gaussian dispersion and reaction model Sirane [START_REF] Soulhac | The model SIRANE for atmospheric urban pollutant dispersion; part I, presentation of the model[END_REF][START_REF] Soulhac | The model SIRANE for atmospheric urban pollutant dispersion; PART II, validation of the model on a real case study[END_REF]

Modeling Chain

The modeling chain consists of these three steps -traffic modeling, emissions calculation, and dispersionreaction modeling -and the conversions between outputs and inputs. In figure 2

Surrogate Modeling Chain Construction

As noted above, the traffic emissions on a geographically finer road network provided as input to the air quality model represent over 47, 000 line sources. In the context of model order reduction, this represents as many parameters, which in the practice of projection-based reduction methods makes the identification of the projection coefficients α AQ n dependent on 47, 000 parameters unfeasible (or impossible). We thus need to reduce the complexity of the problem bu reducing the dimension of the input parameters. To do so we will construct a reduced basis of the traffic emissions, again using PCA.

Reduction of line emissions.

We currently have the full chain parameter vector p T f ull = (p T traf f ic , p T e , p T AQ ), where the outputs of the emissions model consist in a coefficient for each of the links in the road network. These coefficients are then treated as the (very large) input parameter vector for the air quality model. To reduce the dimension of this vector, we will use the same method as in section 2. E(p i traf f ic , p i e ). We can compute the orthogonal projection of any emissions field onto the traffic emissions RB as follows.

E(p traf f ic , p e ) Π N lin E(p traf f ic , p e ) = Ē + N lin n=1 (E(p) -Ē) T Ψ E n Ψ E n = Ē + N lin n=1 α lin n Ψ E n . ( 7
)
For our case study, we chose N lin = 11 to represent 95% of the variability of the emissions solutions. This corresponds to a relative projection error tolerance over the training samples of 2 lin = 0.05. In the model chain, the over 47, 000 line source parameters will henceforth be replaced by the N lin = 11 projection coefficients {α lin n } n≤N lin , and the traffic emissions field for a given parameter approximated by its projection Π N lin E(p traf f ic , p e ) onto the traffic emissions RB. We perform the same reduction over the hourly surface emissions with N surf = 1 and projection coefficient α surf . In figure 3, we can see the largest singular values of the PCA step, and the relative mean projection errors of the training traffic emissions simulations onto the RB {Ψ E n } 1≤n≤N lin , as defined by In figure 4 we can see the first 4 principal components of the traffic emissions RB.

Err N = 1 N t Ntrain i=1 Π N E(p i ) -E(p i ) 2 E(p i ) 2 . ( 8
)
Construction of the air quality meta-model.

We now can write the reduced concentration model parameters p T c = (α T lin , α T surf , p T AQ ). We will construct a metamodel of the air quality model to complete the meta-modeling chain, with reduced full parameters as described in table 1. The choice to build a separate air quality meta-model to complete the chain of meta-models (as opposed to a meta-model of the chain) was for multiple reasons. First, to allow multi-level assessment using traffic flow and air quality measurement data (a possibility particularly pertinent in a study of uncertainty quantification). In figure 5 we compare the parameters α lin n selected by these two methods by plotting the parameter spaces In table 2, we set the ranges of each input parameter which defines the parameter space D.

(
From a training set of N train = 9347, we compute the NO 2 concentration fields c(p f ull ) to construct a reduced basis {Ψ AQ n } 1≤n≤N by PCA, representing the variations of the concentration fields centered around the sample concentration mean c. We set the RB dimension N = 5 to represent 98% of this variability. In figure 6, we can see singular values of the matrix of centered solutions ȲAQ defined in section 2, and the projection errors of the training ensemble of the atmospheric pollution model, as defined by equation [START_REF] Shorshani | Modelling chain for the effect of road traffic on air and water quality: Techniques, current status and future prospects[END_REF].

In figure 7 we see the first 4 principal components of the concentration RB. We can see that the first basis function represents urban background concentration in the denser urban areas. The second seems to represent additional pollution from traffic. The third appears to represent strong wind from the East, while the fourth displays wind from the North.

For any new parameter value, the concentration field can be approximated by the orthogonal projection onto the RB, for projection coefficients

{α AQ n } 1≤n≤N , c(p f ull ) Π N c(p f ull ) = c + N n=1 α AQ n Ψ AQ n . ( 9
)
Finally we use the statistical emulation method described in section 2.2 to construct an emulator of the concentration projection coefficients α AQ n . The full chain can be computed with a single code which applies the traffic-emissions meta-model, the calculation of emissions RB projection coefficients, and the atmospheric pollutant meta-model. This meta-model chain provides outputs on traffic flow, speed, and traffic emissions over the road network, and NO 2 concentrations over a 20m-resolution grid. 12

Results

In this section, we will summarize the results of the method described in section 2 to the case study in section 3 using data over the month of November 2014. Traffic flow measurement data serves as inputs to the model chain for deterministic simulation, and data on pollutant concentration serves to study model and meta-model performance. We will compare the meta-model output to simulations from the full model Sirane, as well as to concentration observation data, and we will assess computational savings.

Meta-model performance

We introduce the following statistical scores commonly used for evaluation of models [START_REF] Mallet | Meta-modeling of ADMS-Urban by dimension reduction and emulation[END_REF] : the normalized mean square error (NMSE), the normalized root mean square error (NRMSE), and the correlation. We define here the output functionals m : Ω → R associated to each of the concentration sensors m, such that the observation data

y obs m (p) = m (c true (t)
) 1≤m≤M . RMSE = 1 M M m=1 (c m -y obs m ) 2 . ( 10
) √ NMSE = 1 M M m=1 (c m -y obs m ) 2 cȳ obs . ( 11
) Correlation = M m=1 (c m -c)(y obs m -ȳobs ) M m=1 (c m -c) 2 M m=1 (y obs m -ȳobs ) 2 . ( 12
) Bias = 1 M M m=1 (y obs m -c m ) (13)
Finally we define the NRMSE as RMSE ȳobs , and the MNRMSE as the mean over all sensors (or grid points) of the NRMSE calculated over the concentration c i at each sensor (or grid point) over the month.

MNRMSE = 1 N grid N grid i=1 N RM SE(c i ) (14)

Comparison with the full model chain

We first analyze the precision of the meta-modeled concentration fields as compared to the full model Sirane.

This will help us understand the ability of the meta-model to reproduce the concentration state and quantify the loss of precision caused the the dimensional reduction.

In figure 8, we see statistical scores spatially mapped over the meta-model domain. The NRMSE shows that the emulated solutions perform well in approximating the urban background concentration levels, but don't capture the highest concentrations along the large highways as well. The correlation map also shows low correlation between the meta-model and full model only along the roadways, where the dimensional reduction has failed to capture the extent of the increased concentrations due to traffic emissions. Finally the bias map shows that the meta-model generally predicts higher concentrations in the denser urban areas when compared to the full model, again matching the trend of the dimensional reduction reducing the sensitivity of the meta-model to sharp spatial variations in concentrations. However, the areas with poor scores remain limited, and we must also consider the significant error that will inevitably be committed by the full model in the next section. In figure 9 We see the relative errors of the full model concentration projected onto the reduced basis {ψ AQ n } 1≤≤N , averaged over the set of deterministic simulations for the month of November 2014. We also see the emulated concentration relative error, averaged over the same set of simulations. While the emulation of the projection coefficients is globally responsible for the most error, we can see that the regions with the highest projection error correspond to high errors in the meta-model as well. This is expected, as the emulated solution can only perform as well as the projected solution. We see that larger errors are located on roads, mostly the large highway and outside the dense urban area. Meta-model error remains below 20% over a large portion of the domain, which shows that much of the spatial variation of the concentration is captured by the reduced order solution.

In table 3 we can see statistical scores comparing the meta-modeled concentration to the full concentration model over all hours of November 2014. We compare both the entire grid (here c m is the concentration at a grid point and M = N grid is the total number of grid points) and at the NO 2 sensor locations. While the dimensional reduction means the meta-model does not fully capture spatial variations of the simulated concentration state, we can see that the relative RMSE errors are satisfactorily low, and the correlation between the two is very high. In figure 10

State Estimation

Comparison with observational data

We next analyze the accuracy of the full model and meta-model compared to observational data on NO 2 concentrations. In figure 11, we see the temporal profile of average NO 2 concentrations at M = 4 sensor locations :

1 M M m=1 c m . We compare observed, emulated, projected and Sirane modeled concentrations of all weekdays in November 2014. We see that the bias in the modeled concentrations underestimating peak concentrations, notably during heavy traffic periods in the mornings and evenings. We also notice a seemingly delayed reaction of the model chain to the pollution increase during the evening peak hour. In [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF] this delay was less evident, suggesting that factors such as the dispersion and reaction parametrizations in the AQ model or the averaging of time scales from 15 minutes to one hour may have an effect. The exploration of this question will require more study of uncertainties in the model chain. We notice the the temporal trend representing morning and evening peak hours in traffic is reproduced by the model chain. We also note that the emulated concentrations are closer to the observations than the full model. This is likely due to the "smoothing" effect of the dimensional reduction causing less sharp concentration variations, as small parts of the modeled concentration fields are not reproduced by the reduced basis.

In table 4 we compute statistical scores over the month of November 2014, comparing the full model simulations Finally, the station Chamalières is located outside the city center, where the model exhibits a higher level of bias.

The performance of the meta-model with respect to observation data is highly satisfactory.

.

In figure 12 While we have seen that the model reduction by statistical emulation causes loss of precision, and the metamodel simulations contain error with respect to the full model, comparing to observation data suggests that this error is not significant with respect to the model error inherent to operational models for urban air quality, and does not reduce the accuracy of the predicted concentrations at sensor locations.

Computational savings

We have seen that the meta-model chain produces satisfactory results when compared to observational data, and determined that the loss of precision due to the dimensional reduction is not higher than the error committed by the full model. Now we will show the computational savings afforded by the meta-model chain. In table 5 we can see the computational times required for a single simulation of the chain by the meta-models or the full models.

The meta-models depend on three reduced basis, representing traffic for the traffic assignment meta-model, road emissions for the reduction of pollution model input dimension, and concentration fields for the pollution metamodel. The initialization of the meta-model chain requires loading these basis and building the RBF emulators.

Once the chain is initialized, it can be run for any number of simulations at very low cost, under 0.1 seconds for a simulation representing a one-hour period. In comparison, the full model chain requires nearly three hours for a single simulation. The offline construction of the meta-models required 6000 traffic model simulations [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF] and 10000 pollution model simulations, which represents a significant computational investment. However, these meta-models are trained over training points {p i } 1≤i≤Ntrain ∈ D representing two years of data, and once constructed are useful for study over multiple years. In the absence of high performance computing machines or clusters, the simulations can be run using a pseudo-parallel technique running one simulation per core on desktop calculation machines.

The Sirane simulations described in section 3 took around one day using this method on multiple machines of 64GB RAM or less. Once the meta-model chain is constructed, the online phase for the simulation given any parameter p ∈ D is very cheap, which makes real-time or many-query contexts possible, for example for use in uncertainty quantification study.

Conclusions

In this work we constructed a meta-model chain by statistical emulation of reduced basis projection coefficients for an urban air quality modeling chain over the agglomeration of Clermont-Ferrand. We used the road traffic meta-model constructed in [START_REF] Chen | Metamodeling of a dynamic traffic assignment model at metropolitan scale 43[END_REF], built a reduced basis representing road traffic emissions, and constructed a second meta-model of NO 2 concentration fields over the agglomeration, substituting thus a low-cost chain of meta-models for a computationally costly modeling chain over a large urban area. This required the dimensional reduction of the inputs to the atmospheric pollution model, and the appropriate sampling of the parameter spaces to construct a satisfactory reduced basis and reduced order modeling scheme. We reduced computation time from over two computational hours per simulation representing an hourly concentration field to under 0.1 second. Results show good precision of the meta-model simulations with respect to the full model chain, and similar accuracy when compared to measurement data. The meta-model can be used in various applications requiring numerous solutions to the model chain, rendering studies requiring many solutions to the model chain computationally feasible. In future work, we will use this low-cost modeling chain in the study of uncertainty quantification and the propagation of uncertainties throughout the meta-model chain.

Figure 1 :

 1 Figure 1: meta-modeling chain over Clermont-Ferrand.

  we can see the traffic flow (veh/h/link) and associated emissions (g/km/s), and NO 2 concentration ( µg m 3 ) simulations at 8a.m. on a Tuesday in November 2014, provided by the traffic meta-model and full air quality model. The task remains to reduce the computational time required to obtain concentration fields by constructing a meta-model for the entire chain.

Figure 2 :

 2 Figure 2: Simulations over Clermont-Ferrand on 18/11/2014 at 8am. Traffic flow (left), NO 2 emissions (center), NO 2 concentration (right).

Figure 3 :

 3 Figure 3: Left : Singular values of the emissions mass matrix. Right : L 2 relative mean projection errors of the LHS training ensemble of road traffic emissions fields onto the RB.

Figure 4 :

 4 Figure 4: Principal components of the emissions mass matrix, first four.

Figure 5 :

 5 Figure 5: Projection coefficients on the traffic emissions basis from LHS performed directly on [α lin min , α lin max ] N lin (blue) compared to the projection coefficients of traffic emissions model outputs E(p traf f ic , pe) (red) over a training ensemble of parameters (p i traf f ic , p i e ) selected by LHS. Left : the parameter space of (α lin 1 , α lin 2 ). Right : the parameter space of (α lin 1 , α lin 4 ).

Figure 6 :

 6 Figure 6: Left : Singular values of the NO 2 concentration field mass matrix. Right : L 2 relative mean projection errors of the LHS training ensemble of NO 2 concentration fields onto the RB.

Figure 7 :

 7 Figure 7: Principal components of the NO 2 concentration field mass matrix.

Figure 8 :

 8 Figure 8: Left : NRMSE (10) of the emulated NO 2 concentration field compared to the projected solution, for parameters over the month of November 2014. Center : correlation (12) over the same set. Right : Normalized bias (13) over the same set.

Figure 9 :

 9 Figure 9: Left : Relative mean error of the projected NO 2 concentration field compared to the full model solution, for parameters over the month of November 2014. Right : Relative mean error of the emulated NO 2 concentration field compared to the full model solution over the same parameter set.

  we see a visual representation of hourly scores of the meta-model solution compared to the full solution at each grid point for simulations corresponding to the month of November 2014. The NMSE (11) remains globally below 0.4, and the RMSE (10) often below 10 µg m 3 . Correlations scores are grouped above 0.75, and the bias distribution is nearly centered around -2 µg m 3 , showing a slightly higher concentration approximation by the meta-model, when averaged over the grid.

Figure 10 :

 10 Figure 10: Scores of the meta-modeled NO 2 concentration field compared to the full model solution, for parameters over the month of November 2014. Top left : NMSE (11). Top right : RMSE (10). Bottom left : correlation (12). Bottom right : bias (13).

Figure 11 :

 11 Figure 11: Mean NO 2 concentrations at 5 sensor locations over weekdays in November 2014. Curves show observations, full model simulations, projected simulations onto the reduced basis, and emulated solutions.

  we see a visual representation of daily scores of the meta-model solution and the full solution compared to NO 2 observations over the month of November 2014. The meta-model shows similar score distributions to the full model, excepting the occasional outlier. The RMSE (10) is below 25 µg m 3 on the majority of days for both the full and reduced simulations, with the bias distribution nearly centered around 10 -15 µg m 3 , showing an underestimation of concentrations by the simulations.

Figure 12 :

 12 Figure 12: Scores of the Sirane NO 2 concentration field compared to the observation data over the month of November 2014. Top left : NMSE (11). Top right : RMSE (10). Bottom left : correlation (12). Bottom right : bias (13).

  use LHS to select N train sample points (p 1 , . . . , p Ntrain ) in the parameter domain D, and compute model simulations from each point to build the training ensemble Y AQ = [c(p 1 ), . . . , c(p Ntrain )] to train the model reduction. As is common practice in PCA applications, we will first compute the ensemble mean c = 1

	Ntrain	Ntrain i=1	c(p i )
	of the training ensemble. PCA is computed on the centered ensemble ȲAQ = [c(p 1 ) -c, . . . , c(p Ntrain ) -c]. The
	eigenvalues {λ k } 1≤k≤N and eigenvectors {Ψ AQ k } 1≤k≤N of the covariance matrix CAQ = ȲAQ T ȲAQ of the
	training ensemble are such that
	Ntrain
	i=1

  indirectly in the formation of NO 2 . The background concentrations represent the imported concentrations of pollutants, that is, concentrations transported from distant locations to the city, and possibly dispersed from previous emissions in the case of stationary solution. We will provide inputs on NO 2 , NO, O 3 , PM 2.5 , and PM 10 . Input data on meteorological conditions and surface emissions sources are also provided. The AQM output is the NO 2 concentration over a grid at ground level, at 20m resolution. Hourly concentration observations are available over two years at 5 stations, or around 90, 000 NO 2 observations for analysis of model simulation outputs.

over a simulation domain of 180 km 2 . Sirane is used as a static model which approximates the solution at a given time of the transport-reaction equations satisfied by the pollutant concentrations. The traffic emissions over a relatively coarse road network are converted to g/s/link on a finer network representing over 47, 000 line sources. For the calculation of NO 2 concentrations, we provide the so-called background concentrations of pollutant species involved directly or

  We compute the emissions solutions E(p traf f ic , p e ) to construct a reduced basis {Ψ E n } 1≤n≤N lin by PCA, representing the variations of the emissions fields centered around Ē = 1

		Ntrain
	Ntrain	i=1

1. We first select a set of training parameters (p T traf f ic , p T e ) by LHS to represent the variations of these parameters in the admissible parameter space D.

Table 1 :

 1 Summary of input parameters to the full meta-model chain.

	Dynamic Traffic Model	Emissions Database	Air Quality Model
	LADTA	COPERT IV	Sirane

α lin 1 , α

lin

2 ) and (α lin 1 , α lin 4 ). We can see that the parameter spaces in red, which correspond to performing LHS on

p traf f ic ∈ R 23 p emiss ∈ R 6 p AQ ∈ R 22

Temporal traffic demand ; link capacity and speed limit ; distance traveled data ; direction coefficient Vehicle fleet data Meteorological data ; traffic, surface and railway emissions ; background pollution

Table 2 :

 2 Model chain input parameter ranges.

Table 3 :

 3 Statistical scores of the meta-model approximation results compared to the scores of the model chain using the full air quality model.

Table 4 :

 4 Statistical scores of the meta-model approximation results compared to observation data. NRMSE[START_REF] Zhang | Real-time air quality forecasting, part II: State of the science, current research needs, and future prospects[END_REF],

	Statistical scores over		Sirane solution c(p)		Meta-model solution ĉ(p)
	one-month simulations	NRMSE	√	NMSE Bias Correlation NRMSE	√	NMSE Bias Correlation
	All 4 stations	0.479	0.499	12.88	0.746	0.461	0.467	10.28	0.728
	Lecoq	0.368	0.351	6.81	0.847	0.364	0.347	6.5	0.846
	Montferrand	0.366		0.33	2.04	0.833	0.377	0.325	-1.6	0.836
	Gare	0.655	0.587	33.13	0.776	0.469		0.52	20.43	0.783
	Chamalières	0.474	0.573	16.05	0.719	0.498		0.55	15.8	0.648
										√	NMSE (11),
	Correlation

Table 5 :

 5 Computation times using the meta-model or full model chain.

	Computational times
	Meta-model simulation	Traffic Emissions	Compute αlin (p m )	Compute αAQ (p m )	Total CPU
	Initialize meta-model chain	-	-	-	24 min
	Emulating α AQ (p m )	0.05 sec	0.006 sec	0.02 sec	0.076 sec
	Full model simulation	Traffic Emissions	-	Compute c(p m )	Total CPU
	Simulation of c(p m )	117 min	-	23 min	140 min

Acknowledgments

This research is supported by the French National Research Agency (ANR, the Agence Nationale de la Recherche), project ANR-13-MONU-0001, ESTIMAIR. The City of Clermont-Ferrand provided traffic flow measurement data. SMTC (Syndicat Mixte des Transports en Commun de l'agglomération clermontoise) provided the traffic network geometry for the agglomeration of Clermont-Ferrand and the static O-D matrix representing spatial traffic demand in the traffic assignment model. The SME NUMTECH provided data necessary for the model SIRANE, including geometrical features of the traffic network, meteorological data, emissions data, and background concentrations. We particularly thank David Poulet for his contributions to the data treatment. The SIRANE model was developed by "Laboratoire de Mécanique des Fluides et d'Acoustique", UMR CNRS 5509, University of Lyon, Ecole Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon I. We thank Lionel Soulhac for providing it and for his guidance. Atmo Auvergne-Rhone-Alpes (https ://www.atmo-auvergnerhonealpes.fr/) provided pollutant concentration measurements used in the calculation of statistical scores.


Mario Sánchez Aguilar

Ana Kuzle
email: kuzle@uni-potsdam.de

Kjersti Waege
email: kjersti.wage@matematikksenteret.no

Morten Misfeldt
email: misfeldt@learning.aau.dk

Mario Sánchez Aguilar

Introduction to the papers of TWG23: Implementation of research findings in mathematics education

Keywords: Implementation research, research findings, replication studies, innovation, change

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction to the papers of TWG23: Implementation of research findings in mathematics education

Origin of the Thematic Working Group 23

During five decades, the field of mathematics education research has generated a multitude of products, such as theoretical frameworks, concepts, didactic designs, solid findings, etc. Although the research community has always been concerned with the theory-practice relationship, it remains an open and a challenging problem how such products could be used and applied in practice. Many mathematics education researchers work on (large) developmental projects that rely heavily on previously documented research results. However, reporting on these projects in mathematics education scientific outlets have proven to be challenging, since such projects do not necessarily fall under the usual paradigms of research in mathematics education. Before CERME 10, this issue was identified. More concretely, there was no forum to discuss important issues related to implementation of research findings. For that reason, TWG23 was founded with the idea of being a forum dedicated to presenting and discussing empirical and theoretical studies focused on elucidating the enablers and general conditions that favor or inhibit the implementation of research products generated in our field in practice as well as reporting on research-based designs themselves. The original concern of the group was "How can we apply and implement stable research findings in «real life»?" and "How can we bring the accumulated research knowledge into practice?"

These questions remained to be in the focus of TWG23 also at CERME 11. We were embracing for papers and poster proposals addressing, from both empirical and theoretical perspectives, issues related to the implementation of research findings. In particular, we wanted to focus on a wide variety of "good examples" of implementation of research findings and products into practice, aiming at improving the teaching and learning of mathematics at all educational levels. Moreover, we were interested in empirical research and theoretical discussions that address the challenges and possibilities of "implementation research", as well as in-depth literature reviews that provide general and updated views of the state of development of this type of research in the field of mathematics education. Thus, the work of the group was organized around the idea of closing the identified "gap" and acting as a bridge between research and practice by focusing on multiple perspectives of implementation research.

Evolution of the TWG 23

At CERME 10, the working group undertook initial attempts to make sense of the construct of implementation. In the call for papers for the TWG23 at CERME 10, the construct "implementation research" was operationalized rather broadly, as a wide range of different kinds of didactical design, from task design, (model) lesson design, teaching modules and courses to design of entire programs at all educational levels. Furthermore, "implementation research" was inclusively treated as research on aspects of developmental projects, intervention projects, as well as research on aspects of the development and use of educational media, such as textbooks, apps, software and learning platforms. However, the call for papers required to make explicit connections between the reported designs and findings or insights from mathematics education research. As a result, many examples of implementation of research findings in mathematics education were presented, including process models, determinant frameworks, classic theories, implementation theories and evaluation frameworks (Jankvist, Aguilar, Bergman Ärlebäck, & Waege, 2017)

Understanding our emerging "paradigm" and its blind-spots and tacit assumptions

It was apparent at CERME 11 that many of the contributions shared certain ideas of what implementation in educational situations is. The ideas were shared both in terms of terminology and in terms of key references. Two key concepts apart from implementation are worth mentioning, namely innovation as described by [START_REF] Rogers | Diffusion of innovations[END_REF] and his famous work on the diffusion of innovations, and central factors of "implementation" as proposed by [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF]. This means that many of us considered "implementation" as innovation that causes a change when enacted in ordinary practice, and that the change is stipulated by relationships between the nature of the innovation, influential factors and key stakeholders.

Furthermore, we internalized the [START_REF] Rogers | Diffusion of innovations[END_REF] idea that innovations spread through different user segments, from innovators that participate in developing and realizing the innovation, over early adopters to a majority who instantly or relatively soon take up the innovation, and then to the late majoritytaking the innovation up lateand laggards who might never change practice.

The five factors that influence innovation proposed by [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF], were actively used in a number of the contributions in the group, as a lens to study the implementation and implementability of initiatives. These five factors are: characteristics of the individual users, organizational and environmental factors, attributes of the innovation, implementation support strategies, and implementation over time.

Apart from referring to the Rogers' and Century and Cassata's ( 2016) "standard models" for implementation (as something that causes a purposeful change, related to factors and stakeholders) in the group discussions, we also spent some time investigating the downsides of looking at changes in the educational sector through this lens.

First of all, we observed that the "standard models" tend to fetishize the "manifestation" of a change as an innovation. Change in the educational sector is typically not caused by one "thing". Rather, it is better studied as an interplay between many factors and forces. The constant focus on "innovation" in our work and discussions might have a problematic downside here. Another problem with the term "implementation as enactment of innovation" is that it overplays the intentional aspects of educational change. Many changes come as organic developments rather than as purposeful implementations.

We also began understanding the implementation paradigm by drawing a parallel to [START_REF] Stein | How curriculum influences student learning[END_REF] phases of curriculum, which was brought to the discussion by Boris Koichu. In [START_REF] Stein | How curriculum influences student learning[END_REF] model, the curriculum can be seen as consisting of different categories, namely written, intended, enacted, and attained curriculum based on what phase of the teaching and learning process one is focusing on. If we replace the written curriculum with (stable) research findings, a new model emerges (see Figure 1). Thus, in this new model published research refers to the research findings available to the teachers. Intended implications refer to the set of objectives to be accomplished in practice on the basis of research findings. Enacted implications refer to various learning activities or experiences of the learners on the basis of research findings in order to achieve the intended implications. Lastly, attained implications refer to the implementation outcomes with respect to student learning on the basis of the intended and enacted implications. Twelve papers and two posters were presented at CERME 11, all of them addressing a wide variety of topics and using different approaches. In this section we use the five categories proposed by [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF] in order to present an overview of the works discussed in this thematic working group. Although not mutually exclusive, to some extent these categories allow us to group the works presented in the TWG23 according to their interests and perspectives.

Inform innovation design and development

The research study by Ioannis Papadopoulos and Nafsika Patsiala falls into this category. They describe a pilot study, aiming at capturing the landscape of problem posing in a Greek grade four classroom. It is intended that the information produced by this pilot study informs and helps to setup a year-long intervention in the classroom aiming to develop the problem posing abilities of the students.

Understand whether (and to what extent) the innovation achieves desired outcomes for the target population

Two papers belong to this category. The first work, by Inga Gebel and Ana Kuzle, presents a problem-solving innovation for students in grades 4-6 has been designed and evaluated. During the evaluation phase, the researchers analyze students' problem-solving solutions with respect to their fluency and flexibility, and present the (dis-)advantages of the innovation in regular mathematics lessons.

The second work is by Helena Gil Guerreiro, Cristina Morais, Lurdes Serrazina, and João Pedro da Ponte. These authors try to understand, in a context of teachers' collaborative group, how emphasizing multiple representations can contribute to the learning of the rational numbers by elementary school students.

Understand relationships between influential factors, innovation enactment, and outcomes

Three of the papers presented in the TWG23 focus on identifying factors that influence innovation enactment. One of them was developed by Rikke Maagaard Gregersen, Sine Duedahl Lauridsen, and Uffe Thomas Jankvist, who focus on identifying the enablers and barriers for the implementation of the so-called Swedish Boost for Mathematics, which is one of the largest initiatives on improving mathematics teaching and learning in the Nordic countries in recent times.

The paper presented by Johan Prytz also analyzes the implementation of the Boost for Mathematics initiative, but this paper adopts a historical and comparative approaches. Through this perspective, the author shows the role that research and models of governance played in the New Math project in the 1960s and 1970s, and later in the Boost for Mathematics project.

Finally, the paper by Dorte Moeskaer Larsen, Mette Hjelmborg, Bent Lindhart, Jonas Dreyøe, Claus Michelsen, and Morten Misfeldt describes a recent attempt to implement on a large-scale an inquiry-based mathematics teaching in Danish compulsory school. In their analysis, they identify critical factors for this large-scale implementation.

Improve innovation design, use, and support in practice settings

Most of the works presented fall into this category. For instance, Tomas Højgaard and Jan Sølberg refer to a longitudinal project called KOMPIS (a Danish acronym that stands for "Competency Goals in Practice"). In particular, they present a two-dimensional content model derived from the KOMPIS project that aims at supporting competence-based curriculum development and teacher planning.

In turn, Morten Elkjaer discusses the design of a dynamic online diagnostic tool aimed at assessing students' mathematical misconceptions in lower secondary school. The design of such online tool is based on the application of research findings on mathematical misconceptions in algebra and numeracy when working with equations.

The implementation of notions from the inquiry-based mathematics teaching perspective (IBMT) is the focus of the work presented by Per Øystein Haavold and Morten Blomhøj. In particular, they discuss how the design of a four-year professional development project called SUM (A Norwegian acronym that stands for "Coherence through inquiry based mathematics teaching"), can support the implementation of research findings related to the IBMT approach.

The large-scale implementation of alternative models for multiplication is the focus the paper by Anna Ida Säfström, Ola Helenius, and Linda Marie Ahl. They make use of Vergnaud's theory of conceptual fields to produce a teaching design where models in the form of iconic representations serve as a means for creating patterns that make multiplicative invariants and structures visible.

The work by Nina Ullsten Granlund reports on a professional development project where preschool teachers are offered theoretical tools -in particular Bishop's theory of mathematical activityso that they can think about their teaching in a structured and explicit way, but respecting the playbased tradition on which Swedish preschool education is based.

The last article located in this category is the one presented by Boris Koichu and Alon Pinto. They report on the TRAIL project (Teacher-Researcher Alliance for Investigating Learning), which is a co-learning project between mathematics teachers and mathematics education researchers. They illustrate how this project favors teachers' adaptation of research procedures and ideas in their classrooms as part of participation in community educational research.

Develop theory

The main emphasis of the works included in this category is in the conceptualization of the implementation of research findings. This focus is clearly reflected in the work of Uffe Thomas Jankvist, Mario Sánchez Aguilar, Jonas Dreyøe, and Morten Misfeldt. Taking as a reference implementation research frameworks from outside the field of mathematics education, they try to outline what an implementation research framework in mathematics education could encompass.

Finally, Andreas Lindenskov Tamborg argues in his paper that in order to synthesize the research results in implementation research as an independent sub-field in mathematics education research, there is a need for a consistent vocabulary. He then proposes to combine [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF] definition of implementation research with theoretical notions developed in the realm of the documentational approach to didactics.

Results of the thematic group discussions

During our sessions, we organized two thematic discussions. The first thematic discussion focused on the topic of "Methodology in implementation research", while the second one concentrated on "Replication studies". Below we mention the main ideas that were addressed during those group discussions.

Methodology in implementation research

As of 2017, we as a community of mathematics educators interested in implementation research, got a forum to talk about and publish our "implementation" research, while analyzing our data from multiple perspectives. Here, different dependent variables can be considered: climate variables, learning variables, system variables or independent variables. Yet, in the wider community this type of research does not necessarily fall under the usual paradigm of "research in mathematics education". As a consequence, it is not highly-recognized and often is not considered for publishing in high ranking peer-reviewed articles, due to different factors. That said, the following question naturally arises: "What would count as a high-quality implementation research (recognized by the community and having impact in high-ranked journals)?" In order to answer this question, we had small group discussions organized around the following questions:

1. What perspective of implementation research does each project follow? Ones identified, what is the purpose of one's perspective? 2. How is innovation within each project measured and analyzed? What are methodological challenges faced by each "implementation" project? 3. Which of perspectives and purposes from (1) are closely aligned with the usual paradigm of "research in mathematics education"? 4. Which of these perspectives would count as high-quality implementation research? If so, why? If not, why not? 5. What implications would the decision from (4) have with respect to our methodology (e.g., instruments, units of analysis, such as products, processes, evaluation of particular factors, analysis)? 6. What challenges do arise in implementation research (e.g., related to data analysis, to ethical issues) and how could one go about them?

The discussions have shown that TWG23 is at its early stage where many fundamental issues are not yet completely clear. Both thematic discussions as well as the paper discussions made us enter a pathway towards formulating a more precise collective definition of the construct "implementation research in mathematics education" (see last section of this paper). Another issue identified was the complexity of implementation research, which consequently implies the variety of theories needed in order to inform it. Given the length limitation of publications in most of the high-ranked journals, it makes is extremely difficult to report on implementation research in its full complexityalthough this is also the case for other types of research. Lastly, in our community implementation of research findings is still not highly recognized by many colleagues. Hence, we discussed the need of a community building and raising awareness with respect to potential of implementation research for the field mathematics education.

Replication studies

The interest in replicating didactic designs and empirical studies has been present in the community of mathematics educators for several years. In the 1970s Phillip M. [START_REF] Eastman | Replication studies: Why so few[END_REF] wondered why there were no more studies of replication in our field, while arguing for relevance of such studies. During the 1980s several reproducibility studies were carried outmainly in the French teaching communityfocused on understanding the conditions that allowed a didactic design to be implemented with enough fidelity in different scenarios, preserving the effects in student learning.

Simply put, a replication study can be seen as the attempted repetition of a study or an experiment that has been published in a peer-reviewed journal or book. However, since in the social sciences there cannot be two identical qualitative studies (i.e., there is no duplication), the development of qualitative replication studies involves maintaining certain variables similar to the original design, such as investigating a similar population, using the same didactic design, or applying the same modes and categories of analysis or coding. This kind of replication is known as conceptual replication (Hüffmeier, Mazei, & Schultze, 2016). Yet, this characterization does not exclude the possibility of developing quantitative approaches to replication studies.

Lately, replication studies are gaining the attention of the mathematics education community (and beyond). An example of this is the working group on replication in mathematics education, which met for the first time at the PME 42 conference in Umeå, Sweden. Another indicator is that some specialized journals are beginning to receive and publish replication studies, for example the Journal for Research in Mathematics Education (JRME). However, despite its apparent importance, replication studies have found some obstacles to establishing themselves as an "accepted" type of study in our field. As expressed by Hugh Burkhardt ( 2013) "Replication, a key element in scientific research, is simply not sexy" (p. 225).

There are several arguments for this lack of "sexiness". One is that the academic system privileges the publication of original and innovative works, and devalues works that are not perceived as novel. Additionally, some published research does not provide sufficient methodological and empirical details that allow them to be replicated [START_REF] Schoenfeld | On replications[END_REF].

With these ideas in mind, the second thematic group discussions revolved around the following questions:

1. Can replication studies be useful for the development of implementation research in mathematics education? If yes, in what way? If not, why not? 2. What characteristics should a replication study have in order to be useful for the development of implementation research in mathematics education?

Although there was no definitive answer to these questions, it did appear a feeling in the group that replication studies could help us identify conditions that allow (or prevent) certain innovations to be reproducible in different scenarios, or even help us investigate the effects of a particular treatment under different conditions or populations, which in turn would allow us to advance the implementation of research results. However, it was also acknowledged that for this type of development to happen in the field of educational mathematics different conditions must be in place. For example, it would be necessary to promote a culture of data sharing where the research protocols, data sets, and other elements on which the research is based are not only public, but also shared among researchers through data repositories.

Toward a new definition of implementation research in mathematics education

In TSG23, we feel that we are witnessing the emergence of a new and vibrant research area within the field of mathematics education. In plan to continue looking for a common ground on what do we understand as implementation research in our future work, because once we have a more precise understanding of what implementation research is, we be able to think more precisely about what its methodological implications and challenges are. A first collective attempt to formulate a chain of definitions of the key concepts of "implementation" and "implementation research in mathematics education" was led by Boris Koichu at the last meeting of the group. At the end of this discussion, the group formulated the following proposal.

Implementation is a change-oriented process of adapting and enacting a particular resource (e.g., an idea, a tool, an innovation, a framework, a theory, an action plan, a curriculum, a policy) that occurs in partnership of two communities, a community of the resource proponents (CRP) and a community of the resource adapters (CRA). These communities differ but can intersect. At the beginning of the process, the CRP has the ultimate agency over the resource. The process of adapting a resource by CRA includes some of the following: (1) constructing an agency over the resource, (2) changes in ways of communicating, and (3) changes in practice. Accordingly, implementation research in mathematics education is research that focuses on aspects of implementation, as specified above, in the context of mathematics education.

As we progress in our work, we will certainly refine our working definition. However, we feel that the above proposal adequately captures many ideas and concerns discussed in the group at CERME 11 and that the journey aimed at further theorizing and characterizing implementation research in our field is worth taking.

Figure 1 :

 1 Figure 1: Phases of implementation of research findings into practice Introduction to the papers and posters presented at the TWG23


Morten Elkjaer
email: me@edu.au.dk

Identification and diagnosis of students' mathematical misconceptions in a dynamic online environment

Keywords: Formative assessment, implementation research, misconceptions, diagnostic test

Identification and diagnosis of students' mathematical misconceptions in a dynamic online environment

This poster aims to present the ideas behind the design of a dynamic online tool to assess students' mathematical misconceptions in lower secondary school. These ideas and this tool are part of an industrial PhD project, which is a collaboration between the Danish School of Education, the doctoral candidate, and the company Edulab.

Edulab is a Danish private company that has developed an online mathematical learning platform for the Danish elementary school. In Denmark, 75% of elementary schools subscribe to this platform. In fact, every day Danish school students answer 1.5 million tasks on this online platform. This creates a unique opportunity to implement didactical research results regarding mathematical misconceptions directly into practice. The industrial PhD project aims to provide Edulab with an evidence-based tool to capture students' misconceptions, understood as concept images differing from the formal concept definition [START_REF] Tall | Concept image and concept definition in mathematics with particular reference to limits and continuity[END_REF]. Thereby providing the teachers in lower secondary school with easily accessible formative assessment of their students in the fields of algebra and numeracy in working with equations. [START_REF] Cai | Making classroom implementation an integral part of research[END_REF] claimed that linking intervention with implementation through iterations of design research is a strong way for researchers to have an impact on practice.

From reviewing the literature on misconceptions, it concludes that the community knows a great deal about mathematical misconceptions and learning difficulties related to algebra, equations and numeracy in lower secondary school (e.g., [START_REF] Booth | Misconceptions and learning algebra[END_REF]. Such research-based findings in combination with an online mathematics platform that communicates directly to the end user (i.e., the students) provides a unique opportunity to shorten the implementation time span between the research findings and the end user as was pointed out by [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF]. The project goal is to investigate: How can an online diagnostic tool for lower secondary school be designed, utilizing existing research findings on mathematical misconceptions in algebra and numeracy when working with equations?

To answer the question the aim is to design, and develop an online diagnostic test. Designing an online diagnostic test aimed at lower secondary school students is going to be an iterative process. The literature review will give rise to a categorization of algebraic misconceptions. The idea is to use the 'Concepts in Secondary Mathematics and Science' projects diagnostic tests [START_REF] Hart | Chelsea diagnostic mathematics tests: Teacher's guide[END_REF] in algebra and number operations in a modernized, digital version in order to locate students holding misconceptions fitting this categorization. Based on the algebraic tasks the students have interacted with on the platform, it should be possible for machine learning algorithms to create a "platform behavior" fitting the categorization. It is worth noting that students having greater usage of the platform are a preferable choice for diagnosis in order to create a more precise platform behavior for the different categories of misconceptions. Design-Based Research (e.g., [START_REF] Barab | Design-based research: Putting a stake in the ground[END_REF][START_REF] Disessa | Ontological innovation and the role of theory in design experiments[END_REF] will pave the way for the knowledge generated through iterations of the work with the diagnostic tests and the machine learning algorithms. The projects hypothesis is that math teachers in Denmark, based on collected data, have a unique opportunity to organize, plan and complete their teaching. The point of the continuous diagnosis of all students using Edulab's platform is to catch if students are using the platform as if they are is holding one or more misconception(s). The teachers are provided with a formative assessment report, generated for each student that possibly holds misconceptions, describing the student's behavior and the misconceptions present. Thereby, providing an in-depth explanation of which algebraic concepts and processes in which the students' concept image is conflicting with the concept definition may give the teacher the opportunity to aid the student avoid further progress in constructing a deeper rooted and conflicting concept image. [START_REF] Lucariello | A formative assessment of students' algebraic variable misconceptions[END_REF] emphasize the importance of teachers having access to formative assessment regarding individual student's misconceptions in order not to impede the student's future learning. This formative assessment report should contain ideas and guidelines for how the teacher further teaches the student in the form of literature-based explanations of the origin of the diagnosed misconceptions, and how the concept image may be changed.

Ultimately, the industrial PhD project is looking to generate knowledge and perspectives on how mathematical online diagnostic testing together with mathematical online formative assessment can be implemented in order to help teachers, and thereby students in Danish lower secondary school.