P Lapeyre
email: paule.lapeyre@promes.cnrs.fr

S Blanco

C Caliot

J Dauchet

M El Hafi

R Fournier

O Farges

J Gautrais

M Roger

Monte-Carlo and sensitivity transport models for domain deformation

Keywords:

We address the question of evaluating shape derivatives of objective functions for radiative-transfer engineering involving semi-transparent media. After recalling the standard Monte-Carlo approach to sensitivity estimation and its current limitations, a new method is presented for the specific case of geometrical sensitivities. This method is then tested on configurations with multiple-scattering and absorbing (nonemitting) semi-transparent medium. A new geometrical sensitivity algorithm is presented with full details in order to extend, on several examples, its implementation in complex geometries.

Introduction

The optimization of engineering processes involves an objective function driven by physical mechanisms. In the field of geometrical design, shape optimization models are required to find the extremum of an objective function, denoted here as J( π ) , that depends on a vector of design parameters π . Among optimization methods, the gradient descent method can be used to find a local extremum, and stochastic methods (genetic algorithms, particle swarm optimization) can be used to look for a global extremum [START_REF] Daun | Geometric optimization of radiant enclosures containing specular surfaces[END_REF][START_REF] Farges | Conception optimale de centrales solaires à concentration: application aux centrales à tour et aux installations" beam down[END_REF][START_REF] Rukolaine | The shape gradient of the least-squares objective functional in optimal shape design problems of radiative heat transfer[END_REF] . 1 In any case, the derivative of J( π ) with respect to π is a valuable piece of information for the optimization of engineering processes [START_REF] Caliot | Validation of a monte carlo integral formulation applied to solar facility simulations and use of sensitivities[END_REF][START_REF] Dauchet | The practice of recent radiative transfer monte carlo advances and its contribution to the field of microorganisms cultivation in photobioreactors[END_REF][START_REF] Piaud | Application of monte-carlo sensitivities estimation in solfast-4d[END_REF][START_REF] Shim | Adjoint sensitivity and uncertainty analyses in monte carlo forward calculations[END_REF] .

Our work address process configurations that have a complex geometry and where the radiative transfers are a major component of heat transfer models, and therefore of the objective function [START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF] . Classically, the Monte-Carlo method is regarded as a method of reference in such cases because it remains unbiased even for configurations with high level of complexity in radiative proper-ties and geometry [START_REF] Dauchet | Calculation of the radiative properties of photosynthetic microorganisms[END_REF][START_REF] Delatorre | Monte carlo advances and concentrated solar applications[END_REF][START_REF] Galtier | Radiative transfer and spectroscopic databases: aline-sampling monte carlo approach[END_REF][START_REF] Villefranque | A path-tracing monte carlo library for 3-d radiative transfer in highly resolved cloudy atmospheres[END_REF] . To carry out the computation, the models are stated in an integral formulation, which is then considered as an expectation, which in turn is estimated by sampling, yielding an unbiased estimate of the expectation as well as its confidence interval. Considering sensitivities, a well-known advantage of the Monte-Carlo method is its ability to estimate such expectation and its derivatives by using the very same sampling, avoiding additional computation time [START_REF] Blanco | Short-path statistics and the diffusion approximation[END_REF][START_REF] Giles | Monte carlo evaluation of sensitivities in computational finance[END_REF][START_REF] Melchers | A fast approximate method for parameter sensitivity estimation in monte carlo structural reliability[END_REF][START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF][START_REF] Roger | Applications of sensitivity estimations by monte carlo methods[END_REF][START_REF] Roger | Monte carlo estimations of domain-deformations sensitivitiies for radiative transfer applications[END_REF][START_REF] Tregan | Convergence issues in derivatives of monte carlo nullcollision integral formulations: a solution[END_REF] . This advantage has however some limitations when sensitivity to geometrical parameters is considered because the integral formulation of the model raises formalization and implementation difficulties [START_REF] De La | Calculs de sensibilités par méthode de Monte-Carlo, pour la conception de procédés à énergie solaire concentrée[END_REF] . Those limitations led us to consider the geometrical sensitivities restarting from the local differential equations, namely not formulating first the physical models as integrals.

Several studies have addressed the question of finding shape derivatives for linear physical models through shape sensitivity analysis [START_REF] Sokolowski | Introduction to shape optimization[END_REF] . In the field of deterministic approach, [START_REF] Rukolaine | Shape optimization of radiant enclosures with specular-diffuse surfaces by means of a random search and gradient minimization[END_REF] propose a generalization of the sensitivity models to the radiative transfer equation with specular and diffuse radiative boundary conditions in order to solve the shape optimization problems. In the present work, we extend the question of evaluating shape derivatives when using the Monte-Carlo method as a method of reference. Towards this goal, we propose a construction of the shape sensitivity models that overcome the limitations mentioned above, still preserv-ing the intuitive aspect of the method and its unbiased prediction power. In short, we consider the geometrical sensitivity as a quantity that is transported.

We show how to formulate its transport model [START_REF] Dauchet | Analyse radiative des photobioréacteurs[END_REF] , as well as how to translate the boundary conditions of the physical problem into the corresponding boundary conditions for the sensitivity transport model. Then and only then, the MC method is used to draw estimates from this new transport model. As a result, we gain that geometrical sensitivities can be solved in configurations that were previously unsolvable by classical differentiation of Monte-Carlo integrals. As a side-effect, we lose the well-known advantage of estimating sensitivities using the same sampling -yet at the boundaries only.

In Section 2 , we will briefly review the existing Monte-Carlo method used to estimate an objective function J ( π ) and its derivative ∂ π J ( π ) and the difficulties raised by geometry deformation. In next section, we will introduce our new sensitivity model. We will use there an academic configuration with an uni-dimensional geometrical parameter π in order to expose it in full details.

Finally, the geometric sensitivity model will be applied on four different configurations, presenting the scope of our development so far for complex geometry, reflective boundary conditions and curve surfaces.

Classical approach for sensitivity in Monte Carlo methods

In this section, we briefly summarize the sensitivity calculation as described in [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF] and [START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF] . Any radiative quantity J expressed as a linear integral function can be estimated by the Monte Carlo method as long as it is stated as an expectation ( Eq. ( 1) of the first insert). Formulated as such, the objective function might depend on the parameter π through the density probability functions, the Monte Carlo weight function and the integration domain. One of the advantages of the Monte Carlo method is that the sensitivities of the objective function can always be defined from the same random variable (same domain definition and same probability density functions) as the objective function ( Eq. (2) of the first insert). Numerically, this means that the objective function and its sensitivities are estimated simultaneously from the same set of the random variable samples.

However, when the integration domain is parameterized by π (e.g. due to geometry deformation), two major difficulties have been pointed out in previous works [START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF][START_REF] De La | Calculs de sensibilités par méthode de Monte-Carlo, pour la conception de procédés à énergie solaire concentrée[END_REF] . The first one lies in constructing the formal expression of the deformation velocity vector inside the domain (see Eq. ( 4) in the first insert). The second one arises from the implementation of this formal expression when dealing with multiple scattering and multiple reflection. These two difficulties jeopardize the extension of the existing method to configurations with complex geometry. Therefore, we present an alternative proposal where the geometric sensitivities are obtained through the differentiation of the radiative transfer model itself.

As a last point, we recall that the typical objective functions for geometrical design of (stationary) radiative transfer, when addressed with Monte Carlo, are integrals over frequency, over positions x and over directions u of linear functions of the monochromatic intensity I( x , u ) . The Monte Carlo algorithms therefore start by sampling the corresponding spaces and the algorithm completes by estimating I( x , u ) at the sampled frequency. The four examples that will be presented below will include such preliminary sampling, but all our formal developments will focus on the expressions for I and ∂ π I .

Classical approach to estimate sensitivities inside a Monte Carlo algorithm

The random variable , established over a multiple dimension domain D according to a probability density function p (γ ) , is identified such as:

J(π ) = D (π ) p (γ , π ) dγ ˆ w (γ , π ) = E ˆ w ( , π ) (1)
where ˆ w is a random variable function called the Monte Carlo weight function. When the integration domain is independent of π it is shown in [START_REF] De Lataillade | Monte carlo method and sensitivity estimations[END_REF] that the sensitivity of J(π ) with respect to π is :

∂ π J(π ) = D p (γ , π ) dγ ∂ π p (γ , π ) p (γ , π ) ˆ w (γ , π ) + ∂ π ˆ w (γ , π ) (2)
∂ π J(π ) is an expectation of a new random variable function ˆ w π of the same random variable .

ˆ

w π (γ , π ) = ∂ π p (γ , π ) p (γ , π ) ˆ w (γ , π ) + ∂ π ˆ w (γ , π ) (3)
When the integration domain does depend on π , it is shown in [START_REF] Roger | Modèles de sensibilité dans le cadre de la méthode de Monte-Carlo : illustrations en transfert radiatif[END_REF][START_REF] Roger | Monte carlo estimates of domain-deformation sensitivities[END_REF] that the sensitivity of J(π ) with respect to π can be written as :

