Modeling and optimizing a road de-icing device by a nonlinear heating

In order to design a road de-icing device by heating, we consider in the one dimensional setting the optimal control of a parabolic equation with a nonlinear boundary condition of the Stefan–Boltzmann type. Both the punctual control and the corresponding state are subjected to a unilateral constraint. This control problem models the heating of a road during a winter period to keep the road surface temperature above a given threshold. The one-dimensional modeling used in this work is a first step of the modeling of a road heating device through the circulation of a coolant in a porous layer of the road. We first prove, under realistic physical assumptions, the well-posedness of the direct problem and the optimal control problem. We then perform some numerical experiments using real data obtained from experimental measurements. This model and the corresponding numerical results allow to quantify the minimal energy to be provided to keep the road surface without frost or snow.


Introduction
De-icing a road pavement is an important issue in many countries subjected to winter weather conditions that have a strong impact on road maintenance and road safety. The use of salt spreaded to ensure the de-icing of pavements can affect the environment close to the road. Also, devices have been implemented to heat the road: electric heating, infrared lamps above the road surface, circulation of a heat transfer fluid in pipes inserted in the road [20,21,22,13,19,10,11]. The concept of heated pavement is not new. One of the most notable is the Serso project built in Switzerland, with a surface of 1,300 m 2 : a heat transfer fluid, which draws its energy in the ground, circulates in pipes inserted in the wearing layer pavement [10]. Due to the energy issue being an important concern in all the world, the concept of a positive energy road has emerged since few years. By way of illustration, considering a mean solar energy on the French territory of 1,400 kWh/m 2 /year, the French national road network captures 196 billion kWh/year, seven times more energy than needed to ensure for example its surface de-icing by heating.
Under the impulse of the European project R5G (Road of 5 th generation), recent research [3,17] has been undertaken on heating the road by circulating a heat transfer fluid within a porous layer of the road. In these works, the pavement structure studied is composed of three asphalt layers where the central one is a highly porous draining asphalt through which circulates a fluid (water) via gravitational flow under slant effect. A demonstrator is developed by Cerema (Centre for expertise and engineering on risks, urban and country planning, environment and mobility), EATP (French school of public works) and IFSTTAR (French institute of science and technology for transport, development and networks). This is an experimental roadway constructed in 2014 at Egletons, France. The road with 50 m in length and 4 m in width is composed of three layers: a wearing course layer of semi-phaneritic asphalt concrete 0.06 m thick; a bonding course layer of 0/14 porous asphalt 0.08 m thick; and a base layer of asphalt concrete with a thickness of 0.05 m. The transversal slope is around 2% with no longitudinal slope. The demonstrator consists of two parts: a control road and a road in which a fluid (water) circulates through the porous asphalt layer (see Figure 1). The second one is equipped with two tanks for the supply and recovery of fluid circulating in the bonding course layer. A pump is used to remove fluid from the downstream tank and reinject it into the upstream tank. Fluid circulation is maintained in the porous asphalt layer by a watertight seal between it and the underlying layer (see Figure 2). The weather is known on the site thanks to meteorological sensors delivering time-real data on humidity and temperature of the air, radiation fluxes (solar and infrared), wind speed, rainfall and snow. Temperature sensors are inserted in the road from surface until one meter deep in each part of the road (the control one and the one circulated by coolant). To address the energy issue in an experimental way, a heat storage in the granitic soil will be studied by implanting 3 vertical geothermal probes of depth 50m, allowing to take energy from the ground to heat the road in winter, and to store in the ground the heat coming from the road in summer. The demonstrator was presented in a French TV show in 2016 [18]. Whatever the technique used (electrical heating, pipes, porous asphalt), the design of the devices needs to quantify the minimal energy to be provided to keep the road surface without frost or snow, for a given winter meteorological scenario. In order to address this question, we consider in this work the optimal control of a generic heating road model. We are interested in heating during the time interval the road by a punctual heating source inserted in the road, which is modeled according to its vertical dimension [8,5,2,12]. The optimization of the heating through the punctual source leads to an optimal control associated to a one dimensional heat type equation. The underlying functional to minimize represents the energy expended by the positive source over the time period. On a mathematical viewpoint, we are then face to a optimal control problem, subjected to constraints both on the source control and on the temperature state, this latter being solution of a partial differential equation with a nonlinear boundary condition of the Stefan-Boltzmann type.
This paper is organized as follows. In Section 2, we describe respectively the bi-dimensional road heating model and its longitudinal one dimensional reduction leading to the aforementioned non linear optimal control problem. Section 3 and Section 4 then study the mathematical properties of the one dimensional model and its associated control problem respectively. Section 5 describes the finite dimensional approximation use to solve the control problem then discuss some numerical experiments. Sections 6 provides three extensions and perspectives and 7 concludes.

