A Limkilde

A Evgrafov

J Gravesen

A Mantzaflaris

Practical isogeometric shape optimization: Parameterization by means of regularization

Shape optimization based on Isogeometric Analysis (IGA) has gained popularity in recent years. Performing shape optimization directly over parameters defining the CAD geometry, such as for example the control points of a spline parametrization, opens up the prospect of seamless integration of a shape optimization step into the CAD workflow.

One of the challenges when using IGA for shape optimization is that of maintaining a valid geometry parametrization of the interior of the domain during an optimization process, as the shape of the boundary is altered by an optimization algorithm. Existing methods impose constraints on the Jacobian of the parametrization, to guarantee that the parametrization remains valid. The number of such validity constraints quickly becomes intractably large, especially when 3D shape optimization problems are considered.

An alternative, and arguably simpler, approach is to formulate the isogeometric shape optimization problem in terms of both the boundary and the interior control points. In order to ensure a geometric parametrization of sufficient quality a regularization term, such as the Winslow functional, is added to the objective function of the shape optimization problem.

We illustrate the performance of these methods on the optimal design problem of electromagnetic reflectors and compare their performance. Both methods are implemented for multipatch geometries, using the IGA library G+Smo and the optimization library Ipopt. We find that the second approach performs comparably to a state of the art method with respect to both the quality of the found solutions and computational time, while its performance in our experience is more robust for coarse discretizations.

Introduction

Isogeometric analysis (IGA) introduced in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] is a Galerkin method that uses splines to approximate both the geometric domain and the solutions to partial differential equations (PDEs). Splines are commonly used in computer aided design (CAD) and IGA is an attempt to bridge the gap between simulation and design [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF]. This makes it beneficial for shape optimization as the optimization can be performed directly over parameters defining the CAD geometry, for example the control points of a spline parametrization, and it opens up the prospect of seamless integration of a shape optimization step into the CAD workflow.

One of the key challenges when using IGA in general, is that one needs a parametrization of the interior of the physical domain, on which the PDE is posed [START_REF] Shamanskiy | Isogeometric parametrization inspired by large elastic deformation[END_REF][START_REF] Hinz | Elliptic grid generation techniques in the framework of isogeometric analysis applications[END_REF][START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF][START_REF] Martin | Volumetric parameterization and trivariate b-spline fitting using harmonic functions[END_REF]. This parametrization is used to pull back the weak form of the PDE to the parameter domain where the basis splines (B-splines) are defined. The choice of parametrization can affect the accuracy of the resulting IGA discretization [START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF][START_REF] Xu | Optimal analysis-aware parameterization of computational domain in isogeometric analysis[END_REF] and, at the very least, the parametrization should be valid (a bijective map), that is, its Jacobian determinant should be non-zero. One approach to constructing a valid parametrization in 2D is to search for the one which has harmonic inverse. In [START_REF] Hinz | Elliptic grid generation techniques in the framework of isogeometric analysis applications[END_REF] this property is reformulated as a nonlinear PDE and the parametrization is found by solving this PDE. In [START_REF] Hinz | The role of pde-based parameterization techniques in gradient-based iga shape optimization applications[END_REF] this PDE based parametrization technique is used for a gradient based shape optimization algorithm with IGA. In [START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF][START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF] the same property is attained by minimizing the Winslow functional [START_REF] Winslow | Numerical solution of the quasilinear poisson equation in a nonuniform triangle mesh[END_REF]. The method can be made more flexible by the use of adaptive splines [START_REF] Falini | Planar domain parameterization with thb-splines[END_REF], that allow to enrich the feasible region near to complex boundaries. In the recent works [START_REF] Shamanskiy | Isogeometric parametrization inspired by large elastic deformation[END_REF][START_REF] Sajavičius | Template Mapping Using Adaptive Splines and Optimization of the Parameterization[END_REF] the approach of parametrizing a complex domain by deforming a given template is explored. In [START_REF] Pan | Low-rank parameterization of volumetric domains for isogeometric analysis[END_REF] the focus is on producing parametrizations with low-rank with respect to the coefficient tensor.

When using IGA for shape optimization the challenge of finding a valid parametrization is even more important, since the shape of the physical domain changes during an optimization process. This means that a valid parametrization needs to be maintained during this process.

To guarantee that the parametrization remains valid during the optimization process, shape optimizations methods based on IGA often rely on constraints on the Jacobian of the parametrization [START_REF] Gangl | Isogeometric simulation and shape optimization with applications to electrical machines[END_REF][START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF]. These can be enforced either by using injectivity cones or by using the spline coefficients of the Jacobian determinant [START_REF] Xu | Parameterization of computational domain in isogeometric analysis: methods and comparison[END_REF]. However the number of constraints needed quickly becomes very large, especially in 3D. Furthermore, when using the coefficients of the Jacobian determinant for the constraints, as we will do in this work, it may be necessary to expand the Jacobian determinant on a finer spline space, which increases the number of constraints even further.

In this work we will compare an existing approach to IGA shape optimization, relying upon reparametrizing the domain, with a simple approach to maintaining a valid parametrization without the use of explicit validity constraints. Namely, one lets the positions of all the control points that define the parametrization and the shape of the domain enter the formulation as independent optimization variables and adds a regularization term to drive the optimization towards a design with a valid parametrization. Such an approach has been considered in the context of shape optimization in mechanics in [START_REF] Scherer | A fictitious energy approach for shape optimization[END_REF], and has to the best of our knowledge, only been considered very briefly in the context of shape optimization with IGA in [START_REF] Friederich | Isogeometric structural shape optimization using a fictitious energy regularization[END_REF]. It remains a question whether this approach performs comparably to state of the art methods, and the aim of the work is to investigate exactly this question.

Different frameworks for shape optimization with IGA exist in the literature. For example in [START_REF] López | Cad-compatible structural shape optimization with a movable bézier tetrahedral mesh[END_REF] a tetrahedral mesh is used to represent the interior of the computational domain, while still representing the boundary using splines. The authors rely on external mesh generator software, while we in this work aim to avoid mesh generation by instead maintaining a parametrization of the interior. In [START_REF] Ghasemi | A level-set based iga formulation for topology optimization of flexoelectric materials[END_REF] the shape is represented implicitly as the level set of a function. The PDE is then posed on a design domain, which remains constant during the optimization, and the level set function enters the formulation by the ersatz material approach. This framework has the advantage that it allows for changes in the topology of the shape, however at the same time the final shape is represented as a level set and therefore post processing is required to represent this shape using splines, which is necessary for importing the shape into CAD software. In our work the optimization is performed directly on the spline representation and thus the result can be readily imported into CAD software after the optimization. Another possible approach is to use the isogeometric representation of geometry combined with the boundary element method (IGABEM) [START_REF] Liu | Shape optimization of sound barrier using an isogeometric fast multipole boundary element method in two dimensions[END_REF][START_REF] Chen | Structural shape optimization of three dimensional acoustic problems with isogeometric boundary element methods[END_REF]. In IGABEM the PDE is reduced to an integral formulation on the boundary of the domain, and therefore maintaining a parametrization of the interior of the domain is avoided altogether. However in this case one has to deal with the standard complications of BEM, in particular the fact that the system matrices are dense, non-symmetric, and costly to compute. Additionally the Greens function for the considered PDE has to be known, which is not always the case.

In this work we illustrate that the simple regularization based approach is able to handle complicated geometries, by comparing its performance to a shape optimization approach based on using a linearization of the Winslow minimization problem as a parametrization strategy, and employing locally refined splines to represent the Jacobian determinant. The method closely resembles the one in [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF]. The main difference is that we use Truncated Hierarchical Basis splines (THB-splines) [START_REF] Giannelli | Thb-splines: The truncated basis for hierarchical splines[END_REF], which possess the partition of unity property, to refine the spline space in which the determinant is expanded locally. This reduces the number of constraints needed compared with tensor product global refinement.

We will apply the two methods and compare their performance on the shape optimization problem of designing electromagnetic reflectors. In this problem we have two metallic reflectors in a dielectric medium and search for a shape that maximizes the electrical energy close to a chosen point. The same problem has been studied with topology optimization in [START_REF] Aage | Topology optimization of metallic devices for microwave applications[END_REF][START_REF] Wadbro | Topology and shape optimization of plasmonic nano-antennas[END_REF][START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF], and with IGA in [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF].

The methods are implemented for multipatch geometries, using the IGA library G+Smo1 and the optimization library Ipopt2 . A benchmark study on the performance of different optimization algorithms in the context of structural optimization can be found in [START_REF] Rojas-Labanda | Benchmarking optimization solvers for structural topology optimization[END_REF]. For a review of implementation aspects of IGA in general we refer the reader to [START_REF] Nguyen | Isogeometric analysis: an overview and computer implementation aspects[END_REF] and for details of the G+Smo library we refer the reader to [START_REF] Jüttler | Geometry+ simulation modules: Implementing isogeometric analysis[END_REF][START_REF] Mantzaflaris | An overview of geometry plus simulation modules[END_REF]. The code used for this work can be found at https://github.com/gismo/shapeopt.

The paper is organized as follows. In section 2 we outline the relevant notation, and in sections 3 and 4 we describe the two methods we are going to compare. In section 5 we apply these methods to the aforementioned shape optimization problem, and discuss and compare the performance of the two approaches. We end the paper with some discussion and conclusions. Some of the more technical details are presented in the appendices.

Preliminaries and notation

Let us consider the following PDE-constrained shape optimization problem: max

Ω∈O ad E(Ω, u), (1a) s
.t. a Ω (u, v) = Ω (v) for all v ∈ V, (1b)
where O ad is a set of admissible shapes, E is the objective and (1b) is the governing PDE in the weak form. Within the IGA framework both Ω and u will be approximated numerically using splines. Namely, we have Ω = G(]0, 1[ d ), with the parametrization

G(ξ) = N g i=1 c i R g i (ξ), (2)
where c i ∈ R d are the control points, N g is the number of control points, and R g i are the basis (B) splines. In this work, unless specifically stated otherwise, we will use tensor product B-splines. The superscript g indicates that the B-splines R g i (that is, their degrees and knotvectors) are specific to the geometry representation. Within the shape optimization framework it will sometimes be necessary to distinguish between boundary and inner control points. We will therefore introduce the notation

c = c b c i , (3)
where c b are the boundary control points and c i are the inner control points. The Jacobian J = ∂G ∂c will also play an important rôle in the forthcoming development.