∂ π J(π ) = D (π ) p (γ , π ) dγ × ∂ π ˆ w (γ , π ) + ˆ w (γ , π ) ∂ π p (γ , π ) p (γ , π ) + ∇ . ˆ w (γ , π ) p (γ , π ) v π p (γ , π ) (4)
where v π denote the deformation velocity. This velocity is constrained only at the integration domain boundaries, and must be extended arbitrarily all over the domain.

An alternative proposal : The sensitivity transport model

In this section, the intensity geometric sensitivity ∂ π I is introduced at the modeling level. That is, instead of derivating I ( π ) stated as an integral, we view the sensitivity as a quantity with its own transport model. This model is obtained by derivating the radiative transfer equation (RTE). In the present work, we consider the case of a scattering, absorbing non-emissive medium (adding medium emission would not convey any additional idea).

Denoting the domain , and the boundary ∂ ( π ), the RTE then reads :

u . ∇ I( x , u , π ) = -(k a + k s ) I( x , u , π ) + k s 4 π p( u | u ) d u I( x , u , π ) for x ∈ I( x , u , π ) = I ∂ (π ) ( x , u , π ) for x ∈ ∂ (π ) and u . n > 0 (5) The RTE differentiation leads to u . ∇ s ( x , u , π ) = -(k a + k s ) s ( x , u , π ) + k s 4 π p( u | u ) d u s ( x , u , π ) for x ∈ s ( x , u , π ) = s ∂ (π ) ( x , u , π ) for x ∈ ∂ (π ) and u . n > 0 (6)
where s ( x , u , π ) = ∂ π I( x , u , π ) . Since the geometrical parameter π only appears at the boundary conditions, it does not affect the operators of Eq. ( 5) . We can therefore notice that the sensitivity transport equation is identical to the intensity transport equation in the domain . It implies that, in the domain, the algorithms used to solve the RTE can be used unchanged to solve the sensi-tivity transport equation. The sensitivity transport and the intensity transport only differ over boundary conditions. Therefore sensitivity can be regarded as a quantity that is transported by the same transport model as the intensity, only changing the sources at boundaries. The question that remains is then to determine the sources of sensitivity at boundaries. We show below that sensitivity boundary conditions involve a coupling with intensity and depend on the radiative properties of the boundary ( Eq. ( 9) ) so that the sensitivity transport can not be solved on its own. In the context of Monte Carlo formulation, this will result in an integral formulation of the sensitivity which includes Eqs. ( 5) and (6) (details in Section 4 ), which requires in turn to solve systems ( 5) and ( 6) in parallel.

General expression for sources of sensitivity at boundary

Let us denote I B the field of the incoming radiative intensity at the boundary. All Monte Carlo approaches rely on an integral formulation of the RTE solution of the form

I( x , u , π ) = D (π ) p (γ | x , u ) dγ ˆ w ( γ , I B ) ( 7
)
where p (γ | x , u ) is the probability density of path γ propagating the radiative sources from the boundary through the medium up to x in direction u . Identically, let us denote s B the field of the incoming sensitivity at the boundary. The Monte Carlo algorithm evaluating the intensity is straightforwardly translated into a Monte Carlo algorithm evaluating the sensitivity based on the following integral formulation in which only the boundary sources are changed:

s ( x , u , π ) = D π (π ) p π (γ π | x , u ) dγ π ˆ w ( γ π , s B ) (8)
The absorbing and scattering medium properties are preserved during the process and therefore the space of optical paths is identical to that of the sensitivity paths π . The main point here is that the sensitivity boundary conditions are transported by the same paths as the intensity boundary conditions and we can therefore use the same set of sampled paths to estimate simultaneously both quantities.

However at this point the boundary conditions of the sensitivity model are unknown. The main objective of constructing the sensitivity model is to formulate those boundary conditions: the sensitivities associated to the geometrical perturbations for all incoming directions at the boundary. Preserving the complete generality of the emission, absorption, reflection and scattering phenomena compatible with the radiative transfers equation, these incoming sensitivities can be written under the form of a linear application L of • the sensitivities s ( x , u , π ) in all outgoing directions u (as for any surface reflection problem, only here it is the sensitivity that is reflected),

• the intensities I( x , u , π ) also in all outgoing direction u ,

• and the black-body intensity at the local temperature

I b (T ( x , π )) . s ∂ (π ) ( x , u , π ) = L s ( x , u , π ) , I( x , u , π ) , I b (T ( x , π )) with u . n > 0 (9)
The main point here is the fact that I appear in this linear application. This implies that via the boundary condition, the sensitivity model is coupled with the intensity model: the sources of sensitivity have to be evaluated, from the outgoing intensity ( Fig. 1 ), before being propagated into the domain. Appendix A and Appendix B provide practical examples of these coupling at the

A roadmap to the sensitivity transport model in the Monte Carlo context

• The sensitivity transport equation is well established and familiar: it represents the same transport phenomena as for the intensity ( Section 3 ).

• Regarding the boundary conditions, as far as our configurations are concerned: sources of sensitivities appear at geometrically modified boundaries ( Section 5 ):

the impact of the volume of semi-transparent medium added or withdrawn when modifying the surface ; the term u . ∇ I of Eq. ( 18) is directly derived from the radiative transfer equation ( Section 5, Appendix A ). sources of sensitivity appear as Dirac distributions where intensity displays discontinuities at the boundary ( Section 5, Appendix A ) everywhere the boundary condition is reflecting for the intensity, it is as well reflecting for the sensitivity

• From an algorithmic point of view:

well-identified difficulties arise from sources expressed as Dirac distribution in location and/or direction ( Section 5 ), the coupling between sensitivity transport equation and intensity transport equation through the boundary conditions is dealt with standard double randomization technique ( Section 4 ).

In configuration 1 , the source of radiation is a black cylinder inside a square cavity filled with absorbing and scattering material. The objective function is the local absorption power and the geometric parameter π acts on the size of the square cavity. This configuration is discussed in full details in Section 3 in order to illustrate the entire process of designing a shape sensitivity algorithm: expressing the sensitivity sources (as function of the incoming intensities); choosing a Monte Carlo algorithm propagating these sources as if they were known (according to the standard practice of Monte Carlo in radiative transfer); coupling it with a Monte Carlo algorithm evaluating I . In the following section (5) three other configurations are successively discussed. Configuration 2 is very similar to configuration 1 but now the square cavity is emitting with a temperature discontinuity. This leads to the definition of a new type of sensitivity sources: Dirac sources at the discontinuity location. Configuration 3 deals with the same family of sensitivity sources but they are now related to reflection at geometrically modified boundaries and the technique is scaled up to realistic thermal solar plants. The Dirac sources will here be at the edge of heliostats of increasing sizes. With configuration 4 we go back to semi-transparent media, still with complex geometries but now concentrating on the question of evaluating sources at altered curved surfaces.

Fig. 2.

On the left, the square cavity geometry filled by a semi-transparent medium, lightened by an emissive cylinder at its center. On the right, the dashed line represents the uniform scaling of the cavity. The solid multiple-scattering path is a typical realization of a radiative path use in a Monte Carlo algorithm evaluating the local absorption power at x 0 . The same path is used to evaluate the sensitivity. The dotted multiple-scattering path illustrates how we translate in Monte Carlo terms the coupling of the sensitivity model with the initial radiative transfer model: the incoming sensitivity at x γ requires the evaluation of the outgoing intensity at the same location and thanks to double randomization, only one single sample of a standard Monte Carlo path is required to complete the sensitivity realization.

Foundational algorithm (configuration 1)

Here, we expose the algorithmic principles of the method, using, for the sake of clarity, a fairly simple configuration with a simple geometrical alteration and homogeneous medium and boundaries.

Let us consider a square cavity containing a gray semitransparent absorbing, scattering, but non-emitting medium. Cavity dimensions are set along the x and y axes and are infinite along the z axis. At the center of the square an infinite isothermal blackbody cylinder at temperature T . The square walls are also black but emit no radiation ( Fig. 2 ). The semi-transparent domain is denoted as , and its boundary ∂ ( π ) is constituted by the cavity boundary R (π ) and the cylinder emitting boundary F . The geometrical parameter π is setting the cavity boundary R (π ) : the deformation is a colinear translation of R (π ) according to the normal vector n R (π ) resulting in an uniform scaling of the cavity. The objective function is the local absorption power J( x 0 , π ) at a location x 0 and we aim at evaluating its derivative ∂ π J( x 0 , π ) with respect to π . Both are integrals over the unit sphere of I and s :

J( x 0 , π ) = 4 π k a I( x 0 , u 0 , π ) d u 0 (10) ∂ π J( x 0 , π ) = 4 π k a s ( x 0 , u 0 , π ) d u 0 (11)
In the Monte Carlo context, these integrals will be evaluated by first sampling a direction u 0 and then evaluating I( x 0 , u 0 , π ) and s ( x 0 , u 0 , π ) . The radiative model is the RTE of Eq. ( 5) with the boundary conditions:

I( x R , u , π ) = 0 I( x F , u ) = I b (T ) (12)

The sensitivity model

With regard to the square cavity deformation and considering the boundary radiative properties, the linear application L describing sensitivity boundary conditions is here

s ∂ (π ) ( x F , u ) = 0 ( 13
)
because F is unaffected when changing π , and