Modelling of the road heating and optimal control
We first describe the modeling of the heating thanks to the circulation of a coolant in a bonding porous layer of the road described in Figure 2. Space variables are x along the sub-horizontal transversal axis of the pavement with slant angle β and y along the upwards sub-vertical axis, perpendicular to x: we refer to Figure 3 which depicts a transversal two dimensional view of the road. The road is assumed to   de  no  te  sp  eci  fic  he  at,  the  rm  al  con  du  cti  vit  y  an  d  po  ros  ity  of  lay  er  i,  sp  eci  fic  he  at  an  d  the  rm  al  con  du  cti  vit  y  of  the  flu  id,  Da  rcy  flu  id  vel  oc  ity  alo  ng  x  ax  is  (he  re  su  pp  ose  d  to  be  un  ifo  rm  wi d  in  Fig  ur  e  4,  bo  un  da  ry  con  dit  ion  s  are  ho  mo  gen  eou  s  Ne  um  an  n  ex  cep  t  for  the  up  str  eam  con  dit  ion  of  po  rou  s  asp  ha  lt  lay  er Fo  r  the  fir  st  on  e,  the  inj  ect  ion  tem  pe  rat  ur  e  of  the  flu  id  is  im  po  sed : Fo r 0 ≤ x ≤ L an d 0 ≤ y ≤ h, on e ha s: , v an d K den ote spe cifi c hea t, the rm al con du cti vit y an d po ros ity of lay er i, spe cifi c hea t an d the rm al con du cti vit y of the flu id, Da rcy flu id vel oc-ity alo ng x ax is (he re sup po sed to be un ifo rm wit h x an d y) an d hy dra uli c con du cti vit y of the po rou s asp ha lt, res pe cti vel y. H 1 an d H 2 rep res ent hy dra uli c hea ds im po sed upstr eam an d do wn str eam of flu id cir cul ati ng in po rou s dra ini ng asp ha lt lay er. We stu dy flo w cir cul ati on in a sat ura ted env iro nm ent , wh ich cor res po nd s to H 1 −H 2 ≥ βL . T de-no tes tem pe rat ure fiel d in the pa vem ent . We ass um e equ ali ty of tem pe rat ure s be tw een the flu id an d dra ini ng asp ha lt at an y po int in spa ce. Th e mo del is dis cre tiz ed usi ng a fin ite diff ere nce me tho d. At the int erf ace s be tw een lay ers , con dit ion s of con tin uit y of tem pe rat ure an d the rm al flo w are im po sed . In the fol low ing , ε, σ, T s , R atm , α, R g an d H v den ote em iss ivi ty, Ste fan -B olt zm an n con sta nt (5. 67 × 10 −8 W /m 2 K 4 ), sur fac e tem pe rat ure ( o C) , atm osp her ic rad iat ion (W /m 2 ), alb edo , glo ba l rad iat ion (W /m 2 ) an d con vec tio n flu x (W /m 2 ), res pe cti vel y. As me nti on ned in Fig ure 4, bo un da ry con dit ion s are ho mo gen eou s Ne um an n exc ept for the up str eam con dit ion of po rou s asp ha lt lay er (x = 0, e 1 ≤ y ≤ e 1 + e 2 ) an d the roa d sur fac e con dit ion (y = 0). Fo r the firs t on e, the inj ect ion tem pe rat ure of the flu id is im po sed : (2) Figure 3: Scheme of pavement structure with its limit conditions (θ f , θ a are the injection temperature of the fluid and the air temperature respectively).
have no longitudinal slant and to be infinite in its third dimension. h e and L denote the height of the road structure and its length, respectively. The hydraulic regime is assumed stationary with hydraulic parameters independent of temperature θ, expressed in Kelvin. Denoting by 1 ≤ i ≤ 4 the indices of the road layers, the thermo-hydraulic model is as follows. For t > 0, 0 ≤ x ≤ L and 0 ≤ y ≤ h e : and K denote specific heat, thermal conductivity and porosity of layer i, specific heat and thermal conductivity of the fluid, Darcy fluid velocity along x axis (here supposed to be uniform with x and y) and hydraulic conductivity of the porous asphalt, respectively. H 1 and H 2 represent hydraulic heads imposed upstream and downstream of fluid circulating in porous draining asphalt layer.
The assumption of a saturated fluid circulation corresponds to H 1 − H 2 ≥ βL. As mentioned in Figure 3, boundary conditions for problem (2.1) are homogeneous Neumann except for the upstream condition of porous asphalt layer (x = 0, e 1 ≤ y ≤ e 1 + e 2 ) and the road surface condition (y = 0). For the first one, the injection temperature of the fluid is imposed : The second one, that is the road surface boundary condition expresses the energy balance between road and atmosphere where the following notations are used : According to [2], the convection coefficient is defined by where the following notations are used : Cp a : thermal capacity (J/kg.K) of the air, ρ a : density of the air (kg/m 3 ), V wind : wind velocity (m/s), C d , C d1 : two convection coefficients (-).
As mentioned in the introduction, we use and study in this work, as a first step, a one dimensional reduced model obtained by fixing the sub-horizontal axis x. The injection temperature term θ f supported on the boundary {0} × (e 1 , e 1 + e 2 ) is replaced by a punctual heating source q inserted in the road (localized at y = y 0 ∈ (0, h e )). For any T > 0, we denote Q T := (0, h e ) × (0, T ). The evolution of the temperature along the road, now modeled according to its vertical dimension (as explained and used in [8,5,2]) is now modeled as follows: Here and in the sequel, θ t and θ y stand for the partial derivative of the function θ = θ(y, t) with respect to the space variable y and the time t respectively. The time positive functions f 1 and f 2 which appears in the nonlinear boundary conditions at y = 0 are defined by: The optimization of the heating through the punctual source q localized at y = y 0 ∈ (0, h e ) leads to the following optimal control problem : for some positive value θ, which represents the minimal temperature required. The functional J represents the energy expended by the power q over the period [0, T ]. (2.7) is thus an optimal control problem subjected to contraints inequality both on the control variable and the temperature state, this latter being solution of a partial differential equation with a nonlinear boundary condition of the Stefan-Boltzmann type. The constraint condition on the state, i.e. θ(0, t) ≥ θ, is not frequent in the optimal control literature as it involves technical developments and is reminiscent of Signorini type conditions (see [9] and also section 6.1). We mention [16] which addresses the controllability of the one dimensional linear heat equation subjected to the boundary condition θ y (0, t) = θ 4 (0, t) + u(t), u being the control. We also mention [1] which addresses the controllability of a linear string submitted to a unilateral obstacle at one extremity.
(3.10) (, ) H denotes the inner product over H while the form a : In order to simplify the notation, we introduce, for almost all t ∈ (0, T ) and g ∈ V , the element Φ(g) = g(0) 3 |g(0)|δ 0 which be can identified with an element in V so that (3.11) Using that φ := ( φ y 2 H + φ(0) 2 ) 1/2 defines a norm equivalent to φ V , the following holds true.
Let us then introduce the following regularity assumptions : Theorem 3.1 Assume the hypothesis (H). Assume moreover that the control q satisfies q(0) = 0 and that the initial condition θ 0 satisfies the compatibility condition at the point y = 0 and y = h e respectively. Then, there exists a unique solution θ of (3.9) and (3.10) such that θ, θ t ∈ L 2 (0, T, V ) ∩ L ∞ (0, T, H).
Moreover, there exists a constant and a constant In particular, this implies that the solution belongs to Proof -The proof follows the arguments developed in [9], Chapter 1, Section 5, based on the monotony of Φ, that is (3.14) The dependance of the emissivity function ε in (3.9) with respect to the time variable requires however additional developments.
Step 1 -Galerkin approximation. We introduce an orthonormal basis {w k } k≥1 of the separable space V and define, for any integer m, by V m the subspace of V generated by {w k } 1≤k≤m . Moreover, for simplicity, we choose w 1 in order that θ 0 = w 1 . Then, we define the solution θ m (t) ∈ V m of the following finite differential system: Consequently, the function θ m : [0, t m ] → V is given by θ m (t) = m k=1 g k (t)w k where the time function {g k } 1≤k≤m solves a differential system for some time t m . The following a priori estimates show that t m > 0 is independent of m.
Step 2 -First a priori estimates. The sum over k = 1, . . . , m of (3.15) first leads to leading, in view of (3.14), to Therefore, Lemma 3.1 implies Taking α 2 small enough so that α(t) − C L α2 2 ≥ 0, we obtain that Step 3 -Second a priori estimates. Taking t = 0 in (3.15), we obtain from which we deduce, following [9], the bound Then, by differentiating (3.15) with respect to time, we get, for and then to with α defined in (3.12). To estimate the remaining boundary term, we multiply (3.15) by g k,t and sum over k leading to the equality and therefore to Using Lemma 3.1, we write where β is defined in (3.12). Using again several times the Young inequality, we obtain the inequality for any α 1 , α 2 , α 3 , α 4 > 0 small enough, from which we deduce that A Gronwall type lemma then allows to deduce that, for all t ∈ (0, T ), Uniform estimates (3.17) and (3.19) then imply that √ cθ m,t is uniformly bounded in L ∞ (0, T, H) with respect to the parameter m. Inequality (3.21) then implies that θ m,t is uniformly bounded in L 2 (0, T, V ) as well, and that there exists a constant C > 0 such that Step 4 -Limit as m → ∞. From (3.22) and (3.19), θ m is bounded uniformly in L ∞ (0, T, V ) and therefore, from (3.15) The previous convergences imply notably that θ µ (0) θ 0 as µ → ∞ in V separable so that θ 0 = θ(0). Writing (3.15) for the subsequence θ µ and using (3.23), we obtain that θ and Q satisfy the equality (3.24) It remains to pass to the limit in the nonlinear term in order to show the equality To do that, we use the equality (3.15): precisely, let φ ∈ L 2 (0, T, V ) and let The monotony of Φ implies that U µ ≥ 0. Moreover, the equality − H and the weak star convergence of θ µ in L ∞ (0, T, V ) (and so in L ∞ (0, T, H)) allows to take the lim sup of U µ . Precisely, from the upper semi- H for all t ∈ (0, T ) and therefore, On the other hand, the integration over (0, T ) of (3.24 The sum of the two previous relations implies that Step 5 -Uniqueness. By contradiction, let θ, θ be two distinct solutions of (3.9) and let w = θ − θ .
The arguments of step 5 allow to show that if θ 1 and θ 2 are solution of (3.9)-(3.10) associated to q 1 and q 2 respectively, the other data being equals, then there exists a constant C > 0 such that Remark 3.2 A similar well-posedness result holds true if the Neumann boundary condition at y = h e is replaced by a non homogenous Dirichlet boundary θ(h e , ·) = θ d ∈ R. This latter, which is physically relevant if the height h e is large enough, will be used in the numerical section 5.
Moreover, the solution enjoys the following comparison principle.
Proposition 3.1 Assume the hypothesis of Theorem 3.1. Let θ andθ the solutions of (3.9)-(3.10) associated to the pair (q, θ 0 ) and (q, The right hand side is positive, while, in view of the monotony of Φ, the term Φ( and using that a(t, m(t), m + (t)) = a(t, m + (t), m + (t)) and that (cm 2 In particular, since the time function f 1 is non negative, we deduce the following property. It results that if the source term and if the initial condition are non negative, then the boundary value problem (3.8) coincides with the initial one (2.4).
4 Analysis of the optimal control problem 4.1 Well-posedness of the optimal control problem Let α > 0 and θ > 0. We assume that the initial condition θ 0 is non negative a.e. in [0, h e ]. In view of the regularity assumption on q in Theorem 3.1, we introduce the following Tychonov regularization of the optimal control problem (2.7) : where the constraint set is given by Consequently, if the control function q belongs to C, then from Corollary 3.1, the weak formulation of (2.4) coincides with the weak formulation of (3.8) and is well-posed. We relax the second inequality constraint from C and introduce the equivalent extremal problem: where ψ K is the indicator function of K, that is it ψ K (q) = 0 if q ∈ K and ψ K (q) = +∞ else with K = {q ∈ H 1 (0, T ) s.t. Lemma 4.1 Let us assume that θ 0 ≥ θ on (0, h e ). If one of the following hypothesis holds true, then K is not empty.
2 According to the physical intuition, if the heat source is located on the surface, then the temperature can be maintained as large as requested. In view of the continuity of the temperature θ(y, t) with respect to the longitudinal axis y, the same result holds true if y 0 is negative but closed enough to 0. Actually, we believe that the same result holds true if the heat source, located on any y 0 < 0, is large enough, the temperature θ being a continuous and monotonous function of q, see Proposition 3.1.
Lemma 4.2 K is a closed convex subset of H 1 (0, T ) and so ψ K is convex and semi-lower continuous.
2 It results that the equivalent optimal control problem (4.26) is well-posed as well. From a practical viewpoint, it is convenient to address this problem with a penalty method. For any parameter > 0, we introduce the penalized extremal problem : where a − = min(a, 0) denotes the negative part of any real a and D defined in (4.27).
Theorem 4.2 For all > 0, Problem (P ) admits a unique solution q . Moreover, q strongly converges in H 1 (0, T ) as → 0 to q, the solution of (4.26).
Proof-The functional J α, is strictly convex, satisfies the property lim q H 1 (0,T ) →+∞ J α, (q) = +∞ and is continuous over H 1 (0, T ) in view of remark (3.1). Moreover, the set D is a closed convex set of H 1 (0, T ) which gives the existence and uniqueness of a solution q of Problem (P ). The strong convergence of q with respect to is the consequence of the strict convexity of J α, . Moreover, in view of Theorem (3.1), this strong convergence implies the convergence of the corresponding solution θ = θ(q ) to the solution θ(q) in C([0, T ], V ).

Optimality system for the penalized extremal problem
We derive in this section the optimality condition associated to the extremal problem (P ). This allows to determine the first order variation of the cost J α, and define a minimizing sequence. Let Theorem 4.3 For any α > 0, > 0, the functional J α, is Gâteaux differentiable on the set D and its derivative at q ∈ D in the admissible direction q (i.e. q ∈ H 1 0,0 (0, T ) such that q + ηq ∈ D for all η = 0 small) is given by where p solves the adjoint problem and θ q solves (2.4). As a consequence, the unique minimizer q in the convex set D of the convex functional J α, is characterized by the optimality condition < J α, (q ), q − q > ≥ 0, ∀q ∈ D.
Remark 4.2 (4.38) is equivalent to non local differential equation Remark that actually we may drop the L 1 -norm in the left hand side: since H 1 (0, T ) ⊂ L 1 (0, T ), the resulting sequence is still in D. Moreover, we may impose the Dirichlet assumptionq n (T ) = 0; (4.38) is unchanged except that the space for the test functions q is then H 1 0 (0, T ).
5 Approximation of the optimal control problem and experiments

Numerical approximation
The approximation of the variational formulations (3.9) and (4.35) is performed using a finite element approximation with respect to the space variable y and a finite difference approximation with respect to the time variable t. For convenience, we replace the Neumann boundary condition at y = h e by a Dirichlet condition (see Remark 3.2). Let [y i , y i+1 ]. We note h = max i |y i+1 − y i |. We then introduce the following conformal finite element approximation V h of V : where P 3 denotes the space of polynomial functions of order 3. We also consider, for some The weak formulation associated to (2.4) (see (3.9-3.10)) is then approximated as follows: where π h : V → V h is the projection operator over V h . Similarly, let N t be a positive integer and (t n ) n=1,·,Nt a uniform subdivision of the time interval [0, T ] such that t 0 = 0, t Nt = T and t n = n∆t for all n and [0, T ] = ∪ Nt−1 n=0 [t n , t n+1 ]. We note by (θ n h ) an approximation of θ h (·, t n ) the solution of the following implicit Euler type scheme: In particular, the nonlinear term θ 4 h (t n+1 ) is approximated as follows : In the sequel, we define Θ h = {θ n h } n=0,··· ,Nt ∈ (V h ) Nt+1 . The extremal problem (P ) is then approximated by the following one : .
Here, Π ∆t (Θ h (0))(t) is the piecewise affine function such that Π ∆t (Θ h (0))(t n ) = θ n h (0) for all n where (θ n h ) solves (5.41) and The extremal point of the functional J α, in D ∆t is determined by using the gradient projection algorithm described in the previous section. In particular, the adjoint problem (4.35-4.36) is approximated as follows:

Presentation of the experimental data and initial condition
We use real data obtained from measurements on the French highway A75 in Cantal (1100 meter altitude) from october 2009 to march 2010. Measurements are made each hour and allow to compute the time functions f 1 and f 2 defined in (2.6). Figure 4 depicts these functions. For the other parameters, we use numerical values obtained from experimental validations described in [2]: precisely, the albedo of the road surface is A = 0.08 (used in f 1 ), while for the emissivity function of the road, we use the constant value ε(t) = 0.92. Even if measurements of ε show a time variability and a sensitivity of the road surface temperature θ(0, ·) w.r.t. ε (see Figure 5), the lack of emissivity measurements on a long period requires to identify a constant value ε = 0.92 leading however to a good agreement between simulated and measured temperatures [2]. Remark however that Theorem 3.1 is valid for time-dependent function.
As described in Figure 3, the road is composed of 4 layers of depth e 1 = 0.06m, e 2 = 0.08m, e 3 = 0.05m and e 4 =h e − e 1 + e 2 + e 3 = 14.81m respectively. The total height of the road structure is h e = 15m. The specific heat function c and the thermal conductivity function k are constant on each layer and takes the following values [2]: Here, we have replaced the homogeneous Neumann boundary condition θ 0,y (h e ) = 0 by a Dirichlet condition θ(L) = θ d ; the reason is that the temperature, at the height h e = 15m under the road surface, is time independent and close to a value equal to 15 degree celsius (equivalently 273.15 + 15 Kelvin). The resolution of the non linear boundary value problem with θ d = 288.15, using a Newton type algorithm and a finite element discretization as above, leads to a continuous function, affine on each layers, and increasing from the value θ 0 (y = 0) ≈ K + 6.29 Kelvin at the road surface to θ 0 (y = h e ) = θ d Kelvin. Moreover, by definition, this solution θ 0 satisfies the assumption (k(θ 0 ) y ) y ∈ V (see H), that is the jump [k(θ 0 ) y · ν] is equal to zero at y = e 1 , e 1 + e 2 and at y = e 1 + e 2 + e 3 . Eventually, we observe that the value θ 0 (y = 0) ≈ 279.43 Kelvin as the initial temperature at the road surface is in agreement with typical measurements in october (we recall that our study is based on measurements from october 2009 to march 2010).

Numerical experiments : Optimal control minimizing J α,
We now discuss some numerical experiments associated to the optimal problem (P ): very precisely, we minimize over D the functional .

(5.44)
It is necessary to adjust the constants in front of the norms in term of the time interval T , which is, expressed in seconds, large: precisely T = 15494400 seconds (equivalently 4304 hours). We use for that the inequality q L 1 (0,T ) ≤ √ T q L 2 (0,T ) and the Wirtinger inequality q L 2 (0,T ) ≤ T /(2π) q t L 2 (0,T ) (taking q(0) = q(T ) = 0). For the same reasons, coefficients must be properly chosen in the computation of the descent direction (see (4.38)): we use here Similarly, the penalty parameter is adjust so that the optimal control leads to values of the same order for the norm (θ q (0, ·) − θ) − L 2 (0,T ) , independently of the value of the parameter α, the value of the height y 0 (punctual support of the control q) and of the minimal temperature required θ.
The time discretization parameter ∆t is given by the data: ∆t=one hour over a period of 4304 hours. Scaled to a time interval of length one, ∆t = 1/4304 ≈ 2.32 × 10 −4 . Concerning the approximation with respect to y, we use a non uniform discretization of the interval [0, h e ]: precisely, we use 50 finite elements in each layers of height 0.06, 0.08, 0.05 and 14.81 meters respectively. In the sequel, all iterative algorithms used to approximate the optimal controls are initialized with a zero source q 0 = 0 and stopped at iterate k with the criterion |J(q k ) − J(q k−1 )|/J(q 0 ) ≤ 10 −10 , where J is the minimized function and q k is the k th iterate of the control.
In order to analyse the specific effect of the L 1 norm minimization in (5.44), we start by the standard minimization of the L 2 (0, T ) norm : inf .

(5.47)
We gather in Table 1 some numerical norms of the optimal control q of (5.46) with respect to the penalty parameter . The punctual control q is located on the middle of the porous layer, i.e. at y 0 = 0.1m. The minimal temperature required along the six months period is θ = 275.15 Kelvin: we therefore consider a slightly larger value than 273.15 (corresponding the 0 degree Celsius) in order to account for the possible errors in the experimental measurements and models. Remark that the initial condition θ 0 , solution of (5.43), satisfies this property. As expected, q L 2 (0,T ) , q L 1 (0,T ) and q L ∞ (0,T ) increase and (θ(0, ·) − θ) − L ∞ (0,T ) decreases when decreases.  Table 1: Numerical norms of the optimal control of problem (5.46) with respect to the parameter ∈ {10 −6 , 10 −4 } -y 0 = 0.1. Table 2 collects some numerical norms of the optimal control q (minimizer of the functional J α, ) with respect to the regularizing parameter α. The penalty parameter is chosen so that, for each α, the L 2 -norm (θ(0, ·) − θ) − L 2 (0,T ) takes approximatively the value 5 × 10 1 , leading to a violation of the constraint θ(0, t) ≥ K + 2 of the order 10 −1 degree (see (θ(0, ·) − θ) − L ∞ (0,T ) in Table 2).  Table 2: Numerical norms of the optimal control minimizing (5.44) with respect to the parameter α ∈ {10 −7 , 10 −3 } -y 0 = 0.1.
As we can see in Figure 6, the L 1 -minimization (smaller α) leads to a sparse control with a higher L ∞ norm (see [6]), but with a smaller L 1 norm. Moreover higher is α, more regular is the control, smaller is the L ∞ norm and higher is the L 1 norm (see Table 2).
In order to avoid large punctual values for the optimal L 1 control, which may be not physical or unrealistic at the practical level, we constrain it to take its values in a given interval [0, λ]. Results in Section 4 may be easily adapted to this new situation. We notice the effect of such boundedness on the sparsity of the control, which tends to decrease when λ decreases, as we can see especially for λ = 200. However, in this latter case, q L 1 (0,T ) = 5.78 × 10 8 , which is much higher than 4.63 × 10 8 or 4.46 × 10 8 corresponding respectively to λ = 285 and λ = +∞. Moreover, the value λ is not large enough to satisfy the constraint at the surface of the road, since we obtain (θ(0, ·) − θ) − L ∞ (0,T ) of order 2.0, much larger to that corresponding to λ = 285 and λ = +∞ (less than 0.2), as we can see in Figure 8. A minimal maximal punctual value of the source term is then necessary to maintain the road out of frost or snow (see Section 5.4 on bang-bang controls).
It is also interesting to study the local effect of control. As can be seen in Figure 9, the more the control is regular, the less it reacts instantaneously with respect to the uncontrolled surface temperature θ q=0 (0, ·).

Bang-bang control
We are now looking for so-called bang-bang controls: there are more convenient on a practical viewpoint as they take only a finite number of distinct values. Precisely, we assume that controls q take the form and minimize the amplitude λ of the control q, assume piecewise constant. A volume constraint through the parameter L is introduced here to prevent the optimal control to be constant q(t) = λ (i.e. s(t) = 1) for all t. Adapting [15], we introduce two penalty parameters and 1 for the constraints s ∈ {0, 1} and θ(0, t) − θ > 0, t ∈ (0, T ) respectively and a scalar Lagrange multiplier µ for the L 1 -constraint on the density function s. This leads to the saddle-point problem with η > 0 small enough and µ k ∈ R, the multiplier determined (by a line search algorithm) so that s k+1 L 1 (0,T ) = LT . In (5.50), p k is the solution of the adjoint problem (4.32) associated to the control function q k = λ k s k .
With the data and discretization used in Section 5.3, Table 3 provides some characteristics of the bangbang controls with respect to the parameter L ∈ [0, 1]. We use the numerical value ( , 1 ) = (10 −6 , 10 −4 ). The projection algorithm is initialized with λ 0 = 0 and s 0 = L. The variation of the bound λ with respect to L is small. Beyond L = 1/4, the optimal λ is almost constant. Actually, the amplitude of the control is driven by the data at the beginning of the time interval (corresponding to the colder period) where a large enough control is required, independently of L. Figure 10 depicts the optimal control corresponding to L = 1/4 while Figure 11 depicts the associated temperature θ(0, ·) at the road surface. Lower values of L lead to larger amplitudes λ and to an increasing number of switching points (from 0 and λ). We observe L 1 -norm values of the same order as in the previous section (e.g. q L 1 (0,T ) ≈ 4.57 × 10 8 for L = 1/10). Actually, for small values of L, we observe that the structure of the bang-bang control is similar to the optimal L 1 norm control (see Figure 6 for small α, e.g. α = 10 −7 ). The extreme value L = 1 (for which the control is active on the whole time period) leads to a value λ = 2.34 × 10 −2 : this value is in agreement with the result of the previous section which indicates that the optimal L 1 control satisfying the additional bound q L ∞ ≤ 200 does not allow to maintain the temperature greater than θ.  Table 3: Characteristics of the bang-bang controls with respect to L.

Command law
A priori, the optimal control we obtain in the previous section is non local in the sense that its value q(t) at the time t ∈ (0, T ) depends on the values of the data functions f 1 , f 2 , on the full time interval. This may not be realistic from a practical point of view. A first way to circumvent, in the spirit of [4] is to split the time interval into a finite number of sub-intervals and then compute optimal control on each subinterval.
A second way is to consider a source q which depends explicitly on the data. According to (4.28), for y 0 = 0, if the source q satisfy the condition q + f 1 (t) − f 2 (t)θ − σε(t)θ 4 ≥ 0, then the corresponding variable θ q satisfies θ q (0, t) − θ ≥ 0 for all t ∈ (0, T ). This suggests to consider, in our context where y 0 is close to zero, the following explicit source for some real δ ≥ 0 large enough, dependent of y 0 . Table 4 gives the L 1 -norm of q and the corresponding value of min((θ(0, ·) − θ) − ) for some values of δ. The value δ = 55 is enough to satisfies the condition θ(0, ·) ≥ 2 o C at the road surface. The corresponding L 1 -norm q L 1 (0,T ) ≈ 7.52 × 10 8 is of the same order as in the previous section.  Nevertheless, we observe that source of the form (5.51) is active on some period where the value of θ(0, ·) is (significantly) above θ. This is due to the irregularity of functions f 1 and f 2 . A third way is therefore to consider source term q which depends explicitly on the variable θ q , for instance as follows: for some reals θ m > θ, δ ∈ (0, T ) and a negative function f which depends only at time t on the temperature θ(s), s ∈ (0, t). Such source is active at time t ≥ δ if and only if two situations occurs: either when the value of θ(0, t − δ) is below θ, either when the value of θ(0, t − δ) is slightly above θ (precisely, in the range [θ, θ m ]) but decreases with respect to the time variable. Figure 12 depicts the source (5.52) associated to the value with θ m = 273.15 + 3 and δ = 1 hour and the corresponding temperature θ(0, ·). We take here simply for f a large enough constant so that (θ(0, ·) − θ) − L ∞ (0,T ) ≈ 0.08 leading to q L 1 (0,T ) ≈ 1.01 × 10 9 and q L ∞ (0,T ) ≈ 8.48 × 10 2 . Although larger than for the optimal controls computed earlier, these values, associated to the source q given by (5.52) explicitly in term of the history of the temperature, lead to promising results. 6 Perspectives 6.1 Control by a road surface heating source: a different mathematical approach If the source q is located at the road surface (i.e. y 0 = 0), we can construct a control by introducing a Signorini boundary condition as follows. Indeed, in the case where q is acting on the road surface, the right-hand side term of the first equation of (2.4) vanishes and the boundary condition at y = 0 becomes: we can eliminate q and replace (6.53) by the Signorini type conditions: which are equivalent to θ y (0, t) − f (t, θ(0, t)) ∈ γ(θ(0, t)) where γ : R → P(R) is the following multivalued and maximal monotone operator: The formal variational problem under conditions (6.54) is: find θ with θ(0, ·) ≥ θ such that for all v ∈ H 1 (0, h e ), he 0 c(y)θ t (y, t)v(y)dy + he 0 k(y)θ y (y, t)v y (y)dy + f (t, θ(0, t))v(0) + γ(θ(0, t))v(0) 0, which is an extension of (3.9) with a multivalued and monotone term γ(θ) added to θ 3 |θ|. Results of Section 3 may be extended to this variational inclusion as made in [9], leading to an existence and uniqueness result for (6.55), and then to a control q = −k(0)θ y (0, t)−f 1 (t)+f 2 (t)θ(0, t)+σε(t)θ 4 (0, t). From a mathematical viewpoint, it would be interesting to study the link between this control of Signorini type and the L 1 -optimal control exhibited in Section 4.

Optimal heating for the 2D diffusion-convection model
We are interested in the 2D model described in (2.1)-(2.3) and in Figure 3. In this context, the energy lost by the fluid is expressed as follows: For wider applications, it will be interesting to address this problem for a 3D model, allowing to take into account a longitudinal slant of the road.

Storage problem
Once the minimum energy has been calculated, there is the question of its origin. Indeed, in a context of energy transition, we must seek a renewable energy to heat the road. In an urban site, one can imagine using for example a heat network. In the case of an isolated road, one can seek to produce energy during the summer period, store it, and use it during winter. In this case, it is necessary to design a heat storage in a minimal manner. Let us denote by C the thermal capacity of the storage and denote by θ s its temperature. Considering the thermal losses of the storage and taking energy in this storage, we can express the time evolution of θ s as follows: C dθ s dt (t) = −A(t)λ s (θ s (t) − θ(y 0 , t)) − p(t, θ s (t)), (6.58) where λ s is the heat exchange coefficient between the storage (at temperature θ s ) and the road (at temperature θ(y 0 , ·)), A is a 0 − 1 function modeling the heat exchange activation and p denotes a law (assumed known) of thermal losses of the storage. In the one dimensional case, we can then consider the following system:                      c(y) ∂θ ∂t (y, t) − ∂ ∂y k(y) ∂θ ∂y (y, t) = (q(t) + A(t)λ s (θ s (t) − θ(y 0 , t))) δ(y 0 ), (y, t) ∈ (0, h e ) × (0, T ), − k(0) ∂θ ∂y (0, t) = f 1 (t) − f 2 (t)θ(0, t) − σε(t)θ 4 (0, t), ∂θ ∂y (h e , t) = 0, t ∈ (0, T ), θ(y, 0) = θ 0 (y), y ∈ (0, h e ), C dθ s dt (t) = −A(t)λ s (θ s (t) − θ(y 0 , t)) − p(t, θ s (t)), t ∈ (0, T ).

Conclusion
We have presented in this work an original approach to evaluate the minimal energy for heating a road in order to keep its surface frost free. With the help of a transient 1D thermal road model, we look for a heat source inserted in the pavement that allows to maintain the surface temperature above a value allowing the melting of snow or ice. Theoretical results are established for the uniqueness and existence of a solution to the direct problem and of optimal controls. For these results, we take into account the non-linear term due to the Stefan-Boltzmann law for which an emissivity can depend on time. We analyze the different controls obtained numerically according to we minimize the L 2 , L 1 or H 1 norm of the control. The L 1 norm has the most relevant physical meaning since it corresponds to the energy for heating the road and that we want to minimize. On the other hand, the control is sparse and leads to L ∞ norm higher than that corresponding to the H 1 or L 2 minimization. Also, for a better application, we study a bang-bang control taking only two values and allowing then to bound the control (the heat power). Some estimated optimal energies are given in Section 5: the total energy needed to keep the road surface temperature over 2 • C during a winter with snow is about 5.10 8 J 139 kWh per m 2 of road, with minimal and maximal values per m 2 respectively equal to 124 kWh and 213 kWh. Moreover the L ∞ norm of the optimal power q ranges in 240-500 W/m 2 . In [14], some experiments of heating roads by electric heating cables show needed power and energy equal respectively to 500 − 750 W/m 2 and 100 − 170 kWh/m 2 . In [10,7] needed power and energy are experimentally evaluated around 400 − 500 W/m 2 and 130 − 350 kWh/m 2 for deicing obtained by the circulation of a coolant in pipes inserted in the road. Even if the meteorology is not the same according to the above mentioned studies, one notices that the power and energy provided by our simulations have the same order of magnitude as those obtained experimentally for systems based on electric heating or coolant circulation in pipes. As mentioned in Section 6, we plan to extend our works to the heating system by circulation of a coolant in a bonding porous layer of the road. We will be able to compare our future numerical results with experimental measurements collected on the Egletons demonstrator described in Section 1.