Similarly to (2), we approximate the state of our system as a pulled back spline

u h = N i=1 u i R i • G -1 , (4)
where R g i , i = 1, . . . , N are B-splines. The expansion coefficients u i , i = 1, . . . , N will be found by solving a system of linear algebraic equations

K c u = f c .
As standard in the Galerkin approach, the elements of the stiffness matrix K c and the load vector

f c are computed as K c,i,j = a G(]0,1[ d ) (R j • G -1 , R i • G -1 ), and f c,i = G(]0,1[ d ) (R i • G -1
). Note that the dependence of K c , f c on the control points is encapsulated in (2) and ( 4).

Already at this point the importance of the geometry parametrization should be apparent. Indeed, at the very least it should be an invertible map, which is used to pull back the weak form of the PDE defined on the physical domain Ω into the parameter domain ]0, 1[ d . In particular, for all ξ ∈]0, 1[ d it is necessary that det(J(ξ)) > 0.3 A sufficient condition, which guarantees the validity of the parametrization, is discussed in Appendix A.

3 Boundary-driven approach to IGA shape optimization

In this section we will consider one possible approach to IGA shape optimization, which follows the ideas developed in [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF][START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF]. Within this framework we formulate the optimization problem in terms of boundary control points c b . The collection of interior control points c i for the geometry parametrization is treated as an implicit function of c b , see Appendix B. Additionally, the parametrization validity constraints det(J) > 0, or a sufficient condition for these (cf. Appendix A), have to be explicitly included into the problem formulation.

In order to compute domain parametrizations of high quality we rely upon minimizing the Winslow functional (Section B.2). However, to avoid solving a non-linear optimization problem at each shape optimization iteration, we construct a quadratic approximation to the Winslow functional around a reference parametrization, and update the reference parametrization when it becomes necessary to do so. Specifically, given a reference parametrization G 0 defined by the control points c 0 , to find a new parametrization we consider the quadratic programming problem

min ∆c i 1 2 ∆c T H(c 0 )∆c + ∇W (c 0 ) T ∆c + W (c 0 ), ( 5
)
where W is the Winslow functional, and H is its Hessian. The minimizer of this problem can be found by solving a linear system H(c 0 )∆c = -∇W (c 0 ).

Using (3), this can be restated as

H c i ,c i ∆c i = -∇W c i -H c i ,c b ∆c b , (6)
where

[H c i ,c i ] ij = ∂ 2 W /∂c i i ∂c i j , [H c i ,c b ] ij = ∂ 2 W /∂c i i ∂c b j and [∇W c i ] i = ∂W /∂c i i .
The new parametrization is then defined by the control points given by c = c 0 + ∆c.

With this in mind, to approximate (1) numerically we solve a sequence of subproblems max

∆c b E(c, u), (7a) s
.t. K c u = f c , (7b)
H c i ,c i ∆c i = -∇W c i -H c i ,c b ∆c b , (7c)
c = c 0 + ∆c, ( 7d
) d ≥ ε, (7e)
c b L ≤ c b ≤ c b U , (7f)
where c 0 is the reference parametrization, ∆c = (∆c b , ∆c i ), and (7e) is the sufficient condition for the validity of the parametrization discussed in Appendix A. Each time we solve the subproblem [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF] we update the reference parametrization. In the model problem considered in this work we saw no further progress after 5 to 10 reparametrizations.

Reference parametrizations can be computed as follows. We minimize the Winslow functional as described in Appendix B, and check if the sufficient condition d > 0 is violated. If it is, then this condition is too strict and should be relaxed. To facilitate this we refine the spline space S det where we compute expansion coefficients d of det J. To reduce the number of constraints resulting from such refinement steps, we utilize local refinement. Specifically, we use Truncated Hierarchical B-splines (THB-splines) as basis functions. Note that it is important here to use the truncated version of hierarchical splines, since the partition of unity property (cf. [START_REF] Giannelli | Thb-splines: The truncated basis for hierarchical splines[END_REF]) implies that the spline control polygon converges locally to function values. Note that other locally refinable splines with this property are available, for example Polynomial splines over Hierarchical T-meshes (PHTsplines) [START_REF] Deng | Polynomial splines over hierarchical t-meshes[END_REF] or Locally Refinable (LR B-splines) [START_REF] Johannessen | Isogeometric analysis using lr b-splines[END_REF]. For a comparison of these methodologies we refer the reader to [START_REF] Li | A survey on the local refinable splines[END_REF].

The refinement strategy we employ is as follows. For all indices i that have negative spline expansion coefficient d i ≤ 0 of det J, we refine the support of the associated basis function R det i . This is repeated until d i > 0 for all i = 1, . . . , N det . 4 In the subproblem [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF] we then put ε = ρ•min i d i with ρ = 0.25, see [START_REF] Falini | Planar domain parameterization with thb-splines[END_REF].

The full optimization loop is illustrated in Figure 1. Note that we have chosen to carry out the spline space refinement described above only when the reference parametrization, and therefore also the subproblem [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF], is updated. This allows us to keep the number of constraints constant when solving [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF] numerically and therefore employ off-the-shelf optimization software. The initial guess for the non-linear Winslow optimization problem is generated using Coons' patches, see Appendix B.

4 Regularization-driven approach to IGA shape optimization

In this section we discuss an alternative approach to shape optimization using IGA, which does not involve explicit constraints on det J. The positions of the inner control points c i enter this formulation as independent optimization variables, in the same way as c b . Consequently, we do not need to explicitly compute a domain parametrization, as this will be part of the outcome of the optimization process.

To this end we add the Winslow functional W as a regularization term to the objective function. Its role is to penalize configurations of control points that result in poor parametrizations. This idea has been used previously in the context of shape optimization in mechanics [START_REF] Friederich | Isogeometric structural shape optimization using a fictitious energy regularization[END_REF][START_REF] Scherer | A fictitious energy approach for shape optimization[END_REF]. Thus for a regularization parameter τ > 0 we consider the optimization problem

min c τ W (c) -E(c, u), (8a) s
.t. K c u = f c , (8b)
c L ≤ c ≤ c U . (8c)
We put W (c) = ∞ if det J ≤ 0 at one of the quadrature points used for the integration when calculating W (c). In this way, when the optimization algorithm enforces e.g. the standard sufficient decrease condition, such a choice ensures that the chosen step will always have positive determinant at the quadrature points. This does not guarantee that it is positive everywhere, but it means that the numerics will not collapse due to a division by zero.

The regularization parameter τ needs to be tuned for the specific problem at hand. If it is too large the minimization will find a design with a small value of the Winslow functional but disregarding the objective E(c, u). If it is too small the optimization will find positions of the control points that have a low objective E(c, u), but with a poor parametrization, which might give a large discretization error of the discretized PDE. The appropriate values of τ would lead to a compromise between these two extreme situations. One simple strategy for choosing such a value is to solve a sequence of problems [START_REF] Compton | Bow-tie antennas on a dielectric half-space: Theory and experiment[END_REF] for decreasing values of τ . 5The two points above constitute the main drawback of this method. Namely that, in contrast to a constraint based method, we cannot guarantee that the parametrization is valid in between the quadrature points. Additionally one has to find a suitable value of the regularization parameter τ , which is problem-dependent.

Case study: Optimization of electromagnetic reflectors

In this section we will consider a 2D shape optimization problem originating from the field of electromagnetism. Our goal is to design a reflector that concentrates electrical energy in a desired area. This problem will serve as a model problem for comparing the two optimization approaches outlined in Sections 3 and 4.

Physical model

We consider a two dimensional scattering problem where a plane wave with frequency f travels in a dielectric (air) and is scattered by two symmetric metallic (gold) reflectors, as depicted in Figure 2. Let cr and µ r denote the complex permittivity and permeability of the medium. Using the first order absorbing boundary condition [START_REF] Jin | The finite element method in electromagnetics[END_REF] at the boundary Γ t of the truncated domain, the electromagnetic field û should satisfy the following PDE:

∇ • 1 cr ∇û + k 2 0 µ r û = 0 in Ω, (9a)
∂(û -u i ) ∂n + (jk 0 + 1 2r t )(û -u i ) = 0 on Γ t . (9b)
In the equations above, k 0 = 2π √ 0 µ 0 is the wave number and 0 , µ 0 refer to the permittivity and permeability, respectively, of free space. The imaginary unit is denoted by j, the radius of the truncated domain is given by r t , and u i is the incident plane wave, given by

u i (x, y) = e -jk0 √ cr µrx .
The objective function of the shape optimization will be given by

E(c, u) = Ω δ|û| 2 dx,
where δ is a Gaussian bell-function with α = 0.1. Thus we aim to focus the incoming energy in the vicinity of the origin (0, 0). The physical parameters that we use are given in Table 1. The complex permitivity of the reflector is calculated as s cr = r,gold -j σ ω 0 . The weak statement of the PDE ( 9) is to find û ∈ H 1 (Ω) such that for all test functions v ∈ H 1 (Ω) the following equality holds

δ(x, y) = e -x 2
Ω 1 cr ∇û • ∇vdx + k 2 0 Ω µ r ûvdx + jk 0 + 1 2r t Γt 1 cr ûvds = 1 cr Γt ∂u i ∂n + jk 0 + 1 2r t u i vds. ( 10
)
Due to the symmetry we only consider the upper half of the geometry shown in Figure 2. To accommodate the change of material parameters between the metallic reflector and the surrounding dielectric medium, we will split the domain into five patches, one for the reflector and the other four for the surrounding air. The layout is shown in Figure 3. 6 Each patch is parametrized using splines of degree p = 2 as described in Section 2. We will use strong patch coupling to enforce C 0 continuity at the patch interfaces, while noting that other alternatives exists, for example weak coupling [START_REF] Hu | Skew-symmetric nitsche's formulation in isogeometric analysis: Dirichlet and symmetry conditions, patch coupling and frictionless contact[END_REF] and strong C 1 coupling [START_REF] Chan | Strong multipatch c1-coupling for isogeometric analysis on 2d and 3d domains[END_REF].

The shape of the reflector is represented via the patch interfaces, which consist of 4 spline curves with C 0 continuity at the 4 corners. Such a description, for example, allows shapes that are only piecewise smooth such as e.g. the classical bowtie antenna [START_REF] Compton | Bow-tie antennas on a dielectric half-space: Theory and experiment[END_REF]. There are other applications of IGA where even more complex geometries are relevant, for example when modelling crack propagation as in [START_REF] Peng | Linear elastic fracture simulation directly from cad: 2d nurbs-based implementation and role of tip enrichment[END_REF][START_REF] Peng | Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth[END_REF].

Using the parametrizations we can pull back the equation [START_REF] Deng | Polynomial splines over hierarchical t-meshes[END_REF] to the parameter domain and apply the Galerkin method to it, which ultimately results in a system of linear algebraic equations; see Appendix C for details. To evaluate the integrals involved we use element-wise Gauss-Legendre quadrature, owing to their immediate availability in the G+Smo library, while noting that more efficient alternatives exists [START_REF] Calabro | Fast formation of isogeometric galerkin matrices by weighted quadrature[END_REF][START_REF] Bartoň | Gauss-galerkin quadrature rules for quadratic and cubic spline spaces and their application to isogeometric analysis[END_REF].

Results with boundary-driven approach

In this section we apply the method described in Section 3 to our model problem. We start with an initial design where the reflector has the shape of a circle. We consider two different spline spaces in which to approximate the PDE [START_REF] Cottrell | Isogeometric analysis: toward integration of CAD and FEA[END_REF], namely where the knotvectors used for representing the geometry are refined uniformly either 3 or 4 times. 7 Both spline spaces have degree p = 2, and the number of degrees of freedoms are N coarse = 2548 and N fine = 9300, respectively. We will refer to these as the coarse and fine meshes.

We will use a tolerance tol = 10 -3 when solving the subproblems [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF] and a fixed number of reparametrizations, namely 10 when using the coarse mesh and 5 when using the fine mesh. We observed that using more reparametrizations did not lead to significant improvements in the design. In our experiments the results with this method are sensitive to the number of quadrature points used when calculating the Winslow functional. To produce the results presented here we use 12 quadrature points per knot interval, to avoid under-integration 8 .

In our implementation we use the interior point solver Ipopt to solve the subproblems [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF]. As we solve a sequence of subproblems a warm start is available, namely the minimizer from the previous subproblem; so an interior point algorithm might not be the best choice of optimization algorithm [START_REF] John | Implementation of warm-start strategies in interior-point methods for linear programming in fixed dimension[END_REF]. However with the following parameter tuning we found the method to work well with Ipopt. Namely, after solving one of the subproblems (7) some of the design bounds (7f) will be active; however, since Ipopt is an interior point algorithm, the starting point for the subsequent subproblem will be pushed away from the boundary as controlled by the parameter bound push. We found that this parameter needs to be lower than the default value since the constraints on det J are quite sensitive and a relatively small perturbation of the control points might violate these constraints, which is undesirable. The default value is 0.1, but in the experiments we set it to 10 -5 instead.

Another key parameter is the barrier parameter mu init. Specifically, we use the monotone strategy, where the barrier parameter is monotonically decreased as the optimization algorithm progresses. However, if this parameter is too large in the beginning of the algorithm we found that it will push the design towards configurations with large det J. To remedy this, instead of the default value 0.1 we use 10 -4 . For more information about the optimization algorithm implemented in Ipopt and its parameters see [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF].

In Figure 4 the design at different stages of the optimization is presented, when using the coarse mesh.

We observe that the design becomes increasingly hard to parametrize as the objective increases each time we change the reference parametrization. In Figure 5 the objective is plotted against the number of iterations. We see that the objective function increases after the reference parametrization is changed, but relatively quickly reaches a plateau. We already use a fairly large tolerance of 10 -3 for the stopping criterion when solving the subproblems [START_REF] Christiansen | A non-linear material interpolation for design of metallic nano-particles using topology optimization[END_REF], however this behaviour indicates that it might help to relax the stopping criterion even further in these subproblems to improve the overall efficiency of the method. However, to allow for a fair comparison between the two methods, we do not investigate this further and use the same tolerance for both methods. The final objective, after 10 reparametrizations when using the coarse mesh, is E h = 1.803. However if we calculate the objective with a mesh that is refined uniformly twice we get E h/4 = 1.556, that is, a 16% difference. This means that there is actually a significant discretization error at this refinement level, which leads to an artificially large objective value.

As described in Section 3 we use local refinement to adapt the constraints on det J to the current design. This is done by using local refinements in the areas where the spline coefficients are non positive. The resulting meshes are plotted in Figure 6 to illustrate where the refinement is needed. We see that it is primarily inside the reflector and near the reflector-air interface, that this Figure 4: The designs at different stages of the optimization process, when using validity constraints and the coarse mesh. The reflector is outlined with a black line, and the control points of this boundary is colored black. The grey lines are parameter lines mapped with the geometry map, to illustrate the parametrization. Figure 6: The mesh (knot lines), used for representing det J during the subproblems (7) for the constraints [START_REF] Falini | Planar domain parameterization with thb-splines[END_REF]. The number of constraints is given by N det and the designs are obtained using the coarse mesh.

Figure 7: The final design using validity constraints and the fine mesh, E h/2 = 1.638 refinement is applied. We see that the number of constraints increases during the optimization process, but by no more than a factor of two. If we were to use uniform refinement the number of constraints would increase by more than a factor of 9.

The final design, after 5 reparametrizations, when using the fine mesh, is shown in Figure 7. We see that it is similar to the design obtained using the coarse mesh. The final objective here is E h = 1.638. After the mesh is refined uniformly the electrical energy is E h/4 = 1.628, which is only a 0.6% difference. The objective function during the optimization process is plotted in Figure 5. We see that the algorithm converges faster when using the fine mesh as no progress was observed after 5 reparametrizations.

Results with regularization-driven approach

In this section we will present the results obtained with the method described in Section 4. With this method we perform the optimization with all control points as optimization variables while using the Winslow functional as a regularization term. Again we use the interior point solver Ipopt with a tolerance of tol = 10 -3 for solving the problem [START_REF] Compton | Bow-tie antennas on a dielectric half-space: Theory and experiment[END_REF]. The regularization parameter is set to τ = 1 8 . One can compare the design evolution shown in Figure 8 with those obtained previously, see Figure 4. The designs obtained using the regularization approach seem to have more regular parametrizations compared to those in Figure 4.

The final objective is E h = 1.684. If we calculate the electrical energy for this design on a twice refined mesh we get E h/4 = 1.546, that is, a difference of 9% . This is less than the 16% we observed when using the boundary-driven method. This increase in accuracy might be due to the parametrization being of higher quality.

In Figure 9 we plot the electric energy E h , the regularization term τ W and the objective function τ W -E h .

When using the method with the fine mesh we get the final design shown in Figure 10. We see that the shape of the reflector is very similar to the one shown in Figure 8d. The main difference is that when using the fine mesh the parametrization is more regular, since the error in the discretization of the PDE is smaller, and therefore the optimization cannot exploit it to the same extent. This is especially notable at the bottom of the reflector were the inner control points where moved away from the point of interest when using the coarse mesh, as seen in Figure 8d. The final objective is E h = 1.545 and when evaluating it on a refined mesh we get the same Figure 8: The designs at different stages of the optimization process when using the regularization approach. The reflector is outlined with a black line, and the control points that controls this boundary is colored black. The grey lines are parameter lines mapped with the geometry map, to illustrate the parametrization. The objective function during the optimization process when using the regularization based approach, with the fine and coarse mesh.

E h , coarse τ W , coarse τ W -E h , coarse E h/4 , coarse E h , fine τ W , fine τ W -E h , fine E h/4 , fine
Figure 10: The final design when using the regularization approach and the fine mesh, E h/2 = 1.545.

result E h/4 = 1.545 with the difference at the 5th digit. The objective during the optimization process is plotted in Figure 9. The behaviour is similar for the two meshes, however the tolerance tol = 10 -3 is reached with fewer iterations when using the fine mesh.

Comparison and discussion

In Table 2 we summarize the performance of the two methods. We report the objective after a fixed number of iterations, in this case after 100 iterations, the final objective computed on 3 different refinement levels and the average running time per iteration. We observe that the average execution time per iteration is the same order of magnitude for the two methods. The main portion of the running times is spent on solving the state equation and computing the gradient of the objective function. The difference in the running time that we observed between the two methods might be due to a different number of function evaluations per iteration needed for trial steps of the algorithm. If we were to consider a larger problem, for example in 3D, the large number of validity constraints would likely lead to an increase in running time for the boundary-driven approach. For the boundary-driven approach the computation of reference parametrizations by minimizing of the Winslow functional only accounts for 2.4% and 0.25% of the total running time for the coarse and the fine grid respectively.

Regarding the quality of the designs we find, we note that shape optimization problems are prone to having many local optima, so it could be that the two methods find two different local optima. Therefore it can be futile to directly compare objective values. That being said, we observe that the boundary-driven approach happens to find solutions with slightly higher (better) objective value. On the other hand, the regularization based approach seems to estimate the objective value more accurately on coarser meshes, since our results were more reliable using this method, probably due to the better quality of the parametrizations it produced.

Also, observe that with the regularization based approach we only need the objective E h , the Winslow functional and their first order derivatives. In addition we can solve a single optimization problem with design bounds as the only constraints. This means that the method is significantly easier to implement.

Method

E iter=100 E h E h/2 E h/4 E h -E h/4 E h/4
Avg time per iteration Linearizations coarse mesh

Conclusion

We described and compared two methods for shape optimization on spline-based representations.

One uses validity constraints to enforce the validity of the geometry parametrization. The other uses a regularization term, and thus avoids both the validity constraints and the need for an explicit parametrization strategy altogether. We demonstrated how this simple approach performed comparably to the more complicated approach in terms of the final design, while requiring similar running times for the 2D problem we considered. The regularization based approach seems to produce more reliable results and it is in addition much simpler to implement, since we only need the objective, the Winslow functional and their first order derivatives. These results are encouraging and we plan to use the regularization based approach for 3D problems, where we expect that the efficiency advantages of the regularization approach will be more prominent, since the number of validity constraints for 3D parametrizations grows quickly. For obtaining competent overall running times, one can couple the approach with low-rank approximation techniques, which have been recently introduced in IGA [START_REF] Mantzaflaris | Low rank tensor methods in Galerkin-based isogeometric analysis[END_REF][START_REF] Scholz | Partial tensor decomposition for decoupling isogeometric Galerkin discretizations[END_REF] and have proved efficient in the frame of PDE-constrained optimization [START_REF] Bünger | A low-rank tensor method for pde-constrained optimization with isogeometric analysis[END_REF]. Now we can use this expansion to derive a sufficient condition to replace the pointwise parametrization validity constraints det(J) > 0, by requiring that

d ≥ ε, (11)
where ε is a small positive algorithmic parameter. This condition guarantees that det J > 0 since B-splines are non-negative and form a partition of unity. This is not a necessary condition, meaning that we might have d i ≤ 0 for some i even though det J(ξ) > 0, ∀ξ ∈]0, 1[ d . But if the spline space S det is refined then the spline expansion coefficients will move closer to values of the spline. So if ε is small enough and det J > 0 then the constraint is likely to be satisfied for a sufficiently refined spline space S det .

B Domain parameterization techniques

In this section we review some techniques for finding a parametrization of the interior given the boundary. In IGA this comes down to finding the position of the inner control points given boundary control points. We do not aim to give a thorough review of all the techniques that are available, as this is out of the scope of this work. We will only introduce the methods that are related, or directly used, by the two shape optimization methods considered in this work

B.1 Coons patch and spring method

Two simple methods for constructing a grid of control points are the Coons patch [START_REF] Farin | Discrete coons patches[END_REF] and the spring method [START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF]. They both produce inner control points that depend linearly on the boundary control points, however they only produce valid parametrizations for geometries that are not too complicated. In this work we will use these methods for finding an initial guess for the optimization based approach that is described below.

B.2 Optimization-based techniques

A more complex, and in general nonlinear, class of parametrization methods consists of optimizationbased methods. Here the geometry map G is chosen such that it minimizes a quality metric w(ξ)

min G ]0,1[ d w(ξ) dξ, (12)
s.t. G| ∂[0,1] d = γ, ( 13
)
where γ is a given boundary curve.

There are several different quality metrics to choose from, cf. [START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF][START_REF] Xu | Constructing analysis-suitable parameterization of computational domain from cad boundary by variational harmonic method[END_REF]. In this work we will consider the Winslow functional which is given in terms of the Jacobian matrix J as

W = ]0,1[ d w(ξ) dξ, (14)
with

w = tr J T J det J . ( 15
)
In 2D the Winslow functional has the nice property that its minimizer has a harmonic inverse [START_REF] Gravesen | Planar parametrization in isogeometric analysis[END_REF]. This guarantees that the minimizer of ( 12) is unique and bijective, i.e. det J = 0. It should be noted that this minimizer is not necessarily a spline, so looking for a spline parametrization on the form (2) by minimizing [START_REF] Farin | Discrete coons patches[END_REF] with the Winslow functional will only guarantee a valid parametrization if the spline space used for the parametrization has high enough resolution.

Within the shape optimization context we also need first and second order partial derivatives of w. The derivative with respect to a parameter α is given by ∂w ∂α = 2(det J) -1 tr J T ∂J ∂α -tr J -1 ∂J ∂α tr(J T J) det J ,

where we used the relation ∂ ∂α det J = det J tr J -1 ∂J ∂α . The second order derivative is given by ∂ 2 w ∂α∂β = 2(det J) -1 tr ∂J ∂α

T ∂J ∂β -2(det J) -1 tr J -1 ∂J ∂α tr J T ∂J ∂β -2(det J) -1 tr J -1 ∂J ∂β tr J T ∂J ∂α + tr(J T J) det J tr J -1 ∂J ∂α tr J -1 ∂J ∂β

+ tr(J T J) det J tr ∂J ∂β T J -1 ∂J ∂α J -1 , (17)
using the fact the fact that ∂ ∂α J -1 = -J -1 ∂J ∂α J -1 and assuming that J depends linearly on α, which is the case when α is a coordinate of a control point. Calculation of ∂w ∂α and ∂ 2 w ∂α∂β can be implemented as an assembly of a linear and bilinear forms within IGA framework. In G+Smo this can, for instance, be accomplished using the gsExprEvaluator class that is typically employed for isogeometric stiffness matrix assembly purposes.

C IGA discretization details

Using the spline parametrizations of patch geometries we pull back the weak form [START_REF] Deng | Polynomial splines over hierarchical t-meshes[END_REF] to the parameter domain, which results in the following equation:

1 cr ∂u i ∂n • G + jk 0 + 1 2r t u i • G v ∂G ∂t dt (18)
where we have u = û • G, v = v • G and where t is the parameter on the boundary. After applying the Galerkin method to [START_REF] Hansen | Regularization tools: a matlab package for analysis and solution of discrete ill-posed problems[END_REF], we arrive at the linear system of linear algebraic equations

Au = (K + M + T )u = f (19)

Figure 1 :

 1 Figure 1: Flowchart of the optimization algorithm.

ReflectorFigure 2 :

 2 Figure2: Sketch of the shape optimization problem. The goal is to find a shape of the reflector that maximizes the field close to a point.

4 •Table 1 Figure 3 :

 413 Figure 3: Patch layout.

  (a) The initial design, E h = 0.200 (b) After 1 reparametrization, 67 iterations, E h = 1.230 (c) After 2 reparametrizations, 143 iterations, E h = 1.484 (d) After 10 reparametrizations, 579 iterations, E h = 1.803

Figure 5 :

 5 Figure 5: The objective function during optimization process for the coarse and fine meshes. The vertical lines indicate where the parametrization is updated. E h/4 is the objective calculated on a refined mesh.

( a )

 a Initial mesh, N det = 380 (b) After 1 reparametrization, N det = 392. (c) After 2 reparametrizations, N det = 448. (d) After 10 reparametrizations, N det = 572.

  (a) The initial design, E h = 0.202 (b) 67 iterations, E h = 1.461 (c) 143 iterations, E h = 1.593 (d) The final design after 432 iterations , E h = 1.683

Figure 9 :

 9 Figure9: The objective function during the optimization process when using the regularization based approach, with the fine and coarse mesh.

2 0

 2 ∇u • J -T ∇v |det J| dξk

Table 2 :

 2 Comparison of the two methods from Section 3 and 4 on two different meshes. E iter=100 is the objective after 100 iterations. E h is the objective computed on the mesh used in the optimization. E h/2 and E h/4 are the objectives computed after refining the mesh uniformly, once and twice. The execution time was measured on a 64 bit HP EliteBook 840 G4 with and Intel(R) Core(TM) i7-7500U CPU, with clock rate of 2.70 GHz.

		1.482	1.803 1.589 1.556	16%	9.68 sec
	Regularization coarse mesh	1.529	1.684 1.561 1.546	9%	11.54 sec
	Linearizations fine mesh	1.495	1.638 1.629 1.628	0.6%	53.5 sec
	Regularization fine mesh	1.482	1.545 1.545 1.545 0.003%	71.0 sec

https://github.com/gismo

https://github.com/coin-or/Ipopt

It is equivalent to require that det(J(ξ)) < 0, however in this work we will use the constraint det(J(ξ)) > 0

In practice we terminate this procedure either if d > 0 or when a maximum level of refinement (7 in our numerical experiments) is attained. The latter termination criterion has not been observed in our experiments.

We should note that the literature on regularization is quite extensive, see for example[START_REF] Hansen | Regularization tools: a matlab package for analysis and solution of discrete ill-posed problems[END_REF] and references therein, and this topic is somewhat beyond the scope of this work.

For automatic generation of patch layouts the interested reader is referred to[START_REF] Zheng | Boundary correspondence of planar domains for isogeometric analysis based on optimal mass transport[END_REF].

Note that while we here use tensor product spline spaces for the approximation of the solution to the PDE, adaptive methods that use local refinement exist, see for example[START_REF] Giannelli | Thb-splines: The truncated basis for hierarchical splines[END_REF].

In the IGA formulation[START_REF] Hinz | Elliptic grid generation techniques in the framework of isogeometric analysis applications[END_REF] we also integrate non-polynomials. We use 3 quadrature points for the mass matrix M and 7 quadrature points for the stiffness matrix K .

Specifically, the smallest spline space containing det J can be obtained by increasing the multiplicity of each knot in geometry map knotvector by (d -1) • p to account for the reduction in the differentiability and increase of the degree.

This could be done using either interpolation or L 2 projection. In either approach one has to solve a linear system to obtain the expansion coefficients. The linear system matrix needs to be inverted only once for a given spline basis and, moreover, the matrix has a Kronecker/separable structure (cf.[START_REF] Mantzaflaris | Low rank tensor methods in Galerkin-based isogeometric analysis[END_REF]) therefore the solution can be obtained extremely fast.

where u = (u 1 , . . . , u N ) T and K, M , T and f is given by

Note that the values of ε cr and µ r are set to the properties of gold in patch 4 and for air in the other patches. Owing to the restriction of the IGA library we are utilizing, we further reformulate the system of complex algebraic equations ( 19) as

.


Changing God's Law: The Dynamics of Middle Eastern Family Law Nadjma Yassari, ed.

Baudouin Dupret

This volume, edited by Nadjma Yassari, a leading authority in the field of Islamic and Middle Eastern law and especially family law, is the outcome of a conference held in Hamburg in 2013. It has the strengths and weaknesses of this type of publication. Regarding the latter, one can mention a certain degree of heterogeneity, which cannot be totally masked by the four-part structuration of the book. In addition, there is the very partial coverage of the very broad topic--change in family law-considered in an area, the Middle East, that is also broad and used in its widest sense, from Morocco to Pakistan. With regard to the strengths, one should emphasize the exceptional quality of the contributors, on the one hand, and the general framing of the question addressed by the volume, on the other. The latter follows the three pillars of the research group organized in the Hamburg Max-Planck Institute: interdisciplinarity, comparison, and attention to procedural law. As is often the case with collections, our review will proceed piecemeal.

In her introduction entitled "New family codes in Middle Eastern countries: reforms that are faithful to Islamic tradition?", Marie-Claire Foblets raises a series of relevant questions regarding, first, the emergence of what she calls a "new family law" and, second, the role of judicial interpretation. Foblets seems to imply that, initially, family codes did not depart from classical fiqh. On the contrary, there are good reasons to think that the adoption of a single code applicable to all Muslim citizens within the borders of one single (nation) state and enforceable through tightly organized jurisdictions structured in the three-level mode which is typical in civil law states did already constitute a legal revolution, in form as well as in content. This does not mean that codified family law was un-Islamic, but that its nature and economy was radically different from that of the normative organization that prevailed before the installment of the modern state. In that respect, the break is not between old and new family law, but between fiqh and positive law. Foblets gives the example of the reference to fiqh in case of gaps in legislation. However, looking closer at what judges do when filling these gaps, one observes that they either use templates or rely upon case law of their respective supreme court, made available on state-centralized databases. Actually, one can even contend that fiqh has become merely one source of law among many, albeit an important one. As a consequence, the conservative or progressist quality of national family laws is not primarily dependent on classical fiqh provisions, but on the social and political orientations of contemporary legislatures, and on the moral orientation of positive-law trained judges.

Part one of the volume consists of three chapters addressing the law, in general, and family law, in particular, at a very wide scale (the Middle East) or at a country level (Pakistan and Tunisia), under the quite broad umbrella of "continuities and changes". Chibli Mallat, in Chapter 2, extends his efforts to draw the history of "Middle Eastern law", taken as a millennium-old category embracing local and contingent developments such as Islamic law, to see, through case law rather than legislation, the main changes one can observe in the status of women after the 2011 upheavals. In legal history, it is rather well established that there are at best merely occasional connections between political and legal revolutions. Mallat explains these connections, in the frame of the Middle East, by the fact that equality already constitutes the horizon of legal reformism. From a distance, one might add that these upheavals, which failed to transform into revolutions, were a game of power much more than a struggle for law, something that proved even truer with the counter-coups that quickly followed. Shaheen Sardar Ali, in Chapter 3, addresses the issue of family law reform and the plurality of legalities in Pakistan. She rightly remarks that "family law is now being discussed as a multilayered and complex subject," and not as a mere continuation of classical fiqh. Ali ascribes the responsibility for gender inequality in Muslim contexts to "a moral and anthropological context that privileged male authority and the male voice" rather than to holy scriptures which, on the contrary, "ensured women a legal status and distinct (if not equal) rights in all spheres of life." This is an assertion that should be qualified, considering the deep heterogeneity of the Sunna and even the Qur'an. Nonetheless, she convincingly states that the nature of legal reforms was closely dependent on the political regime that brought it about. She also emphasizes that, whatever the conservative or progressist spirit of legal reform, it always took place in the frame of the tutelage/protection (wilâya/qiwâma) nexus. This is also called in many places the complementary rather than equal roles of men and women-and it is a frame that prevails among state officials, and especially the judges and their case law, and across the many strata of the population. This notwithstanding, she also shows that present law in Pakistan is the heir of Anglo-Muhammadan law more than classical fiqh. In Chapter 3, Monia Ben Jémia addresses the relationship between family law and human rights in posttransition Tunisia. Interestingly, she demonstrates that even though family law was reformed long ago in Tunisia, the 2011 revolution and especially the 2014 constitution gave a new impetus to its progressist interpretation and actual implementation. Her conclusion remains open: "it remains to be seen whether equality will persist as a pillar of the Tunisian legal system or whether the ambiguity of the constitutional provisions on equality in family relationships will pave the way for a reversal on the basis of religious reference."

Part two of the collection, addressing legislation, is also made of three chapters. In chapter 4, Lena-Maria Möller comparatively examines family laws in the Gulf area. She shows that Qatar, Bahrain and the United Arab Emirates started to codify issues related to the family with the dawning of the third millennium. She emphasizes that it constitutes a break, as codification is a state enterprise aiming to legislatively shape the law, something that strongly contrasts with the way classical Islamic doctrine was elaborated. Interestingly, she stresses the fact that Gulf laws were deeply influenced by Egyptian law, since most jurists in the region are originating from Sanhuri's homeland. In chapter 5, Nora Alim and Nadjma Yassari review the law-making process in Egypt in the context of family law, dividing it into three reformative methods: the codification of substantive family law, the inserting of stipulations in marriage contracts, and the development of procedural law. As to the third method, they relate it to both the French influence and the legacy of siyâsa justice, the norms applied by the ruler parallel to shari'a courts. Such procedural method is particularly useful as it allows the reformation of the law without explicitly touching its content. The authors convincingly show that this approach put the reform of family law in Egypt in a stalemate from which one can see no way out unless further reforms in substantive law are achieved. In chapter 6, Nadjma Yassari explores the financial relationship between spouses in Iranian law. As is the case with the chapter on Gulf laws, Yassari's contribution offers an in-depth study of a topic which is weakly addressed in most literature. Sketching the general evolution of Iranian family law, she shows how the postrevolutionary legislature strengthened the financial situation of divorced women "by reforming existing instruments and introducing new concepts." However, such strengthening reflects a gendered conception of the society, the family, and the role of women, which does not suit the 21 st century Iranian context, thereby risking the disconnection of the law from the society and running against the conception prevailing among young families that marriage is "a partnership that the spouses shape together."

Part three addresses the judiciary at work in matters related to the family. Chapter 7, in French, deals with Tunisia and the powers of judges in matters of family law. Salma Abida stresses that, although the legislative texts are liberal and modernist, after the post-independence adoption of the family code and the unification of the judiciary, their implementation and thus their impact is largely dependent on the judges, especially since the Tunisian legislature issued ambiguous provisions leaving room to judicial discretion. As a whole-and as could be expected-the courts' position has fluctuated between conservative and progressist approaches. Nathalie Bernard-Maugiron's contribution, in chapter 8, is a thoroughly documented attempt to describe judges at work in matters of divorce in Egypt. As Abida in the former chapter, Bernard-Maugiron emphasizes the fact that no legislative reform can be really assessed without a close look at its implementation and especially at judges' work of interpretation. After having quickly reviewed what she calls "law in the books", i.e. statutory law, she proceeds to the examination of case law in the specific domains of divorce on the ground of moral and physical harm and of divorce for polygamy. She demonstrates that reforms in divorce, while modest in their scope, substantially improved the legal status of women and their capacity to free themselves from unhappy marriages. However, even more interestingly, she shows that the role of progressist legal reforms and open-minded judges must be assessed against the background of a society within which the women's exercise of their rights is socially despised and faces therefore severe limitation. Chapter 9 addresses the practice of personal status law in Israel. It can serve as a backdrop to the other studies collected in this volume. Imen Gallala-Arndt offers an interesting summary of Israel's legal history, giving some religion related clues as to why Israel has not adopted any written constitution so far. She also proposes an overview of tribunals having jurisdiction in family matters, a legal domain in which the country follows the principle of the personality of laws. Finally, she describes some of the conflicts arising between religious and civil courts, especially the Supreme Court. In conclusion, she argues that the civil courts have a corrective impact on the personal status system and religion-based discriminatory rules. However, she also shows that this corrective action had an influence on the perpetuation of the religion-dominated personal status system.

The fourth and last part of the volume deals with the concept of "party autonomy", defined in the preface as the capacity of spouses to "abrogate dispositive marital law through marital agreements and thus reach tailored contractual agreements." Chapter 10 draws a grand picture of marriage contracts in Islamic history. Amira Sonbol pointedly remarks that it is paradoxical that any change in the structure of the family in its modern and thus state-centered meaning is perceived as an attack against Islam. Actually, this paradox perfectly reflects the transformation of fiqh into positive law. One can have more reservations regarding the picture she paints of pre-Islamic marriage or of marriage in Islam, as it bears the features of strong essentialism, and one can even doubt the possibility of drawing such picture. What is said about modern times should have been taken into consideration regarding ancient times. Specifically, concepts such as marriage have an history, they appear, disappear, live, and transform; and their ontology is historical. The state has transformed the understanding of marriage and family in a radical way, but probably not more radical than, for instance, the emergence of fiqh literature or the Ottoman empire's Hanafi-dominated judicial system. In chapter 11, M. Siraj Sait examines the issue of matrimonial property in Islamic contexts, whose systematically-claimed personal character contrasts with the paucity of legal sources on the subject. Addressing legal reform which, despite strong social and institutional resistances, point to equitable matrimonial property regimes, his contribution challenges the perceived incompatibility of classical Islamic law with concepts of marital property and the supposed inability of Muslim societies to recognize the contribution of both spouses.

In sum, one could easily subscribe to M. Siraj Sait's concluding quotation, according to which there is "no room for doubt that a common ground could be found between Islamic law and gender equality," if such statement was not marred by an essentialist posture. Islam, the Sharia, fiqh, Islamic law, and related terms have no essential meaning, at least in the perspective of the legal, social, human and historical sciences, and it is precisely for this reason that they can be made compatible with issues such as gender equality, democracy, human rights, and the like. This volume is valuable for the extended and documented information it provides the reader with. One can find in it much material that is needed to produce the comprehensive analysis its heterogeneous character could not achieve.


Farah Gabsi
email: farah.gabsi@univ-lorraine.fr

Frédéric Hamelin

Nathalie Sauer
email: nathalie.sauer@univ-lorraine.fr

Joseph Julien Yamé

Frederic Hamelin
email: frederic.hamelin@univ-lorraine.fr

Joseph J Yame
email: joseph.yame@univ-lorraine.fr

ℓ_ϵ-regularized economic model predictive

Keywords: Energy-efficient buildings, Model Predictive Control, Regularization

come

INTRODUCTION

In the context of intelligent buildings, modern centralised automation systems are often used to improve their energy efficiency. "Building Automation and Control Systems" (BACS) are generally based on a dynamic model of buildings. Depending on their complexity and/or performance, they may also include a precise description of the most energyintensive equipment (heating, ventilation and air conditioning (HVAC) systems [START_REF] Rawlings | Economic MPC and real-time decision making with application to large-scale HVAC energy systems[END_REF],...), the price of electricity or the behaviour of occupants. In addition to this optimized energy management, thermal comfort inside the building is usually a factor taken into consideration, which leads to a global control problem [START_REF] Gabsi | Hygrothermal modelling and MPC-based control for energy and comfort management in buildings[END_REF].

Model Predictive Control (MPC) is one of the most used advanced control strategies in this context, mainly due to its ability to achieve economic objectives, taking into account a simplified dynamic model and different constraints [START_REF] Godina | Model predictive control home energy management and optimization strategy with demand response[END_REF], [START_REF] Serale | Model predictive control (MPC) for enhancing building and HVAC system energy efficiency: Problem formulation, applications and opportunities[END_REF]. The modelling method influences the actual practice of MPC in buildings because of its cost and scalability [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF], (Gabsi et al., 2018a), [START_REF] Zhuang | A new simplified modeling method for model predictive control in a medium-sized commercial building: A case study[END_REF]. Economic Model Predictive Control (EMPC) [START_REF] Zong | Challenges of implementing economic model predictive control strategy for buildings interacting with smart energy systems[END_REF], [START_REF] Rawlings | Economic MPC and real-time decision making with application to large-scale HVAC energy systems[END_REF], [START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF] is becoming increasingly popular because of its interest in considering more general economic cost functions than traditional quadratic cost functions.

In recent years, the theory of LASSO (Least Absolute Selection and Shrinkage Operator), particularly used in signal processing, has led to the emergence of new predictive control strategies called " asso MPC" [START_REF] Gallieri | l asso MPC: Smart regulation of over-actuated systems[END_REF], [START_REF] Rao | Sparsity of linear discrete-time optimal control problems with l 1 objectives[END_REF] or RMPC for "Regularized MPC" [START_REF] Amy | Regularized MPC for power management of hybrid energy storage systems with applications in electric vehicles[END_REF]. By using penalty criteria in 1 -norm that favor some kinds of sparse controls, it becomes possible, for example, to limit the number of active control inputs in an over-actuated system [START_REF] Gallieri | Model predictive control with prioritised actuators[END_REF] or to prioritize actuator actions and efficiently distribute control effort [START_REF] Amy | Regularized MPC for power management of hybrid energy storage systems with applications in electric vehicles[END_REF]. It is also possible to consider certain control applications that require the use of piecewise constant or impulsetype control signals, with as few changes as possible [START_REF] Pakazad | Sparse control using sum-of-norms regularized model predictive control[END_REF]. In the same way, a binary regularization term can be introduced in order to penalize differently the power variations of actuators depend-ing on whether they are in normal operation, startup or shutdown [START_REF] Cojocaru | A binary-regularization-based model predictive control applied to generation scheduling in concentrating solar power plants[END_REF]. Finally, RMPC is also relevant for reducing data packet size in Network Control Systems (NCSs) [START_REF] Nagahara | Sparse packetized predictive control for networked control over erasure channels[END_REF].

These different problems can also be solved by considering penalties 0 [START_REF] Aguilera | Quadratic model predictive control including input cardinality constraints[END_REF], [START_REF] Aguilera | Quadratic MPC with l 0 -input constraint[END_REF], which still have a better parsimonious capacity but make the optimization problem NP-Hard because non-convex. To obtain a continuous but still relatively sparse control, some works use the CLOT (Combined L-One and Two) norm [START_REF] Challapalli | Continuous hands-off control by clot norm minimization[END_REF], which is a convex combination of 1 and 2 -norms and thus allows to benefit from the advantages of each of them.

In this paper, a new predictive control strategy regularized by ε -norm penalties is presented. By a judicious choice of the regularization terms, this approach allows in particular to control the solicitations of certain equipments (HVAC,...), for which too frequent starts/stops are critical and are the most energy inefficient method of operating. It also makes it possible to control the number of active control inputs at any time or to avoid inappropriate control scenarios.

The paper is organized as follows. Section 2 first specifies the objectives of the MPC by defining a criterion combining both economic and thermal comfort aspects. The constraints for using different conventional equipments are also specified. Section 3 makes these various objectives explicit in the form of a regularized functional. The proposed strategy is applied to the thermal regulation of buildings in section 4. In particular, CRAN's "Eco-safe" platform is used to highlight the practical value of the proposed approach. Finally, a conclusion and perspectives are presented in Section 5.

PROBLEM STATEMENT

This section first defines a criterion for thermal comfort in a multizone building. To satisfy such a criterion while minimizing the energy consumed, a functional is defined as part of the synthesis of a predictive control. Constraints linked to the number of actuators used at any time as well as the variability of the control scenarios are added in order to optimize, among other things, the lifetime of the equipment (heat pump (HP) systems, double flow controlled mechanical ventilation (CMV),...).

Thermal Comfort

We consider an air volume z i (hereinafter referred to as zone z i ) delimited by n surfaces (walls, windows, ceiling, floor). Thermal comfort is ensured within this zone if the operative temperature T O p,z i belongs to a comfort temperature range defined by T C ±2 K, where K denotes Kelvin degree.

The operative temperature T O p,z i can be approached by:

T O p,z i ≈ T MR,z i + T z i 2 (1)
where T z i is the ambient air temperature in zone z i and T MR,z i is the mean radiant temperature defined by:

T MR,z i = ∑ n j=1 S z i ,s j × T z i ,s j ∑ n j=1 S z i ,s j (2)
with T z i ,s = T z i ,s j 1≤ j≤n representing the temperature of each surface in contact with zone z i . The contact surface area is assumed to be S z i ,s j .

As for the comfort temperature, [START_REF] Mccartney | Developing an adaptive control algorithm for europe[END_REF] determines it on the basis of studies carried out in situ in buildings. It is a simple linear regression model that fits the filtered temperature T RM of the outside air:

T C = 0.049 T RM + 9.2 if T RM ≤ 283.15 K 0.206 T RM -34.85 if T RM > 283.15 K (3)
with T RM a temperature that changes daily (D) according to the average outdoor temperature T DM of the previous day (D -1):

T RM (D) = 0.8T RM (D -1) + 0.2T DM (D -1) (4)

Thermal Model

Before defining the cost function and control constraints, it is necessary to determine a dynamic model reflecting the thermal behaviour of the building by integrating the various equipment as well as all influential disturbances. According to (Gabsi et al., 2018a), the dynamic thermal behavior of a zone z i delimited by n surfaces (Σ j ) can be represented by the following descriptor time-varying discrete-time system with regular pencil:

                     E z i x z i (k + 1) = A z i x z i (k) + F z i T z i (k) + q z i (k)+ b Sol z i q Sol z i (k) 1 -δu VB z i (k) + b TTW z i (T Ext (k))u TTW z i (k)+ b HP z i (T HP ,φ HP )u HP z i (k)+b CMV z i (T Ext (k),T z i (k))u CMV z i (k)+ . . . (if other equipment is to be considered) T z i T z i ,s (k) = C z i x z i (k)
(5) with:

• x T z i (k) = T z i , T T s i (k);
• T s i (k): the core temperature of each of the surfaces surrounding zone z i ; • T z i (k): the ambient air temperature in each of the zones adjacent to zone z i (possibly including the outside air temperature); • q z i (k): the algebraic value of all incoming and outgoing heat fluxes in z i ; specifically, q Sol z i (k) reflects incoming short-waves solar radiation.

• u HP z i (k), u CMV z i (k), u TTW z i (k), u VB z i (k), . . .: the control inputs in z i associated with the start-up of HP, CMV, tilt-turn windows (TTW), venetian blinds (VB), . . . • 1δu VB z i (k): the solar heat gain coefficient (SHGC [START_REF] Cho | Solar heat gain coefficient analysis of a slim-type double skin window system: Using an experimental and a simulation method[END_REF]) with 0 ≤ δ ≤ 1; u VB z i depends on whether the blinds are raised (u VB z i = 0) or lowered (u VB z i = 1); • T Ext (k): the outdoor temperature; • T HP (vs. φ HP ): the air temperature (vs. speed) at the heating/cooling system's supply line. When a building consists of N contiguous zones z i , the models (5) of each zone can be aggregated, which leads us to consider the time-varying discrete-time system defined as follows:

                   x z (k + 1) = A z x z (k)+ F T T Ext (k)+ F q q z (k)+ N S ∑ ξ=1 B ξ x z (k), q Sol z (k), T Ext (k) u ξ (k)   T z T z,s T Ext,s   (k) = Cx z (k) (6)
with:

• T z (k) = T z i (k) 1≤i≤N ; • T z,s (k) = T z i ,s (k) 1≤i≤N ;
• u(k) = u ξ (k) 1≤ξ≤N S : the control vector defined from the control inputs of all zones z i .

It reflects all possible control scenarios within the multi-zone building (no action, HP (on/off), CMV (on/off), automatic tilt-turn windows (open/close), venetian blinds (open/close),...). Each u ξ (k) element is equal to either 0 or 1, these two values corresponding respectively to the switching off or switching on of the ξth control input;

• x T z (k) = T T z , T T s (k): the state vector defined from vectors x z i (k) and T s (k) = T s i (k) 1≤i≤N ;
• q z (k): the disturbance vector grouping all heat fluxes into and out of the N zones z i ; the incoming short-waves solar radiation is specifically reflected by q Sol z (k). As an example, the diagram in Fig. 1 represents all the temperatures involved in the thermal modelling of a two-zone building. This simplified sketch, which shows a 2D horizontal cross-sectional view, assumes zones with no ceiling and no floor. 23 temperature nodes are useful for modeling this building:

T z ∈ R 2 , T z,s ∈ R 8 , T Ext,s ∈ R 6 , T s ∈ R 7 .

EMPC -Economic Cost function and Constraints

The principle of predictive control [START_REF] Rockett | Model-predictive control for non-domestic buildings: a critical review and prospects[END_REF] is to optimize a cost function to describe the control objectives over a forecast time horizon N P . At each instant k, an optimal control sequence {u * (k + j))} 1≤ j≤N P is calculated to minimize this function and only the first element u * (k + 1) is applied to the system. The economic objective function J MPC (u, x) that we propose in the context of the thermal regulation of a multizone building is as follows:

J MPC (u, x) = min u N P ∑ j=1 T O p (k + j) -T C (k + j) 2 Ψ j + u(k + j) 2 e u j + ∆u(k + j) 2 e ∆u j (7)
with:

• T T C (k + j) = T C,z 1 (k + j) . . . T C,z N (k + j) ∈ R N and T T O p (k+j)= T O p,z 1 (k+j). . . T O p,z N (k+ j) ∈ R N the estimated operative temperature in each zone z i according to model (6); • T O p (k + j) -T C (k + j) 2 Ψ j = N ∑ i=1 ψ z i , j |T O p,z i (k + j) -T C,z i (k + j)| 2 • u(k + j) 2 e u j = N ∑ i=1 N S ∑ ξ=1 e ξ,z i , j u ξ,z i (k + j) 2 ; • ∆u(k + j) 2 e ∆u j = N ∑ i=1 N S ∑ ξ=1 e ξ,z i , j u ξ,z i (k + j) -u ξ,z i (k + j -1)
• According to models ( 5) and ( 6), all control inputs can be grouped according to the vectors u z i (k), u(k + j) and u defined as follows:

     u T z i (k) = u HP z i (k) u CMV z i (k) u TTW z i (k) u VB z i (k) u T (k + j) = u T z 1 (k + j) . . . u T z N (k + j) u T = u T (k + 1) . . . u T (k + N P ) (8)