s ∂ (π ) ( x R (π ) , u , π ) = - k s u . n R 4 π p( u | u ) d u I( x , u , π ) (14)
We will discuss the expressions of such sensitivities for incoming directions at geometrically modified boundaries with more generality in Section 5 and Appendix A . The main point is that they include the intensity in outgoing directions. Here the reason is the following: when enlarging the cavity (for instance), a new layer of semi-transparent material is introduced and this layer can interact with outgoing radiation and scatter it back into the cavity.

Propagating the sensitivity sources with a standard Monte Carlo approach

The sensitivity boundary condition is coupled with the outgoing intensity at all points on the geometrically modified boundaries. This means that we have to evaluate this intensity to set the boundary condition. Let us temporarily assume that the intensity is known and therefore that the sensitivity is also known for all incoming directions. We can then consider these incoming sensitivities as sources, as in a standard radiative transfer problem, and think of Monte Carlo algorithms propagating them. Here we want to evaluate s at a location x 0 and direction u 0 . We know that s is solution of the very same RTE as I and the standard Monte Carlo approach would therefore be to use a reverse algorithm, exactly as for the estimation of I( x 0 , u 0 , π ) . This means that the algorithms evaluating I( x 0 , u 0 , π ) and s ( x 0 , u 0 , π ) start with the very same steps, and therefore can share the same sampled paths:

• A multiple scattering paths γ is sampled, starting at x 0 in directionu 0 , finding successive scattering locations x i according to Beer law for pure scattering, scattering directionsu i according to the single scattering phase function, until reaching one of the two absorbing solids at location x γ in directionu γ .

• Two Monte Carlo weights are then computed using the same path:

• one for the intensity problem that is simply the value of the incomming intensity at x γ in direction u γ attenuated by continuous absorption along the path length l γ , i.e.

I( x γ , u γ , π ) exp(-k a l γ ) ,
• one for the sensitivity problem constructed the very same way but with the incomming sensitivity, i.e.

s ( x γ , u γ , π ) exp(-k a l γ ) .

Coupling the sensitivity Monte Carlo algorithm with an estimate of the outgoing intensity

Since the sensitivity sources are unknown at this stage, the algorithm described above is incomplete for its sensitivity part: we would need to evaluate the intensity at each location x γ where paths meet boundaries, and in all directions u . We could embed a full new Monte Carlo computation of I for such each location x γ , but this is unnecessary due to one fundamental property of Monte-Carlo methods: the double randomization. Since the sensitivity is expressed as an expectation of a linear function of the intensity that is itself expressed as an expectation, both can be nested into one expectation [START_REF] Dauchet | Addressing nonlinearities in monte carlo[END_REF] . Hence, the sensitivity can be estimated with no bias using one single Monte Carlo algorithm to estimate the latter, combining sensitivity and intensity. This is the same idea as in reverse multiple scattering or multiple refection Monte Carlo algorithms. In such very standard algorithms, when a collision location has been sampled, at a surface or in the volume, fully computing the reflection or scattering sources would require to start a Monte Carlo algorithm with a large number of sampled paths for each unknown intensity in each incoming direction. But instead, only one reflection or scattering direction is sampled and a multiple scattering path is initiated for this direction only. Here the process is less intuitive because it combines sensitivities and intensities, but it is strictly similar:

• a path γ is sampled, as described above, from ( x 0 , u 0 ) to ( x γ ,u γ ) ; • the required sensitivity s ( x γ , u γ , π ) is defined in Eq. ( 14) as an integral over directions u and the algorithm therefore samples u (just as sampling the phase function in a multiple scattering algorithm);

• a new path γ is started in directionu as if evaluating I( x γ , u , π ) with a reverse Monte Carlo algorithm.

The resulting algorithm is described in 15) and ( 16) ) to a uniform scaling of R (π ) . By comparison with previous simulation examples, the intensity is here integrated over the sensor's surface and its response frequency range. The sensitivity results show that the additional integration process has no influence on the sensitivity estimation. The sensor of square surface S sensor is located at location x 0 ( Fig. 2 ) and oriented toward thee x direction. The geometrical change of the cavity is the same as in Fig. 3 . The wavelength range is [ λmax when γ hits the square, u is outgoing, and γ hits the cylinder. Fig. 3 displays the simulation results and illustrate the stability of sensitivities estimated by the model. For comparison, finite differences were calculated from Monte-Carlo estimation of the absorbed radiative power density for different positions of the upper surface. The latter results illustrate the typical difficulties encountered when evaluating the sensitivities from finite differences, that is the difficulty to obtain converged results. In this regard, if δL is chosen too small with respect to L , the variance of the finite difference estimate becomes too large and the results are inaccurate.

On the contrary, if δL is too large, the finite difference converges but toward a value that departs from the true sensitivity. This difficulty is well known when gradients are estimated by differentiation (it is outlined in the gradient-based Kiefer-Wolfowitz method presented in [START_REF] Daun | Geometric optimization of radiant enclosures containing specular surfaces[END_REF] ). The figure also shows that the local absorption power becomes less sensitive to a variation of δL L when the size of the cavity increases. Table 1 summarizes the results of the local absorption power and its shape sensitivity for different parameter sets. The closer the probe location gets to the R (π ) boundary the higher the sensitivity. It also shows a practical limit to the evaluation of the sensitivities in the case of high optical thicknesses. Apart from this limit, there are no convergence issues.

Surface and spectral integration

The methodology can be extended to a spectral and surface objective function, for example to simulate the response measured, at the position x 0 within the cavity, by a sensor of surface S , and for a given frequency range. As long as the measure settings are independent of the deformation parameter, the question of evaluating the sensitivity as presented in Section 3 and the corresponding methodology ( Section 4 ) are not impacted. In terms of Monte-Carlo, evaluating the sensor response sensitivity to a perturbation of the square surface is equivalent to sample a position on the sensor surface, a wavelength within the specific frequency range and a direction within the measured hemisphere prior to sample the

Table 1

Absorbed radiative intensity density and its sensitivity results for a fixed value of L . The probability density function for isotropic scattering is p U ( u i ) = 1 4 π and probability density function of scattering extinction is p (l i ) = k s exp (-k s l i ) . Location x 0 on the y axis is determined by y 0 L = 0 . 5 and the ratio between the cylinder radius and L is set at r L = 0 . 125 . The number of realizations of the Monte-Carlo weight functions is N = 2 × 10 6 , (%) represents the relative error of ∂ π J ( π) and J ( π). sensitivity path ( Eqs. ( 15) and ( 16) ).

k a L k s L x 0 L ∂ π J( x 0 ,π ) L 4 π k a I b max σ L 4 π k a I b max ∂ π J (% ) J( x 0 ,π ) 4 π k a I b max l σ 4 π k a I b max J (%)

J( x

0 , π ) = S d S λmax λ min d λ 2 π k a (λ) I( x 0 , u 0 , λ, π ) d u 0 (15) ∂ π J( x 0 , π ) = S d S λmax λ min d λ 2 π k a (λ) s ( x 0 , u 0 , λ, π ) d u 0 ( 16
)
The typical integrated radiative objective functions do not introduce additional difficulties regarding the sensitivity. Fig. 4 displays the corresponding simulation results.

The different types of sensitivity sources

In this section we focus on the sensitivity sources. So far we have seen that beside the coupling at the boundary condition, the sensitivity model can be estimated in the domain exactly like the intensity model, using standard Monte Carlo approaches to radiative transfer. The following examples will all make use of most simple forward Monte Carlo algorithms. More attention will be devoted to the boundary condition definitions, in particular to the Dirac sources due to boundary discontinuities and to sources associated to curved surfaces.

We will see that depending on the configurations, the sensitivity sources can take diverse forms, namely that the linear application expressing the boundary conditions can be composed of distinct linear operators. In Appendix A , formal developments are provided and reported here by Eq. ( 17) for a translation colinear to the normal vector :

s ( x ∂ (π ) , u , π ) = -∂ y I( x ∂ (π ) , u , π ) v + ∂ π I p (x, z, u , π ) (17)
This formulation establishes a relation between the shape derivative s ( x ∂ (π ) , u , π ) and the material derivative ∂ π I p (x, z, u , π ) in this particular transformation [START_REF] Sokolowski | Introduction to shape optimization[END_REF] . Using this relation for blackbody radiative boundary conditions, we identify two types of contributions to the sensitivity in incoming directions:

s ( x ∂ (π ) , u , π ) = - u . ∇ (I( x ∂ (π ) , u , π ) u . n v + u . ∇ T I( x ∂ (π ) , u , π ) u . n v with u . n > 0 (18)
The first contribution type (first term in Eq. ( 18) ) corresponds to the effect of medium diffusion and medium absorption in the infinitesimal layer along the surface added or suppressed after the surface deformation (cf Fig. 5 a). If a surface is shifted toward the outside, its displacement creates a new layer in its vicinity (the interior volume is increased). The intensity exiting the surface in direction u will be attenuated by this layer (absorption and scattering) and furthermore the intensity incoming on this new layer will have a probability to be scattered in direction u . When translated into sensitivities, these decreases and increases of the incoming intensity at the boundary lead to exactly the same terms as the collision terms in the radiative transfer equation ( -

(k a + k s ) I( x , u , π ) + k s 4 π p( u | u ) d u I( x , u , π
) ), which allows to summarize them using the transport operator, which explains the term u . ∇ I in Eq. ( 18) (transport equals collision in the stationary RTE).

In configuration 1 , only this first contribution term was encountered. The second term of Eq. ( 18) corresponds to the effects of non homogeneous radiation in incoming directions at a geometrically modified boundary. We use the notation ∇ T I for the surface gradient, also called tangent gradient, of the incoming intensity at the boundary. Fig. 5 b illustrates this effect on the intensity at location x ∂ (π ) in direction u . Let us assume that the medium is transparent: before the deformation, the intensity in direction u was exactly the intensity emitted or reflected by the surface at that point. After deformation, this intensity comes from another location on the surface which explains the tangential gradient in Eq. ( 18) . This tangential gradient plays a very specific role when dealing with surface discontinuities, for instance when two adjacent surfaces have different temperatures ( Fig. 6 a), or different reflection propor even when they have the same reflection properties but their normals are not the same. In all three cases, the intensity in any incoming direction is discontinuous when modified along the surface and the corresponding Heaviside function leads to a Dirac when applying the tangential gradient. This creates sources that are not continuously distributed along the surface but concentrated on the discontinuity lines: we then have to deal with Dirac sources of sensitivity.

The same idea of a discontinuity of intensities at modified boundaries allows to think the sensitivity sources associated to shading effects. This is illustrated in Fig. 6 b with a curved surface: at point p some directions are blocked by the surface and if the point is displaced along the y axes, it will block new directions and consequently become a local source of sensitivity. The same shading effect occurs with adjacent plane surfaces with different normal orientations as soon as the discontinuity p is geometrically modified (see 6 c).

For configuration 2 , we handle formally this specific familly of sensitivity sources, starting from Eq. ( 18) and we give rigorous definitions of Dirac sources of sensitivity. For further cases, the formal statement of those sensitivity sources at discontinuities is still a work in progress. We do not address them yet using the general sensitivity boundary condition stated in Eq. ( 18) but they are estimated on an ad hoc basis by analyzing geometrically the differential implications of surfaces deformations. This is the case for the curved and reflective surface examples of configuration 3 and configuration 4 (see Appendix B ).

In the following of this section, configuration 2 gives an example of an edge boundary discontinuity treated from the general sensitivity boundary condition statement Eq. ( 18) in an academic case: the same square cavity as in part 4 but with a different deformation. In configuration 3 , the intensity edge discontinuity is extended on a thermal solar configuration with reflective surfaces. It is showed that this particular type of sensitivity sources can have a heavy impact on the algorithm when considering complex geometry. We recognize here similar difficulties as reported in another work about shape sensitivities [START_REF] Li | Differentiable monte carlo ray tracing through edge sampling[END_REF] where optimizing the algorithm process demanded great effort s in the case of image synthesis transparent medium configurations. The boundary condition statement of configuration 3 has been done on an ad hoc basis in Appendix B . Finally configuration 4 gives an example of a shading discontinuity in the case of curved surface.

Configuration 2

Let us consider another square cavity which dimensions are set along the x and y axes and are infinite along the z axis. The medium inside the square, in the domain , is semi-transparent. The four walls of the square cavities are black-bodies and only the upper wall is emitting ( R up ) : The geometrical parameter π induces an extension of the height of the square cavity by a displacement of the upper-wall in a direction colinear to the upward normal vector n R up ( Fig. 7 ). The objective function is the radiative flux crossing the bottom-wall ( F ) and we want to estimate its sensitivity to π . The sensitivity source is only at the modified boundary, i.e. along the upper wall. It has four contributions associated to:

⎧ ⎪ ⎨ ⎪ ⎩ I( x F , u , π ) = 0 or x F = (x, 0 , z) I( x R l , u , π ) = 0 or x R l = (0 , y, z) I( x R r , u , π ) = 0 or x R r = (L, y, z) I( x R up , u , π ) = I b (T ) or x R up = (x, L, z) (19)
• extinction of radiation by the medium in the vicinity of the modified boundary, • scattering sources in the same vicinity,

• temperature discontinuity at the left side of the modified boundary, • temperature discontinuity at the right side of the modified boundary.

The last two contributions are Dirac sources (involving a Dirac distribution at 0 and L ):

s ∂ (π ) ( x R (π ) , u , π ) = - k a + k s u . n R up I b (T ) + k s u . n R up 4 π p( u | u ) d u I( x , u , π ) + I b ( x , u , π , T ) u . e x u . n R up δ(x -0) H ( L -x ) H( u . e x ) -I b ( x , u , π , T ) u . e x u . n R up δ(x -L ) H ( x -0 ) H(-u . e x ) (20)
The second contribution was already discussed in Configuration 1 . The first contribution did not appear because the square surface was neither emitting nor reflecting (no incoming radiation). Here it appears because the geometrically modified boundary is a black-body at temperature T . Both the first and second contributions are distributed all along the upper wall. The other two contributions are at the singular infinite corner lines (dots in Fig. 7 b). Because of the Heaviside terms H( u . e x ) and H(-u . e x ) , for a given direction u , only one of the two is active: that corresponding to the corner from which u is a valid incoming direction. They correspond to the fact that from a given location inside the domain, when considering radiation coming from the corner exactly, increasing π (i.e. raising the top surface) implies that the nonemitting lateral face is viewed, whereas lowering it implies that the emitting top surface is viewed. The Dirac sources of sensitivity translate this discontinuity in distribution terms.

We implemented a very standard forward Monte Carlo algorithm in which these sources are propagated from the upper-wall into the domain along multiple-scattering paths, attenuated exponentially with respect to the absorption optical thickness of the path. Only the paths reaching the bottom-wall contribute to the sensitivity estimate. Fig. 8 displays the results for both the bottomwall radiative flux as function of the top surface height. In the same graph the sensitivity to the height is plotted twice: as estimated using the present sensitivity Monte Carlo algorithm and evaluated from finite differences. Using the sensitivity model convergence is less than 1% for 10 6 sampled paths, in contrast with the finite difference estimates that are much less accurate as explained in Section 4 .

Configuration 3

In this configuration shape sensitivities are implemented for a concentrating solar tower application [START_REF] Farges | Conception optimale de centrales solaires à concentration: application aux centrales à tour et aux installations" beam down[END_REF] . The configuration is illustrated on Fig. 9 . The mirrors are specular reflecting flat squares and the surrounding medium is transparent. The intensity boundary conditions for all incoming directions are stated on each of the following geometrical objects: the sun S, the tower T , the thermal receiver R and the heliostats H:

⎧ ⎪ ⎨ ⎪ ⎩ I( x S , u , π ) = I s I( x T , u , π ) = 0 I( x R , u , π ) = 0 I( x H , u , π ) = ρ I( x H , u spec , π ) (21)
The incoming intensities are null at all the other limits (sky and ground). The simulation objective is to evaluate the radiative flux at the receiver and understand the impact of the heliostat field design in terms of the following three optical phenomena: an optical path from a point in the solar disk to a point on one of the heliostats may

• be shaded by a surface located between the sun and the heliostat (the shading phenomenon); • be blocked by a surface located between the reflection location and the receiver (the blocking phenomenon); • miss the receiver due to optical errors during the reflection (the spillage phenomenon).

We here study the influence of the variation of the heliostat size (which is modified according to the parameter π , Fig. 10 ). The As they are the only deformed surfaces, in terms of sensitivity modeling the heliostats mirrors are the only sensitivity sources. On each heliostat the incoming sensitivity has five contributions:

• one associated to the reflection of sensitivity sources coming from other heliostats; • four associated to the extension of each of the heliostat edges.

Using a Cartesian coordinate system ( O, x, y ) centered on the heliostat and aligned with edges,

s ∂ (π ) ( x H(π ) , u , π ) = ρs ( x , u spec , π ) H( π 2 -x ) H( π 2 -y ) H π 2 + x H π 2 + y + ( ρI( x , u spec , π ) -I( x , u , π ) ) 2 δ π 2 -x H π 2 -y ×H π 2 + x H π 2 + y + ( ρI( x , u spec , π ) -I( x , u , π ) ) 2 δ( π 2 -y ) H π 2 -x ×H π 2 + x H π 2 + y + ( ρI( x , u spec , π ) -I( x , u , π ) ) 2 δ π 2 + x ×H π 2 -y H π 2 -x H π 2 + y + ( ρI( x , u spec , π ) -I( x , u , π ) ) 2 δ π 2 + y ×H π 2 -x H π 2 + x H π 2 -y (22)
The first contribution was not encountered in configuration 1 and configuration 2 because none of the surfaces were reflecting. Here for the sensitivity model, the reflecting heliostats simply reflect the impacting sensitivities exactly as they reflect intensity in the initial radiative transfer model. It is not always the case. When the geometric parameter creates a translation of the reflective boundary, other terms appear (not shown). See, for instance, the generalization of reflective sensitivity boundary conditions that have been done in [START_REF] Rukolaine | Shape optimization of radiant enclosures with specular-diffuse surfaces by means of a random search and gradient minimization[END_REF] for deterministic method. The four other terms are very similar to the Dirac sources of configuration 2 : they translate in sensitivity terms the intensity discontinuity at the limit of modified boundary.

As far as the sensitivity Monte Carlo implementation is concerned, we used a forward approach, which was not the case for the initial radiative transfer algorithm. Indeed, when evaluating the flux collected by the receiver, standard concentrated solar algorithm do not start at the sun: the first step is the sampling of a location on the heliostat surfaces and only then a location on the solar disk is sampled to represent the solar radiative sources in a backward Monte Carlo approach. Here for the sensitivity, we stick to a direct sampling of the sensitivity sources themselves, without any intermediate statistical step. This means that we start by sampling a location along all the edges of all the heliostats. The sensitivity sources are indeed only at the edges: the first contribution in [START_REF] Shim | Adjoint sensitivity and uncertainty analyses in monte carlo forward calculations[END_REF] , spread over all the heliostat surface, is only a reflection of sensitivities emitted elsewhere (at the edge of an heliostat) and propagated to the mirror surface. Therefore in a forward Monte Carlo algorithm for sensitivities, paths will be started at the edges and refection will simply occur when a path encounters a mirror. The only subtle point is that the sources at each edge location involve the knowledge of two intensities: one in the direction u spec before reflection, one in the considered direction u corresponding to radiation propagating at the edge location without impacting it, therefore coming from another heliostat. These two intensities are not known and double randomization is used (as already illustrated in Sections 4 and 5.1 ), meaning that two radiative paths are initiated at the sampled edge location. Apart from this double randomization technique, that will be required for all sensitivity algorithms involving either reflection or scattering, the sensitivity algorithm is fully standard, the sources being propagated along straight lines, reflected by the heliostats, until they reach the receiver or are stopped by an absorbing surface.

Fig. 11 displays simulation results for an heliostat field including 250 heliostat set up with a radial staggered layout [START_REF] Siala | Mathematical formulation of a graphical method for a no-blocking heliostat field layout[END_REF] . The thermal power sensitivities are estimated for an heliostat size range between 0.8 to 1.2 m with a step of 0.02 m. The results presented in the Fig. 11 start at the value of 0.8 m. Indeed when there is no overlap of the heliostats (blocking, shading) the collected power is proportional to the size of the heliostats. This results in a constant sensitivity of the collected power to the size of the heliostats. In this figure we therefore only present the sizes of heliostat mirrors which induce overlap.

The results have been cross validated with finite differences and show an exact match. The standard error is significantly smaller (0.03% for 10 6 samples) than that of the finite difference estimate, as explained in Section 4 . Once the double randomization scheme was clarified, there were no particular difficulties associated to the implementation itself and no specific convergence issues were identified.

Configuration 4

In this configuration we consider a sphere, described in parametric terms, inside a surpershape geometry, descried as a large ensemble of triangles. The sphere is stated as an isothermal emitting black-body at temperature T and the boundary of the supershape is also a black-body but not emitting. The medium between the supershape and the sphere is semi-transparent and non emitting.

The objective function is the radiative power received by a part of the supershape boundary: that is all the supershape surface above the slice in the Fig. 12 b. The goal is to estimate its sensitivity to the sphere radius. The intensity model is stated from the RTE with Lambertian boundary condition at the sphere surface :

I( x Ssphere , u , π ) = I b (T ) with u . n sphere > 0 I( x Sshape , u , π ) = 0 with u . n shape > 0 ( 23
)
with all the normal vectors oriented toward the medium. The radiative power emitted by the sphere surface is 4 π r 2 σ T 4 . This quantity is attenuated during its transport in the semi-transparent medium. Since only the sphere radius r is influenced by the geometrical parameter π only the sphere surface will be emitting sensitivity. The sensitivity source has three distinct contributions associated to:

• extinction of radiation by the medium in the vicinity of the sphere, • scattering sources in the same vicinity,

• shadowing of tangent radiation.

Using standard notations of polar and azimutal angles, θ and φ for incoming directions, θ and φ for outgoing directions,

s ∂ (π ) ( x R (π ) , θ , φ, π ) = - k a + k s cos (θ ) I b (T ) + k s cos (θ ) 2 π d φ π 2 0 d θ p(θ , φ| θ , φ ) I( x , θ , φ , π ) + 1 r σ T 4 π -I( x , θ , φ, π ) δ(θ -π 2 ) cos (θ ) sin (θ ) (24)
The three contributions are distributed all over the surface and the distinction between them is essentially angular. The first two terms of Eq. ( 24) corresponds to the effect of medium diffusion and absorption in the additional infinitesimal volume. The third one only exist in the tangent directions (see the term in δ(θπ

2 )

in Eq. ( 24) ) and come from the intensity discontinuity between the boundary emission and the out-coming intensity at the same point, in that direction. We implemented a very standard forward Monte Carlo algorithm in which these contributions are propagated throughout the medium according to a multiple scattering algorithm, with continuous extinction by absorption along the path, until reaching on of the absorbing boundaries. The path only contributes to the sensitivity if it ends at the surface of interest (the surface above the slice in Fig. 12 b). Fig. 13 displays the results of the Monte Carlo estimates validated against the finite difference estimates. We encountered no algorithmic difficulty and there was no convergence issue. A perspective for curved parametric surfaces is the formal identification of the boundary condition by the methodology presented in Appendix A , which we could not yet achieved ( Eq. ( 24) was only constructed in an ad-hoc manner. This would open the way to a similar handling of any smooth surface.

Conclusion

As Monte Carlo starts from expressing the objective function as an integral (an expectation), the first approaches to the evaluation of sensitivities inside a Monte Carlo algorithm consisted in taking the derivative of this integral. This lead to a solution that although rigorously justified in full general terms (the deformation velocity method), encounters implementation difficulties when dealing with geometrical sensitivities. Here we presented an alternative approach to geometrical sensitivity estimation. It is shown to allow the handling of geometrical sensitivity problems that were technically assumed unsolvable with the deformation velocity method. The last implementation examples deal in particular with both a discretized complex geometry and curved parametric surface. However, this development comes with an algorithmic cost. While the deformation velocity method made it possible to compute the objective function and its sensitivity at the same time, using exactly the same sampled paths, the sensitivity model method leads to a different and more complex algorithm for the objective function sensitivity than for the objective function itself. Since both sensitivity and intensity share the same transport model, the objective function and its sensitivity are addressed simultaneously in the first part of the algorithm, until a geometrically modified boundary is encountered. But then an additional simulation time is required to deal with the coupling with intensity in the sensitivity boundary condition.

So far, there are still several restrictions due to formal developments that are still ongoing. We addressed only one example of parametric curved sensitivities (the sphere of increasing radius) and our only example involving reflection maintains the modified boundary in its origin plane (only its extension is changed). We also did not yet encounter the convergence difficulties associated with the Dirac sources at surface discontinuities such as those reported in [START_REF] Li | Differentiable monte carlo ray tracing through edge sampling[END_REF] for image synthesis applications (in complex geometries without participating media). But they are likely to occur and it will certainly be very useful to make use of the importance sampling solutions implemented in this reference to bypass such convergence issues. Up to know, considering the restricted class of the configuration that we explored, the sensitivity simulations displayed no specific convergence difficulties.

As a perspective, we consider the extension of the sensitivity model to coupled heat transfer models (conduction, convection and radiative transfers). For the moment, we have been working within the specific framework of the radiative transfer equation. Our proposition of the sensitivity model will be examined in the context of the Monte-Carlo method used to treat heat coupled transfer. At this embryonic stage of the work, the difficulty is to estimate boundary conditions gradients. This project has received funding from the Occitanie region and University of Perpignan Via Domitia, from the programme investments for the future of the French National Agency for Research (ANR) under award number ANR-10-LBX-22-01-SOLSTICE. We acknowledge support from the ANR, grant HIGH-TUNE ANR-16-CE01-0010, ( http: //www.umr-cnrm.fr/high-tune and ANR, grant MGC-RAD) and from Region Occitanie (Projet CLE-2016 ED-Star).

Fig. 1 .

 1 Fig. 1. Schema corresponding to the sensitivity boundary condition ( Eq. (9) ).

  boundary. These examples correspond to the four configurations studied in next sections.

Fig. 3 .

 3 Fig. 3. Absorbed radiative intensity density J ( π) and its sensitivity ∂ π J ( π) to the uniform scaling of R (π ) boundary. The phase function in the semi-transparent medium is isotropic p U ( u i ) = 1 4 π and the probability density function of scattering extinction is p (l i ) = k s exp (-k s l i ) . Cavity dimensions are set by a variation δL of L influencing the optical thickness k e L . Single scattering albedo is uniform on the cavity k s k e = 0 . 5 . Estimations of the absorbed radiative intensity density and its sensitivity are obtained from 2.10 6 realizations of the corresponding Monte-Carlo weight function.

  Fig. C.16 . The sampling of u is isotropic. If u is oriented toward the inside of the cavity (scalar product u . n R (π ) > 0 ), the radiative intensity value is known: it is null since the square surface is neither emitting nor reflecting (box B in Fig. C.16 ). So the path γ evaluating intensity is only required when u is an outgoing direction (dotted path in

Fig. 4 .

 4 Fig. 4. Sensor's response J ( π) and its sensitivity ∂ π J ( π) ( Eqs. (15) and (16) ) to a uniform scaling of R (π ) . By comparison with previous simulation examples, the

2 ; 2

 2 λmax ] where λmax (T ) is the wavelength of the Plank function maximum at the cylinder temperature T = 10 0 0 K. The dependencies with wavelength of the absorption and scattering coefficients are k a (λ) = k s (λ) = (1 /λmax ) λ.

Fig. 2 )

 2 Fig. 2 ). Box C in Fig. C.16 displays an example of the only cases where the sensitivity Monte-Carlo weight has a non null value:

Fig. 5 .

 5 Fig. 5. The two different effects of an infinitesimal modification of the boundary on the incoming intensity at position x ∂ (π ) and direction u .

Fig. 6 .

 6 Fig. 6. The two types of discontinuities generating sources of sensitivity. The image 6 a represent the discontinuity induced by the edges connecting two plane surface. The image 6 b represent the discontinuity induced by the shading. Finally the image 6 c represent the combination of both types of discontinuities according to the radiative intensity direction.

Fig. 7 .

 7 Fig. 7. Configuration 2, on the left the square cavity geometry filled by semi-transparent medium and the deformation of the upper wall. On the right the locations of the sensitvity sources.

Fig. 8 .

 8 Fig. 8. Radiative flux on the bottom-wall J ( π) and its sensitivity ∂ π J ( π) to the height of the up-wall. The phase function in the semi-transparent medium is isotropic p U ( u i ) = 1 4 π and the probability density function of scattering extinction is p (l i ) = k s exp (-k s l i ) . Cavity dimensions are set by a variation δL of L influencing the optical thickness k e L . Single scattering albedo is uniform on the cavity ks k e = 0 . 5 . Estimations of the radiative flux and its sensitivity are obtained for 1.10 6 realizations N of the corresponding Monte-Carlo weight function.

Fig. 9 .

 9 Fig.9. Schema of the concentrating solar tower and of the heliostat field with an illustration of the heliostat size deformation. The sun surface is named S, the tower τ , the heliostats surface H and the receiver surface R .

  enlargement of a mirror symbolized in Fig. D.17 by the arrow D, can lead to

Fig. 10 .

 10 Fig. 10. Deformation of an heliostat mirror with respect to the parameter π. The parameter influences the size of the mirror area surface.

Fig. 11 .

 11 Fig. 11. Thermal power sensitivity to the size of the mirrors of an heliostat field. The medium is transparent and the deformation the width of the mirrors increases by π. The size is represented by the width of the square heliostat mirrors.

Fig. 12 .

 12 Fig. 12. Complex geometry academic configuration representing an emitting sphere inside a supershape ( Fig. 12 a) and the receptive part of the supershape surface ( 12 b).

Fig. 13 .

 13 Fig. 13. Radiative flux received by a part of the surpershape surface ( Fig. 12 b) J ( π) and its sensitivity ∂ π J ( π) to the sphere radius. The phase function in the semitransparent medium is isotropic p U ( u i ) = 1 4 π and the probability density function of scattering extinction is p (l i ) = k s exp (-k s l i ) . Cavity dimensions are set by a variation δr of r . Single scattering albedo is uniform on the surpershape k s k e = 0 . 5 . Estimations of the received radiative flux and its sensitivity are obtained for 1.10 6 realizations N of the corresponding Monte-Carlo weight function.

Declarations of Competing Interests

none The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Sensitivity boundary condition statement for a black-body surface

In this appendix we state the sensitivity boundary condition in the case of a black-body plane surface R (π ) . The

The intensity at any point inside the domain is a function of five independent variables:

and its shape derivative:

We state I p ≡ I p (x, z, u , π ) , the radiative intensity at the boundary R (π ) , in the incoming direction u . n R (π ) > 0 , such as:

With v the deformation velocity vector. We are interested on the shape derivative of the intensity at the boundary, that is the derivative of I p with respect to the geometrical parameter π :

In the case of black-bodies walls the radiative intensity at the boundary, in the incoming direction, is equivalent to Planck's distribution for every π of R (π ) . So the intensity at the boundary I p is insensitive to a variation of the shape parameter π and its derivative is null:

Furthermore we can write:

And deduced the shape boundary condition:

Example

Here a short illustration is given of the radiative boundary condition stated from the Example 5.1 The radiative boundary condition is written for u . n > 0 :

If we input this statement on the tangent gradient terms of Eq. (A.12) we will have to estimates the derivatives of Eq. (A.14) with respect to the x and z coordinates. Since the boundary condition does not depend on z variable only one derivative have to be estimated:

From Eqs. (A.15) and (A.12) we can write the general sensitivity boundary condition:

Appendix B. Treating the discontinuities from an ad hoc basis

This appendix is meant only for the identification of discontinuities sensitivity sources, related to the boundary intensity function I p . In this appendix those sources are treated separately from the previous boundary conditions statement Eq. (A.12) . and in a second part on the sphere configuration.

Heliostat mirror : reflective plane surface

We consider one heliostat surface. The The radiative intensity on a reflective heliostat surface (B.15) for x > 0, y > 0 and u . n > 0 can be stated as:

The derivation of the radiative intensity with respect to the helio- stat size parameter leads to :

2

The configuration is symmetric for the x < 0 and y < 0 values.

Appendix C. Sensitivity to the square size algorithm


Keywords: Artificial intelligence for all, Computational thinking, Free reusable online resources, MOOC. 1. Context and Objectives

Understanding Intelligently Artificial Intelligence :

a citizens' open formation.

Submitted tor http://eduai.ist.tugraz.at using this format : 2 pages abstract < 7th june , if ok 10 pages paper before 10th july. Topics : K-12 AI educational initiatives and projects, AI teaching concepts and materials, parental involvement in AI-literacy.

Fréderic Alexandre1 , Robert de Barretin 8 , Jade Becker 7 , Marie-Hélène Comte2 , Martine Courbin 2 , Sonia Cruchon3 , Aurélie Lagarrigue 2 , Bastien Masse 4 , Sophie de Quatrebarbes 5 , Julie Stein 6 , Claude Terosier 7 , and Thierry Viéville 1 .

Our whole society is and will be deeply impacted by digital science and this takes a new qualitative and quantitative turn with what is named artificial intelligence (AI). We must allow everyone to master, thus understand how all this works. This means computational thinking discovery and machine learning initiation. Directly inspired by the Finnish initiative to train 1% of its population on these subjects [START_REF]A free online introduction to artificial intelligence for non-experts[END_REF] and after our success in providing a hybrid formation on computational thinking for teachers and citizens not familiar with computer science [START_REF] Mariais | EIAH'19 Wokshop : Apprentissage de la pensée informatique de la maternelle à l'Université : retours d'expériences et passage à l'échelle[END_REF] where more 45000 persons have been reached, yielding a satisfaction level higher than 90%, we have built and now operate a citizen training in AI in the broad sense, intended to a large public beyond the school domain, with the goal of giving rise to an ubiquitary citizen university in digital science and culture [START_REF] Atlan | Apprentissage de la pensée informatique : de la formation des enseignant•e•s à la[END_REF].

Method and Production

We produced a hybrid , performative and contaminating training: Freely accessible online training course with free of charge attestation of attendance ; online exchange forum, to share and help each other, when training on these subjects; online and in real life hangouts to extend this sharing and serve local projects on these subjects, all content items being also reusable to derive and build other ressources. Hybrid : Because online ressources are designed to derive in situ activities. Performative : Because students are taught to be their own teacher, while helping each other. Performative, in the sense that an adult benefits from this training and then in his family or in his business, association, structure shares what he has discovered with his or her children or the seniors, or the other way round. Part of the training exercises is to be able to transmit what has been learned.

With respect to other offers in the field [START_REF]A free online introduction to artificial intelligence for non-experts[END_REF], the contents are less technical (e.g., machine learning programming in Python is proposed only as an extension), but not only a verbal introduction, while effective activities, including unplugged activities (i.e., without a computer using low-tech everyday tools), are proposed to really understand how-it-works. We not only talk about AI, we make people manipulate.

Contents allow us to gain a first understanding of what AI means

and what it is not, understand the principles of machine learning and the crucial role played by data, taking into account the societal issues of AI. Learners are invited to question themself, beyond preconceived ideas, and really understand the ground of AI; perform activities and manipulate AI programs to build their own vision; acquire a minimal scientific and technical culture (concepts, history of humans and ideas). As an outcome they are able to discuss the subject, question applications, choose some framework to contribute to the development (or not) of AI applications.

3.

First results and Analysis

After one month of activity [results will be updated in the final version], more than 10000 persons joined (35% female, 64% male, 1% not binary) with about 300 persons done yet, more than 1000 persons enjoying at least a part of this modular formation. They are mainly active (55% on activity, 15% retired, 11% students, 8% job researcher), mostly with university level (78% at least bachelors in any field, 10% PhD), but rather beginners in the field (57% full beginners, 40% intermediate non expert), working 10 to 20 hours in average (about 50% of the persons spend from 2 up to 5 hours per week, during about 3 weeks, while 25% spend less and 25% more), and above 90% of person having their expectation satisfied (40% fully satisfied). Detailed results [will be] available as supplementary material. On the forum, more than 1500 persons have been or are active while more than 3500 are reading the about 200 discussions, more than half being on the course contents (e.g., strong versus weak AI, symbolic versus numeric methods, societal issues, …). More than 10 online hangouts of 30 to 200 persons have been already organized by one of the 70 partners of the Class´Code project [4], during this covid-19 confinement.

Qualitatively, beyond very nice and encouraging positive feedbacks, through the constant evaluation of the course by the users, we are invited to take several improvements into account : (i) extend the existing formation with more operational tutorials, (ii) manage some technical weakness because we integrate state of the art recent external resources, some still to be consolidated, (iii) complete the existing contents to better help the learner progression. And we will, since this is a long term open course (probably more than 2 years). A few persons did not find what they expected, either because they are not interested in redoing activities with young person (on any age ;)) to better understand or because they were looking for professional formation.