• Weightings e ξ,z i , j and e ξ,z i , j make it possible to specify the cost in euros of each of the possible actions on the system. The e-terms will reflect the energy cost of starting a HP or a VMC for a period of time while the e -terms will reflect the opening or closing of tilt-turn windows or venetian blinds.

As for ψ z i , j , it reflects the importance attached to

T O p,z i (k + j) being close to T C,z i (k + j).
Remark 1: The minimization of J MPC (u, x) (7) requires the calculation of the operative temperature T O p in each zone z i . According to the definition (1), this calculation requires not only the knowledge or estimation of the temperature of each surface in contact with these zones, but also their prediction over the forecast time horizon N p . In a properly insulated building, since the ambient air temperature in each zone z i is generally not very different from the mean radiant temperature, it is common to consider T z rather than T O p in the definition of criterion J MPC (u, x).

In the long term, the economic cost of a control scenario is not reduced to the mere expression of criterion (7). Indeed, the latter does not take into account abnormal wear and tear or premature failure of actuators, which can generate significant additional costs.

In the building context, for example, it is recognized that increasing the number of on/off cycles of a HP compressor increases its electrical consumption but also its wear and tear. Therefore, when synthesizing the control scenario {u * (k + j)} 1≤ j≤N P over a horizon N P , it is important to reduce as much as possible this number of on/off cycles for this equipment.

In order to minimize criterion J MPC (u, x) while taking into account this last remark, a regularisation term is added.

REGULARIZED EMPC (REMPC)

Regularized Criterion

It is proposed to consider a regularized criterion such as:

J λ,Ω (u) = (1 -λ)J MPC (u, x) + λΩ (u) (9)
with J MPC (u, x) defined by ( 7) and where 0 ≤ λ ≤ 1 is a regularization parameter. The additional term λΩ(u) in the criterion amounts to regularizing the solution through a penalty of the latter. In the context of the problem presented before, λΩ (u) must be a penalty that favours the parsimony of the first derivative of u. Based on Tikhonov's regularization method [START_REF] Engl | Regularization of inverse problems[END_REF], one possible technique is to include a linear operator R in the regularization term Ω (u), and solve the following problem :

u * λ,p =arg min u∈R nu (1-λ)J MPC (u, x)+λ Ru min(1,p) p ( 10
)
where • p is the p -norm w i p :=

n w i ∑ i=1 |w i | p 1 p
for a vector w i ∈ R n w i . The power min(1, p) of the regularization term Ru min(1,p) p makes it possible to consider by continuity the ∞ -norm of Ru for p → ∞:

lim p→∞ Ru min(1,p) p = lim p→∞ ∑ i |(Ru) i | p 1 p = Ru ∞ (11)
and the 0 -norm of Ru for p → 0:

lim p→0 Ru min(1,p) p = lim p→0 ∑ i |(Ru) i | p = Ru 0 (12)
Equation ( 12) is particularly interesting because, from a theoretical point of view, the parsimony of Ru is measured using its 0 -norm corresponding to the total number of non-zero elements:

Ru 0 = # (i| (Ru) i = 0) (13)
The linear transform R can take different forms:

• 0 th -order regularization favouring solutions with a small norm:

R = R 0 = I (14)
• 1 st -order regularization. It consists in focusing a priori on the low oscillating nature of the solution, and thus penalizing rapid variations:

R = R 1 =         -1 1 0 0 0 0 0 0 -1 1         (15)
We notice that at 0 th -order, the product Ru represents a discretization of the vector u, while at 1 st -order it is a discretization of its first derivative.

p Penalization

There are several types of penalty functions [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF]. The Ridge regression corresponds to a penalty of type 2 -norm. As we will see below through a simple example, this function has the particularity of not cancelling the coefficients of Ru but rather reducing them and making them tend towards 0. This is a "shrinkage" of coefficients. The Lasso regression, introduced by Tibshirani [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], is a regression penalized by the 1 -norm of the coefficients of Ru, which favours parsimony. Fused-Lasso is a variant which allows to take into account the spatiality of the variables [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]. The objective is to have close estimates for the same variable when they are "close in time". This is made possible by penalizing the 1 -norm of the difference of this variable in two successive instants.