On our side, we are really delighted to participate on these researcher-citizen exchanges and shares, with very interesting questions, sometimes on crucial subjects, always moderated (no troll, yet ;)). We are less happy about the fact we have mainly reached rather highly educated persons, while we target popularisation. When a person joins, the content is always considered as accessible, and the hard point is to join. We hope that in real life hangouts and also community approaches will help, and work on this aspect [START_REF] Corieri | Enjeux dans la création d'une communauté d'enseignants engagés dans l'apprentissage de l'informatique. Didapro 8 -DidaSTIC -Colloques francophones de didactique de l'informatique[END_REF][START_REF] Romero | When a Master of Sciences on EdTech becomes an International Community[END_REF].

Conclusion and Perspectives

We are proposing a "maker" approach to AI popularization between technical formation allowing to use AI at work, such as [START_REF]A free online introduction to artificial intelligence for non-experts[END_REF], and more superficial presentation of AI only in "words".

Inria, Mnémosyne team, https://team.inria.fr/mnemosyne

Inria, Learning Lab https://learninglab.inria.fr and Science Outreach Department ,

4 minutes 34 https://www.4minutes34.com

Acknowledgements The Digital Direction for Education of the French Ministry, the French Digital University Engineering and Technology , Microsoft, EducAzur , La Compagnie du Code , LINE laboratory , and all Class´Code partners .