An even more natural penalization than Ru 1 is to consider a constraint Ru ε ε (with 0 ≤ ε 1), which not only contracts the value of the different elements of Ru but also forces certain elements u i to be strictly zero for λ large enough thanks to the shape of the isolines of Ru ε ε . By way of illustration, we consider the criterion:

J λ,p (u) = (1 -λ)J 1 (u) + λ Ru min(1,p) p (16) with R = I (14) and J 1 (u) = (u 1 -2 + u 2 ) 2 + (u 2 -0.5) 2 .
First, we observe that ∀p ≥ 0, we have: p Between these two extreme values, the trajectory of Γ λ,p followed by the minimum u * λ,p of J λ,p (u) as a function of λ is represented in the u 1u 2 plane in red in Fig. 2. The ellipsoids and the filled contour plot in the background of these figures are isolines of J 1 (u) and Ru min(1,p) p respectively. We can see on each of the subfigures (Fig. 2a-2d) that the shape of the trajectories Γ λ,p is very different according to the values of p.

In particular, the LASSO selection (Fig. 2c) results in a more parsimonious solution than the Ridge selection (Fig. 2b), which tends to make the coefficients very small without cancelling them. More generally, for p > 1, the trajectory Γ λ,p tends towards u * ∞,p for λ increasing but without being "attracted" by the axes u 1 = 0 and u 2 = 0. On the other hand, as soon as p ≤ 1, we observe that this convergence towards u * ∞,p takes place along one of the axes u 1 = 0 or u 2 = 0. In this context (p ≤ 1), it is interesting to point out the difference in result between a regularization term in ε -norm and in 1 -norm. It appears through this example that the axis of attraction can be different (u 2 = 0 for 1 -norm (fig. 3c) vs. u 1 = 0 for 0.2 -norm (fig. 3d)). Given the shape of the isolines of J 1 (u), the preferred solution (in the sense of minimizing the criterion J 1 (u)) is the one associated with the lowest p-value. Another phenomenon appears for p < 1; the trajectory Γ λ,p followed by u * λ,p becomes discontinuous as the value of p decreases. This is due to the non-convexity of the term Ru p p which increases for small values of p and makes the criterion J λ,p (u) non-convex.

APPLICATION TO THERMAL COMFORT IN MULTIZONE BUILDINGS 4.1 Regularized Criterion

With regards to the regularization term Ω (u), the following specifications are formulated:

• the number of on/off cycles for HP and CMV in all zones z i should be limited. In view of the previous paragraph, the minimisation of norms

N ∑ i=1 R 1 u HP z i ε ε and N ∑ i=1 R 1 u CMV z i ε
ε reflects this dual objective with:

u HP z i = u HP z i (k) | ûHP z i (k + 1) . . . ûHP z i (k + N P ) T u CMV z i = u CMV z i (k)| ûCMV z i (k+1). . . ûCMV z i (k+N P ) T
• the number of effective actuators should be limited at each time of prediction k + j and in each zone z i . This specification translates as minimizing the term

N ∑ i=1 R 0 u z i (k + j) ε
ε with u z i (k + j) the control vector defined by (8);

• the control inputs take values in the discrete set {0, 1} (on/off or open/closed). The minimization of the norm R 0 (u (u -1)) ε ε meets this objective by defining by u (u -1) the element-wise product of two vectors u and (u -1). For illustration purposes, the trajectory of Γ λ,p followed by the minimum u * λ,p of J λ,p (u) = (1λ)J 1 (u) + λ (u (u -1)) min(1,p) p as a function of λ is represented in the u 1u 2 plane in red in Fig. 3. For this example, we observe that Γ λ,p tends towards 1 1 T = on on T , which is the on/off control associated with the lowest value of J 1 (u). For this type of regularization, the choice of a ε -norm is also justified because the appearance of the isolines J 1 (u) is not modified by considering J λ,0 (u) knowing that u (u -1) 0 = dim(u) = 2 for all u except for axes u 1 = 0 and u 2 = 0. This property is interesting because the sole purpose of this regularization term is to impose discrete values on u and not to modify the shape of the isolines;

Therefore, Ω (u) is defined as:

Ω (u) = α 1 R 0 NB i=1 (u -u i 1) ε ε + N ∑ i=1 α 2 R 1 u HP z i ε ε + α 3 R 1 u CMV z i ε ε +α 4 N ∑ i=1 N P ∑ j=1 R 0 u z i (k + j) ε ε (17)
with:

• u and u z i (k + j): augmented control vectors (8);

• NB i=1 (u -u i 1) = (u -u 1 1) . . . (u -u NB 1)
where is the element-wise product of two vectors and (u i ) 1≤i≤NB represent some constant values that can be taken by at least one of the control inputs. For on/off controls, we have NB = 2, u 1 = 1 (on) and u 2 = 0 (off);

• u TTW z i = u TTW z i (k) | ûTTW z i (k+1) . . . ûTTW z i (k+N P ) T ; • u VB z i = u VB z i (k) | ûVB z i (k+1) . . . ûVB z i (k+N P ) T .
An α 2 (or α 3 ) that tends towards 0 will generally cause permanent stress on the associated actuator, which will affect its lifetime. Conversely, an α 2 (or α 3 ) that tends towards infinity ensures a low stress on the actuator but thermal comfort performance can them become poor according to criterion J MPC (u, x) because the solution obtained then becomes too far from the optimal non-regularized solution.

Minimisation of the regularized criterion

The closed-loop REMPC requires an optimal solution to (9) at each step. It is usually very difficult to quickly find optimal solutions (and prove their optimality) for non-convex problems. First of all, it is important to mention that the minimization of the non-regularized economic cost function J MPC (u, x) constrained by the time-varying model ( 6) does not present an analytical solution mainly because of the input matrix B ξ x z (k), q Sol z (k), T Ext (k) , which depends on the chosen control scenario. Thus, even if the non-convexity of the ε -norm makes NP-Hard the regularized optimization problem, the complexity of the solving algorithm is not significantly increased with the regularization term. Although solvers such as "Interior Point OPTimizer" (IPOPT) [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] can efficiently find local solutions to nonlinear programming problems, these solutions are not particularly suited to our problem.