Philippe Maury

Jean Benoit Thambo

Alice Maltret

Nicolas Combes

Sebastien Hascoet

Nicolas Derval

Magalie Ladouceur

Philippe Acar

Pascal Amedro

Frédéric Anselme

environment required for the practice of ablation in children and in congenital heart disease

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Pierre-François Marteau
email: pierre-francois.marteau@univ-ubs.fr

Warp Edit Time

2007 Distance

Valoria

Time Warp Edit Distance

Keywords: Dynamic Time Warping, Elastic Distances, Stiffness Control, Time Series matching, Timestamped Data, Event Data, Linear Segment Matching

come

I. Introduction

At the junction of symbolic edit distances [START_REF] Chen | On the marriage of Lp-norm and edit distance[END_REF] [3] [START_REF] Mäkinen | Using Edit Distance in Point-Pattern Matching[END_REF] [7] [START_REF] Wagner | The String-to-String Correction Problem[END_REF] and dynamic time warping measures [START_REF] Das | Finding Similar Time Series[END_REF] [9] [START_REF] Keogh | A simple dimensionality reduction technique for fast similarity search in large time series databases[END_REF] [START_REF] Velichko | Automatic recognition of 200 words[END_REF] we propose a family of Time Warp Edit Distance (TWED) that we denote γ λ δ , to refer to the two parameters that characterize the family, namely the gap penalty λ and the stiffness parameter γ . We first define γ λ δ , and then we prove successively that whenever 0 0 > ≥ γ λ , 1. Proposition 1: γ λ δ , is a distance metric.

2. Proposition 2: γ λ δ , is upper bounded by twice the L1 distance.

3. Proposition 3: γ λ δ , is an increasing function of γ λ and . [START_REF] Keogh | A simple dimensionality reduction technique for fast similarity search in large time series databases[END_REF]. Proposition 4: Upper-bounding the distance between a time series and its piecewise constant polygonal approximation Further details and experiments on γ λ δ , are described in [START_REF] Marteau | Time Warp Edit Distances with Stiffness Adjustment for Time Series Matching[END_REF].

II. Definitions

Let U be the set of finite discrete time series:

{ } { } Ω ∪ ∈ = + N p A U p / 1
, where Ω is the empty time series (with null length). Let p A 1 be a time series with discrete index varying between 1 and p. Let i a' be the i th sample of time series A.

) ( ) ' b ( ) B , A ( ) ' b ' a ( ) B , A ( ) ' a ( ) B , A ( Min ) B , A ( q q p , q p q p , p q p , q p , 1 1 1 1 1 1 1 1 1 1 1 1 1        → Λ Γ + → Γ + Λ → Γ + = - - - - γ λ γ λ γ λ γ λ δ δ δ δ with λ λ + = → Λ Γ + = → Γ + = Λ → Γ - - - - ) ' b , ' b ( d ) ' b ( ) ' b , ' a ( d ) ' b , ' a ( d ) ' b ' a ( ) ' a , ' a ( d ) ' a ( q q q q p q p q p p p p 1 1 1 1
where d is any distance on R k+1 . In practice, we will choose

) t , t ( d . ) b , a ( d ) ' b , ' a ( d b a Lp LP γ + =
where γ is a parameter which characterizes the stiffness of the elastic distance γ λ δ , , and λ any positive constant element in R that corresponds to a gap penalty.

The recursion is initialized setting: for some n>0. Then for all ) , (

{ } { } . ' b ' a p ,.., i ), ' a , ' a ( d ) B , A ( q ,.., j , ) ' b , ' b ( d ) B , A ( ) B , A ( i k k k i , j k k k
1 1 q p B A in U U × such that n m = , as ) B , A ( ), B , A ( q p q p , , 1 1 1 1 1 1 1 - - - γ λ δ δ γ λ and ) B , A ( q p , 1 1

P2: identity of indiscernibles

For all ) , (

1 1 q p B A in U U × , if q p B A 1 1 = then p=q and { } i i b a p i ' ' , ,..., 1 = ∈ ∀
. It is easy to show by induction on p that if

q p B A 1 1 = then 0 0 1 1 1 = ≤ ≤ ∑ = p i i i q p ) ' b , ' a ( d ) B , A ( ,γ λ δ leading to . ) B , A ( q p , 0 1 1 = γ λ δ
Now consider the backward proposition P'2:

q p q p B A ) B , A ( , 1 1 1 1 0 = ⇒ = γ λ δ . P'2 is proved by induction on m=p+q. P'2 is true for m=0. Suppose P'2 true for all { } 1 ,.. 0 - ∈ n m
for some n>0. Then for all ) , (

1 1 q p B A in U U × such that n m = and 0 1 1 = ) B , A ( q p ,γ λ δ
we have necessarily:

) ' b , ' a ( d ) ' b , ' a ( d ) B , A ( ) B , A ( q p q p q p q p , , 1 1 1 1 1 1 1 1 - - - - + + = γ λ δ δ γ λ .
We verify that the cases where

λ δ δ γ λ γ λ + + = - - ) ' b , ' b ( d ) B , A ( ) B , A ( q q q p , q p , 1 1 1 1 1 1 or λ δ δ γ λ γ λ + + = - - ) ' a , ' a ( d ) B , A ( ) B , A ( p p q p , q p , 1 1 1 1 1 1 are impossible since ) ' , ' ( 1 - q q b b d and ) ' , ' ( 1 - p p a a d
are strictly positive (the reason being that time stamps are strictly increasing).

Thus, 0

1 1 1 1 = - - ) B , A ( q p ,γ λ δ and 0 1 1 = + - - ) ' b , ' a ( d ) ' b , ' a ( d q p q p leading to 1 1 1 1 - -= q p B A and q p b a ' ' = .
Finally p=q and necessarily

q p B A 1 1 =
.

P3: Symmetry

Proof: Since the distance d on the sample space T S × is symmetric, it is easy to show that

) B , A ( q p , 1 1 γ λ δ is symmetric for all ) , ( 1 1 q p B A in U U × induction on m=p+q.

P4: Triangle inequality

For all ) , , (

1 1 1 r q p C B A in , U U U × × ). C , B ( ) B , A ( ) C , A ( r q q p r p , , , 1 1 1 1 1 1 γ λ γ λ γ λ δ δ δ + ≤ Proof:
Ω Ω + Ω Ω ≤ = Ω Ω γ λ γ λ γ λ δ δ δ 0
and the induction hypothesis holds.

(H4): Suppose P4 is true for all { } 1 ,.. 0 - ∈ n m for some n>0. Let ). C , B ( ) B , A ( r q q p , , 1 1 1 1 γ λ γ λ δ δ + = Σ
Then for all ) , , (

1 1 1 r q p C B A in U U U × ×
such that m=n, we have basically 9 different cases to explore for the decomposition of

) B , A ( q p , 1 1 γ λ and ) C , B ( r q , 1 1 γ λ δ : 1 st Case: if      → Γ + = → Γ + = - - - - ) ' c ' b ( ) C , B ( ) C , B ( ) ' b ' a ( ) B , A ( ) B , A ( r q r q , r q , q p q p , q p , 1 1 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then
. inequality r triangula the verifies and applies, ) 4 ( since

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 d H ) C , A ( ) ' c , ' a ( d ) ' c , ' a ( d ) C , A ( ) ' c , ' b ( d ) ' c , ' b ( d ) ' b , ' a ( d ) ' b , ' a ( d ) C , B ( ) B , A ( r p , def r p r p r p , r q r q q p q p r q , q p , γ λ γ λ γ λ γ λ δ δ δ δ ≥ + + ≥ + + + + + ≥ Σ - - - - - - - - - - - - 2 nd Case: if      → Γ + = → Λ Γ + = - - - ) ' c ' b ( ) C , B ( ) C , B ) ' b ( ) B , A ( ) B , A ( r q r q , r q , q q p , q p , 1 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then ) C , A ( , ) ' c , ' c ( d ) C , A ( d , ) ' c , ' c ( d ) C , B ( ) B , A ( ) ' c , ' b ( d ) ' c , ' b ( d ) C , B ( ) ' b , ' b ( d ) B , A ( r p , def r r r p , r r r q , q p , r q r q r q , q q q p , applies (H4) ) inequality r triangula the verifies ( γ λ γ λ γ λ γ λ γ λ γ λ δ λ δ λ δ δ δ λ δ ≥ + + ≥ + + + ≥ + + + + + = Σ - - - - - - - - - - - - 3 rd Case: if      → Γ + = Λ → Γ + = - - - ) ' c ' b ( ) C , B ( ) C , B ( ) ' a ( ) B , A ( ) B , A ( r q r q , r q , p q p , q p , 1 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then ) C , B ( ) ' a , ' a ( d ) B , A ( r q , p p q p , 1 1 1 1 1 1 γ λ γ λ δ λ δ + + + = Σ - - ) C , A ( H ) ' a , ' a ( d ) C , A ( r p , def p p r p , 1 1 1 1 1 1 applies ) 4 ( γ λ γ λ δ λ δ ≥ + + ≥ Σ - - 4 th Case: if      → Λ Γ + = → Λ Γ + = - - ) ' c ( ) C , B ( ) C , B ( ) ' b ( ) B , A ( ) B , A ( r r q , r q , q q p , q p , 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then λ δ δ γ λ γ λ + + + = Σ - - ) ' c , ' c ( d ) C , B ( ) B , A ( r r r q , q p , 1 1 1 1 1 1 ) C , A ( H ) ' c , ' c ( d ) C , A ( r p , def r r r p , 1 1 1 1 1 1 applies ) 4 ( γ λ γ λ δ λ δ ≥ + + ≥ Σ - - 5 th Case: if      → Λ Γ + = → Γ + = - - - ) ' c ( ) C , B ( ) C , B ( ) ' b ' a ( ) B , A ( ) B , A ( r r q , r q , q p q p , q p , 1 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then ) C , A ( ) H ( ) ' c , ' c ( d ) C , A ( ) ' c , ' c ( d ) C , B ( ) B , A ( r p , def r r r p , r r r q , q p , 1 1 1 1 1 1 1 1 1 1 1 1 applies 4 γ λ γ λ γ λ γ λ δ λ δ λ δ δ ≥ + + + ≥ + + + = Σ - - - - 6 th Case: if      → Λ Γ + = Λ → Γ + = - - ) ' c ( ) C , B ( ) C , B ( ) ' a ( ) B , A ( ) B , A ( r r q , r q , p q p , q p , 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then ) C , A ( ) H ( ) ' c , ' c ( d ) C , A ( ) ' c , ' c ( d ) C , B ( ) B , A ( r p , def r r r p , r r r q , q p , 1 1 1 1 1 1 1 1 1 1 1 1 applies 4 γ λ γ λ γ λ γ λ δ λ δ λ δ δ ≥ + + + ≥ + + + = Σ - - - - 7 th Case: if      Λ → Γ + = → Γ + = - - - ) ' b ( ) C , B ( ) C , B ( ) ' b ' a ( ) B , A ( ) B , A ( q r q , r q , q p q p , q p , 1 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then . applies 4 ) inequality triangle the satisfies ( 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ) H ( ) C , A ( ) ' a , ' a ( d ) C , A ( d ) ' a , ' a ( d ) C , B ( ) B , A ( ) ' b , ' b ( d ) C , B ( ) ' b , ' a ( d ) ' b , ' a ( d ) B , A ( r p , def p p r p , p p r q , q p , q q r q , q p q p q p , γ λ γ λ γ λ γ λ γ λ γ λ δ λ δ λ δ δ λ δ δ ≥ + + ≥ + + + ≥ + + + + + = Σ - - - - - - - - - - - - 8 th Case: if      Λ → Γ + = Λ → Γ + = - - ) ' b ( ) C , B ( ) C , B ( ) ' a ( ) B , A ( ) B , A ( q r q , r q , p q p , q p , 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then: ) C , A ( ) H ( , ) ' a , ' a ( d ) C , A ( ) C , B ( ) ' a , ' a ( d ) B , A ( r p , def p p r p , r q , p p q p , 1 1 1 1 1 1 1 1 1 1 1 1 applies 4 γ λ γ λ γ λ γ λ δ λ δ δ λ δ ≥ + + ≥ + + + = Σ - - - - 9 th Case: if      Λ → Γ + = → Λ Γ + = - - ) ' b ( ) C , B ( ) C , B ( ) ' b ( ) B , A ( ) B , A ( q r q , r q , q q p , q p , 1 1 1 1 1 1 1 1 1 1 γ λ γ λ γ λ γ λ δ δ δ δ , then: . applies 4 1 1 1 1 1 1 1 1 ) H ( ) C , A ( ) C , B ( ) B , A ( r p , r q , q p , γ λ γ λ γ λ δ δ δ ≥ + ≥ Σ - - So property P4 holds for all m in { } n ,.. 0
. By induction P4 holds for all m in N and so P4 holds for all ) , , (

1 1 1 r q p C B A in U U U × × .
IV. Proposition 2: γ λ δ , is upper bounded by twice the Lp distance.

Proposition 2: ), Y , X ( D ) Y , X ( U Y , X , LP , ⋅ ≤ ∈ ∀ > ≥ ∀ 2 0 0 2 γ λ δ γ λ
whenever X and Y have the same length.

Proof: let us consider the sequence of editing operations consisting in m match operations, where m is the length of the X and Y. This sequence has a cost equal to twice the Lp-distance between the two time series X and Y. Since γ λ δ , is equal to the cost of the optimal sequence of edit operations, the result follows.

V. Proposition 3:

γ λ δ , is an increasing function of γ λ and Proposition 3:

) , ( ) , ( , ' ' 0 , 0 ' , ' , 2 Y X Y X U Y X γ λ γ λ δ δ γ γ λ λ γ λ ≤ ∈ ∀ ≥ ∀ ≥ ∀ > ≥ ∀ Proof:

VI. Proposition 4: Upper-bounding the distance between a time series and its piecewise constant polygonal approximation.

We define r p A , 1 as a Piece Wise Constant Approximation (PWCA) of time series p A 1 containing 0 1 ≥ r constant segments and p samples. This approximation can be obtained using any kind of solution (from heuristic to optimal solutions), let say the optimal solution similar to the one proposed in [START_REF] Perez | Optimum polygonal approximation of digitized curves[END_REF]. Proposition 4:

) r p ( T ) r p ( ) X , X ( U X [, p ; [ r , , r r , p , p - ⋅ ∆ ⋅ + - ⋅ ≤ ∈ ∀ ∈ ∀ > ≥ ∀ 2 1 0 0 1 1 1 γ λ δ γ λ γ λ

VII Conclusion

We have proposed a family of time warp edit distances for time series matching. This family involves two parameters: the stiffness parameter that controls the elasticity of the distance and the gap penalty that is part of the cost involved in the insert or delete operations.

We have shown that the proposed measure:

-is a metric on the set of discrete and finite time series.

-is upper bounded by twice the Lp-distance.

-is an increasing function of the gap penalty and the stiffness parameter.

Further more, the distance between a time series and its piecewise constant approximation is upper bounded by an expression that only depends on the lengths of the times series, the number of segments of its approximation and the two parameters γ λ and .

  m=p+q. Non-negativity of γ λ δ , is proved by induction on m. P1 is true for m=0 by definition of γ λ δ , and the induction hypothesis holds. Suppose P1 is true for all

  Let us consider one of the optimal sequences of editing operations evaluated with the elementary operation costs we get a cost for this sequence that is lower than ) . The result follows.

  is the time difference average between two successive samples inside the piecewise constant segments of the approximation.The proof for this proposition is straightforward: let us consider the sequence of operations consisting in r match operations for the end extremities of the piecewise constant segments and (p-r) delete operations for the set of samples in r p X , 1 that are not end extremities of the piecewise constant segments. In this sequence each match operation has in average the cost delete operation has a λ fixed penalty and a penalty proportional to

  We will prove P4 by induction on m=p+q+r.

	P4 is true for m=0 since	,	(	,	)	,	(	,	)	,	(	,	)