To solve this problem, we use a recently developed iterative algorithm [START_REF] Gabsi | Hygrothermal modelling and MPC-based control for energy and comfort management in buildings[END_REF], the objective of which is to estimate the optimal control scenario with a controlled computation load. The idea is to keep at each time of prediction k + j only a limited number of scenarios among all those that are possible. For this purpose, main component analysis is performed on a limited number of points judiciously chosen in the variables space (T O p,z 1 (k + j), T O p,z 2 (k + j), . . . , T O p,z N (k + j)) in order to determine a suitable basis for the representation of all possible realizations of T O p,z i (k + j) with i = 1, . . . , N. By prioritizing the information, this makes it possible to replace the set of all these realizations by a smaller subset S j whose cardinality is set a priori. This procedure is repeated N p times in an iterative manner for j ranging from 1 to N p . Of course, the larger the number |S j |, the better the approximation of the optimal control scenario. This prioritization technique allows to consider not only a non-convex regularized criterion J λ,Ω (u) (9) but also important values for the forecast time horizon N P .

Application to the "Eco-Safe" Platform

The "Eco-Safe" platform consists of six zones and a corridor (Fig. 4). Each zone z i = i is equipped with several actuators and sensors to ensure a certain level of thermal comfort. The "research and development" room (R&D, zone z 1 ) and the "handling" room (zone z 5 ) are frequently occupied by students. They have a respective surface area of 51 m 2 and 34 m 2 . Four other rooms, of 17 m 2 each, are used for the storage of different materials. The temperature of zones z 1 and z 5 can be controlled by two reversible air-to-air HPs that can generate hot and cold air. A CMV also allows air exchange between these two zones. A weather station integrating several sensors communicates with the platform and allows to know the different characteristics of the outside air. To illustrate the predictive control strategy presented above, we used a dynamic model of this platform established in [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF]. The latter integrates the thermal behaviour of all partitions/walls/windows in the zones as well as the different energy sources (solar radiation, occupants in the premises, HPs). It is possible to define 9 control scenarios:

• Scenario 0: no action on the system;

• Scenario 1: the CMV is switched on to circulate air from zone z 5 to zone z 1 ;

• Scenario 2: the CMV is switched on to circulate air from zone z 1 to zone z 5 ;

• Scenario 3: a HP is switched on in cooling mode (constant air flow, outlet temperature equal to 291 K) in zone z 1 ;

• Scenario 4: a HP is switched on in cooling mode in zone z 5 ;

• Scenario 5: both HPs are switched on in cooling mode in zones z 1 and z 5 .

• Scenario 6: a HP is switched on in heating mode (constant air flow, outlet temperature equal to 313 K) in zone z 1 ;

• Scenario 7: a HP is switched on in heating mode in zone z 5 ;

• Scenario 8: both HPs are switched on in heating mode in zones z 1 and z 5 .

The state-space model ( 6) associated with the plat-form is:

         x z (k + 1) = A z x z (k) + F T T Ext (k) + F q q z (k)+ 8 ∑ ξ=0 B ξ (x z (k), T Ext (k)) u ξ (k) T z (k) = Cx z (k) (18) with: • T z (k) = T z 1 , T z 2 , T z 4 , T z 5 , T z 7 T (k); no surface temperature is measured in zones z i ;
• u(k)= u ξ (k) 0≤ξ≤8 : the control vector formed of 0 and 1. It is defined as u ξ (k) = 1 and u ξ (k) = 0 for all ξ = ξ when scenario u ξ is implemented at time k;

• q z (k) = T Ext , Occ, Sol West , Sol East , T HP T (k);

• T Ext (k): the outside temperature;

• Occ(k): the number of occupants in zone z 1 ;

• Sol West (k) and Sol East (k): the energy provided by solar radiation for zones z 1 , z 2 , z 3 and z 4 , z 5 , z 6 resp.;

• T HP (k): the temperature of the air forced through the HPs, T HP (k) = 291 K (cooling mode) or 313 K (heating mode).

Since the size of state vector x z is very large (72dimensional), the calculations required to develop a predictive control would become very complicated and time-consuming. The definition of a reduced order model is therefore necessary, based on a balanced state-space realization. In addition, as the purpose of this application example is mainly to demonstrate the interest of the regularisation terms, for the sake of simplicity and in accordance with Remark 1, the ambient air temperature in each zone z i is considered into the economic cost function J MPC (u, x) (7), instead of the operative temperature:

J MPC (u, x)= min u N P =12 ∑ j=1 T z (k+j)-293 2 Ψ j + u(k+j) 2 e u j
(19) The reduced state vector is defined to correspond to the ambient air temperatures measured in each zone z 1 , z 2 , z 4 , z 5 and z 7 . Horizon N p is chosen equal to 12 due to the sampling period which is 5mn. The forecast time horizon is therefore one hour depending on the platform's heating and cooling capacities. The weighting matrix Ψ j is diagonal and time-variant. The ith term of diagonal (Ψ j ) i,i is zero if zone z i is unoccupied at the time of prediction k + j. Otherwise, this term, which is chosen equal to ψ z i > 0, makes it possible to give more or less importance to the energy criterion (second part of ( 19)) compared to the performance criterion (first part of ( 19)). For this illustration example, we will assume that only zone z 1 is occupied during standard working hours. The term ∆u(k + j) 2 e ∆u j is not taken into account in the relationship (19) because the platform does not have automated closing/opening of tilt and turn windows or venetian blinds. As for the diagonal matrix e u j , it is defined according to the energy consumption of each of the scenarios (which is assumed to be constant no matter the day and time), namely:

e u j =
1 12 ×diag 0 30 30 10 3 10 3 2000 10 3 10 3 2000

The regularisation term Ω (u) (17) takes the following form:

Ω(u)=α 1 R 0 u (u-1) ε ε +α 2 R 1 u HP z 1 ε ε + R 1 u HP z 5 ε ε +α 2 R 1 u CMV z 1 ε ε + R 1 u CMV z 5 ε ε (20)
with a very large value for α 1 . Regularization parameter λ (9) is taken as 0.5.

Fig. 7a to 7d represent over the period from 18 to 20 February 2019 the evolution of the ambient air temperatures T z 1 (k) and T z 5 (k) respectively associated with the control scenarios in Fig. 8a to 8d. The cumulative energy cost (in kWh) for each of these scenarios is shown in Fig. 9a to 9d respectively. The weather conditions present during this period at the platform location (Nancy, France) are reflected in Fig. 5a, 5b and 6a. The first two show the daily evolution of the solar energy entering the platform zones according to their orientation. The third figure shows the evolution of the outside temperature in K. Fig. 6b finally shows the number of occupants in zone z 1 (we will assume that each person emits 80W of internal heat gain). Fig. 8 represents as a function of time the value of the index ξ associated with the non-zero element of the vector u(k). A blue cross on a line ξ of these figures at time k therefore reflects the implementation of scenario ξ on the platform (u ξ (k) = 1). For example, Fig. 8c shows that the non-zero element of u(k) over the entire first day is u 0 (k). This control scenario is associated with parameters ψ z 1 = 1 and α 2 = 3. A simple analysis of these different figures shows:

• the influence of parameter ψ z 1 . The larger this parameter is compared to the terms of e u j , the greater the proportion of criterion J MPC (u, x) related to comfort performance increases at the expense of energy consumption. However, the lower this parameter is, the more important the energy aspect becomes compared to thermal comfort. • the interest of the regularization term related to α 2 . By comparing Fig. 8a and 8b with Fig. 8c and 8d respectively, it is clear that an increase in α 2 significantly reduces the number of switching (on/off) of VMCs and HPs. A direct consequence is the almost complete disappearance of the oscillations observed on T z 1 (k) and T z 5 (k) (Fig. 7d vs. Fig. 7b).

• the interest of the regularization term in α 1 which allows to privilege discontinuous optimal values for u during the optimization of the regularized criterion. For this illustration example, only values 0 and 1 are allowed for all elements of u.

CONCLUSION

This paper aims to optimize the energy efficiency of multizone buildings by implementing a regularized economic model predictive controller (REMPC). More precisely, the objective is to maintain thermal comfort in occupied zones while minimizing energy consumption.

To achieve this long-term overall objective, an economic cost function was first defined and control specifications were added via a regularization criterion. The first regularized term concerns the limitation of the frequency of on/off cycles for particular actuators such as HPs or VMCs. Indeed, the very frequent starting and stopping of this type of equipment makes inefficient their energy operating mode and can especially lead to severe damages. Other criteria related to the number of actuators used at a given time or to the discrete nature of certain control variables were also taken into account in the regularized criterion. An analysis showed the importance of choosing an appropriate p -norm to define these regularization terms. It has been shown that a ε -norm is to be preferred. This control strategy was tested on a platform simulator [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF] located in the CRAN laboratory (Nancy/France) and gave very encouraging results for on-site implementation.

Figure 1 :

 1 Figure 1: Thermal modeling of a two-zone building.

Figure 2 :

 2 Figure 2: Trajectory Γ λ,p (in red) followed by u * λ,p as a fonction of λ in the u 1u 2 plane.

Figure 3 :

 3 Figure 3: Trajectory Γ λ,p (in red) followed by u * λ,p as a fonction of λ in the u 1u 2 plane.

Figure 4 :

 4 Figure 4: The "Eco-Safe" platform.

  Sol West (k); z i = z 1 , z 2 , z 3 . Sol East (k); z i = z 4 , z 5 , z 6 .

Figure 5 :

 5 Figure 5: Energy (W/m 2 ) provided by solar radiation for zones z i .

  Outside temperature T Ext (K).

Figure 6 :

 6 Figure 6: Disturbances q(k) (in addition to solar radiation).

  ψ z 1 = 20, α 2 = 10.

Figure 7 :

 7 Figure7: Ambient air temperatures T z 1 (K, in solid blue line) and T z 5 (K, in red dashed line).

  (a) ψ z 1 = 1, α 2 = 0. (b) ψ z 1 = 20, α 2 = 0. (c) ψ z 1 = 1, α 2 = 3. (d) ψ z 1 = 20, α 2 = 10.

Figure 8 :

 8 Figure 8: Control scenarios.

  ψ z 1 = 1, α 2 = 3.

  ψ z 1 = 20, α 2 = 10.

Figure 9 :

 9 Figure 9: Cumulative energy cost (kWh).

ACKNOWLEDGEMENTS

This work has financial support from the Contrat de Plan Etat-Région (CPER) 2015-2020, project "Matériaux, Energie, Procédés".