Martin Bromberger

Mathias Fleury
email: mfleury@mpi-inf.mpg.de

Simon Schwarz
email: sschwarz@mpi-inf.mpg.de

Christoph Weidenbach

SPASS-SATT: A CDCL(LA) Solver

Keywords: Linear Arithmetic, Integer Arithmetic, SMT, Preprocessing

come

Introduction

SPASS-SATT (v1.1) is a sound and complete CDCL(LA) solver for quantifierfree linear rational and linear mixed/integer arithmetic. It is a from-scratch implementation except for some basic data structures taken from the SPASS [START_REF] Weidenbach | SPASS version 3.5[END_REF] superposition theorem prover. It is available through the SPASS-Workbench [START_REF]The SPASS Workbench[END_REF]. We participated with SPASS-SATT in the main track of the 13th International Satisfiability Modulo Theories Competition (SMT-COMP 2018) and ranked first in the category QF LIA (quantifier-free linear integer arithmetic) [START_REF]SMT-COMP 2018 results for QF LIA (main track)[END_REF] and second in the category QF LRA (quantifier-free linear rational arithmetic) [START_REF]SMT-COMP 2018 results for QF LRA (main track)[END_REF]. This system description explains the main features that led to the success of SPASS-SATT. We do not only describe the relevant techniques, but also show their specific impact on dedicated groups of examples from the SMT-LIB by experiments.

By far not all techniques presented in this system description are unique features of SPASS-SATT. The techniques that appeared first in SPASS-SATT are the unit cube test and bounding transformations explained in Section 2. Concerning preprocessing, SPASS-SATT is the first SMT solver implementing the small-clause-normal-form algorithm, see Section 4. Further important techniques implemented in SPASS-SATT have already been available in other SMT solvers such as CVC4 [START_REF] Barrett | CVC4. In CAV[END_REF], MathSAT [START_REF] Cimatti | The MathSAT5 SMT Solver[END_REF], Yices [START_REF] Dutertre | Yices 2.2[END_REF], and Z3 [START_REF] De Moura | Z3: an efficient SMT solver[END_REF], but not all in one tool: (i) the implementation of branch and bound as a separate theory solver and a number of improvements to the simplex implementation such as a priority queue for pivot selection, integer coefficients instead of rational coefficients, dynamically switching between native and arbitrary precision integers, and backing-up versus recalculating simplex states, all in Section 2, (ii) decison recommendations, unate propagations, and bound refinements for the interaction between the SAT and theory solver, in Section 3, and (iii) preprocessing techniques for if-then-else operators and pseudo-boolean inequalities, in Section 4. Although these techniques are contained in existing SMT solvers, not all have been described in the respective literature. The paper ends with a discussion of future extensions to SPASS-SATT in Section 5.

The benchmark experiments with SPASS-SATT consider the 6947 SMT-LIB benchmarks for quantifier-free linear integer arithmetic (QF LIA) [START_REF] Barrett | The Satisfiability Modulo Theories Library (SMT-LIB)[END_REF]. For the experiments, we used a Debian Linux cluster and gave SPASS-SATT for each problem one core of an Intel Xeon E5620 (2.4 GHz) processor, 8 GB RAM, and 40 minutes. The results are depicted as scatter plots and in each of them we compare the default configuration (i.e., without any command line options) of SPASS-SATT (horizontal axis) with an alternative configuration of SPASS-SATT (vertical axis). (The SMT-COMP results were obtained with the default configuration; by default all presented techniques are turned on.)

2 SPASS-IQ: An LA Theory Solver SPASS-SATT's theory solver, called SPASS-IQ, decides conjunctions of linear arithmetic inequations. It is divided into two main components: a simplex implementation for handling linear rational arithmetic and a branch-and-bound implementation for handling linear mixed/integer arithmetic.

However, the division between the two components is in all truth not that strict. The branch-and-bound implementation is more of a supervisor for the simplex implementation. To be more precise, the branch-and-bound implementation coordinates the search for a mixed or integer solution, but the majority of the actual search/calculation is still done by the simplex implementation. For most QF LIA benchmark instances (4894 out of 6947 instances), this supervision is not even necessary; i.e., SPASS-SATT solves these instances with just the simplex implementation as its theory solver. This means that SPASS-SATT's efficiency on the QF LIA benchmarks also highly depends on the efficiency of our simplex implementation and not just on the extensions and optimizations to our branch-and-bound implementation.

The simplex implementation inside SPASS-IQ is based on a specific version [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF] of the dual simplex algorithm [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]. The overall efficiency of our simplex implementation is heavily influenced by the efficiency of the data structures that we use. Our most important data structure features are: 1) Priority Queue for Pivot Selection: Instead of iterating over all basic variables when searching for violated basic variables, we collect the basic variables in a priority queue as soon as they become violated.

2) Integer Coefficients Instead of Rational Coefficients: We avoid rational coefficients in our simplex tableau by multiplying each equation in the tableau with the common denominator of the equations coefficients. As a result, each basic variable also has a coefficient, but all coefficients are integers. This transformation roughly halves the cost of most tableau operations because we do not need to consider rationals which are typically represented by two integers (the numerator and the denominator).

3) Dynamically Switching between Native and Arbitrary-Precision Integers: We use the arbitrary-precision arithmetic library FLINT to represent our integers [START_REF] Hart | FLINT: Fast Library for Number Theory[END_REF]. It dynamically switches between native C integer and arbitrary-precision types. 4) Backup vs. Recalculation: In contrast to Dutertre and de Moura's version of the simplex algorithm, our simplex backtrack function recalculates a satisfiable assignment instead of loading a backup of the last satisfiable assignment. SPASS-IQ's second set of decision procedures revolves around an implementation of the branch-and-bound (BnB) algorithm [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]. Most SMT solvers implement branch and bound through a technique called splitting-on-demand [START_REF] Barrett | Splitting on demand in SAT modulo theories[END_REF], which delegates some of the branch-and-bound reasoning to the SAT solver. In order to keep more control over the branch-and-bound reasoning, we decided against splitting-on-demand and implemented branch and bound as a theory solver separate from the SAT solver. This also made it easier to complement branch and bound with other decision procedures:

The first two extensions that we discuss here are simple rounding (turn off with -LASR 0) and bound propagation (turn off with -LABP 0) [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF], which are both classical additions to most branch-and-bound implementations. For simple rounding, we round any rational solution computed during the branch-andbound search to the closest integer assignment and check whether this is already an integer solution. For bound propagation, we propagate new bounds from ex-isting bounds at every node in our branching tree. Although both techniques are very popular, we could only measure a minor impact on SPASS-SATT's performance on the QF LIA benchmarks. With simple rounding we solve only one instance faster and with bound propagations we solve only 10 additional instances. In part, this is due to our next two extensions that make simple rounding and bound propagation in many cases unnecessary.

The next extension we discuss is the unit cube test (turn off with -C 0). It determines in polynomial time whether a polyhedron, i.e., the geometric representation of a system of inequalities, contains a hypercube parallel to the coordinate axes with edge length one [START_REF] Bromberger | Fast cube tests for LIA constraint solving[END_REF][START_REF] Bromberger | New Techniques for Linear Arithmetic: Cubes and Equalities[END_REF]. The existence of such a hypercube guarantees a mixed/integer solution for the system of inequalities.

The unit cube test is only a sufficient and not a necessary test for the existence of a solution. There is at least one class of inequality systems, viz., absolutely unbounded inequality systems [START_REF] Bromberger | Fast cube tests for LIA constraint solving[END_REF][START_REF] Bromberger | New Techniques for Linear Arithmetic: Cubes and Equalities[END_REF], where the unit cube test is also a necessary test and which are much harder for many complete decision procedures.

The plot in Figure 1(a) shows that SPASS-SATT employing the unit cube test solves 56 additional benchmark instances from the QF LIA benchmarks and solves 705 instances more than twice as fast. 3 Moreover, the unit cube test causes only a minor overhead on problems where it is not successfully applicable.

The final extension that we discuss are bounding transformations (turn off with -B 0). Branch and bound alone is an incomplete decision procedure and only guarantees termination on bounded problems, i.e., problems where all variables have an upper and a lower bound. For this reason, we developed two transformations that reduce any unbounded problem into an equisatisfiable problem that is bounded [START_REF] Bromberger | A reduction from unbounded linear mixed arithmetic problems into bounded problems[END_REF]. The transformed problem can then be solved with our branch-and-bound implementation because it is complete for bounded problems.

The plot in Figure 1(b) shows that SPASS-SATT employing the bounding transformation solves 169 additional benchmark instances from the QF LIA benchmarks and solves 167 instances more than twice as fast. 4 Moreover, the bounding transformation causes only a minor, almost immeasurable overhead on problem instances where it is not successfully applicable.

CDCL(LA): SAT and Theory Solver Interaction

SPASS-SATT uses at its core a CDCL(LA) implementation that combines our CDCL (conflict-driven-clause-learning)-based SAT solver SPASS-SAT with our LA theory solver SPASS-IQ. The result is a decision procedure for ground lineararithmetic formulae in clause normal form. In this section, we quickly explain how our theory solver and SAT solver interact. To this end, we list in Figure 3 the main interface functions of our SAT solver and theory solver and show through RatCheck() GetPhase(L): selects the phase P for the decision literal L. LearnClauses(): Adds a set S of clauses to N that correspond to unate propagations and bound refinements; returns ∅ if there are no clauses to learn. RatCheck(): determines a rational solution for the asserted literals; returns if a rational solution exists; otherwise, returns ⊥. IntCheck(): determines an integer solution for the asserted literals; ends CDCL(LA) and returns sat if an integer solution exists; otherwise, returns ⊥. GetConflict(): returns a clause C that explains the theory conflict. LABacktrack(D): removes all asserted literals that were added after decision level D; recalculates a rational solution for the remaining asserted literals. Fig. 2. CDCL(LA) as implemented in SPASS-SATT a flow graph how they interact. The main focus of this section, however, is to explain in which way our implementation of CDCL(LA) differs from more general frameworks for CDCL(T), also called DPLL(T) [START_REF] Barrett | Splitting on demand in SAT modulo theories[END_REF][START_REF] Ganzinger | DPLL(T): fast decision procedures[END_REF][START_REF] Nieuwenhuis | DPLL(T) with exhaustive theory propagation and its application to difference logic[END_REF][START_REF] Nieuwenhuis | Solving SAT and SAT modulo theories: From an abstract Davis-Putnam-Logemann-Loveland procedure to DPLL(T )[END_REF].

Assert(L) IntCheck() IncLvl () LABacktrack(D) GetConflict() LearnClauses() GetPhase(L') Decide(L',P) UnitPropagate(N) AnalyzeConflict(C) SelectDecLit() SATBacktrack(C') Assert(L) L ⊤ ⊥ ⊥ ⊥ ⊤ ∅ ⊤ ⊥ P S L' ∅ ∅ C C' D SAT
There are three key points that we have changed compared to the more general frameworks for CDCL(T). First of all, we rely on "weakened early pruning" [START_REF] Sebastiani | Lazy satisability modulo theories[END_REF], i.e., we only use a weaker but faster check to determine theory satisfiability for partial (propositionally abstracted) models. We do so because IntCheck(), i.e., checking for an integer solution, is too expensive and not incrementally efficient enough to be checked more than once per complete (propositionally abstracted) model. As a compromise, we at least check the partial model with RatCheck(), i.e., we check for a rational solution, before we add a(nother) decision literal to the model with Decide(L,P ).

As our second key change, we let the theory solver select via GetPhase(L) the phase of the next decision literal L, i.e., whether Decide(L,P ) will add the positive or the negated version of L to the model. We call this technique a decision recommendation (turn off with -p 0). 5 Finally, we use theory reasoning via the function LearnClauses() to find and learn new clauses implied by the input formula. The reasoning techniques we use for this purpose are unate propagations and bound refinements as proposed in [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF].

In Figure 3(a), we examine the impact that decision recommendations have on SPASS-SATT's performance on the QF LIA benchmarks. With decision recommendations it can solve 129 additional problems. Moreover, it becomes more than twice as fast on 389 problems, but only twice as slow on 58 problems. The benchmark family that is impacted the most by decision recommendations is convert with 116 additionally solved problems (see Figure 3(b)). Although SPASS-SATT frequently and regularly performs unate propagations on the QF LIA benchmarks, we are unable to observe any consistent benefit or drawback from this interaction technique. The impact of bound refinements is also relatively minor and SPASS-SATT solves only 24 additional problem instances if they are activated. 4 Preprocessing

Many real world applications can be encoded as linear integer arithmetic formulas, and some of those applications are too specialized to be efficiently handled by our rather general CDCL(LA) implementation. To resolve this, we have complemented SPASS-SATT with several specialized preprocessing techniques. 6Complex input formulas are typically transformed into CNF by a Tseitinstyle renaming [START_REF] Tseitin | On the complexity of derivations in propositional calculus[END_REF] using a static criterion which subformula to replace by a fresh propositional variable. SPASS-SATT includes the small clause normal form algorithm [START_REF] Nonnengart | Computing small clause normal forms[END_REF]. Instead of a static criterion, the number of clauses with or without a renaming is compared and a fresh propositional variable is introduced only if a renaming eventually yields fewer clauses. This results in a more compact CNF with strictly fewer additional propositional variables. To this end we extended the small clause normal form algorithm to ITE formulas. For an ITE formula (ite t 1 t 2 t 3 ), simplified by the below techniques and potentially contained in some formula f , we compare the number of clauses generated out of replacing the formula with a fresh variable P in f and adding P ↔ [(t 1 → t 2 )∧(¬t 1 → t 3 )] with a direct replacement of the ITE formula by the before mentioned conjunction of two implications. This test can be carried out in constant time after having once filled respective data structures. The set up of the data structures needs one run on the overall formula, i.e., can be computed in linear time [START_REF] Nonnengart | Computing small clause normal forms[END_REF].

SPASS-SATT also has five specialized preprocessing techniques for if-thenelse expressions (ITE). Our first technique, if-then-else reconstruction (turn off with -PPIR 0), rebuilds if-then-else operations that were already preprocessedaway by the creators of the input problem. The reconstruction then allows us to apply simplifications missed during the creation of the input problem. To this end, we check whether the first conjunctive layer of our formula f := (and . . . t i . . . t i . . .) contains any pair of clauses t i , t i that match the clauses added by the standard if-then-else elimination, i.e., t i ∈ T i and t i ∈ T i \ {t i }, where T i = {(or t i1 (= y i t i3 )), (or (not t i1 ) (= y i t i2 )), (or (= y i t i2 ) (= y i t i3 ))} and y i is an arithmetic variable. If we find such pairs t i , t i , then we remove them from f and replace all remaining occurrences of y i in f with (ite t i1 t i2 t i3 ).

The next three techniques are all dedicated to so-called constant if-then-else expressions (CITEs). A CITE is either a leaf, i.e., an arithmetic expression that can be simplified to a number a ij ∈ Q, or a branch, i.e., an if-then-else expression (ite t 1 t 2 t 3 ) where t 2 and t 3 are again CITEs.

The first CITE technique is called shared monomial lifting (turn off with -PPIL 0) and we use it to increase the number of CITEs in our formula. It traverses the subterms in our formula in bottom-up order and transforms all subterms t := (ite t 1 (+ q q ) (+ q q)) into (+ q (* 1 (ite t 1 (+ q ) (+ q)))). (We assume here for simplicity that the shared part q appears after the unshared parts q and q. In reality SPASS-SATT has to find and extract the shared parts.)

The second technique is called CITE simplification (turn off with -PPIS 0) and it simplifies atoms (o t 1 t 2 ), where t 1 and t 2 are CITEs and o is one of the operators <=, =, >=. To be more precise, the technique essentially pushes the comparison operator o recursively down the CITE branches and greedily simplifies any branch to true or false if possible. For more details, see Section 4 in the paper by Kim et al. on efficient term-ITE conversion [START_REF] Kim | Efficient term-ite conversion for satisfiability modulo theories[END_REF].

The third technique is called CITE bounding (turn off with -PPIB 0) and eliminates the remaining CITEs. Thanks to it, we only introduce one new variable for each topmost CITE instead of one variable for each CITE branch. Moreover, this new variable describes only a small set of values (often equivalent to just the CITE leaves) because it is bounded as tightly as possible.

As its first step, CITE bounding creates one new integer variable x j for each topmost CITE expression t j and replaces the occurrences of t j in our formula f with x j . As its second step, CITE bounding extends f to the formula f := (and f f 1 f 1 . . . f m f m ), where (i) f j is equivalent to t j except that all leaves a ij of t j are replaced by the equations (= x j a ij ), (ii) f j := (and (>= x j a minj ) (<= x j a maxj )), and (iii) a minj is the smallest leaf in t j and a maxj is the largest leaf in t j . As its last step, CITE bounding replaces all occurrences of x j in f with ( * a gj x j ), where a gj is the greatest common divisor of the leaves a ij in t j .

The final ITE technique handles nested conjunctive if-then-else expressions (AND-ITEs). An AND-ITE is a series of nested if-then-else expressions that can be simplified to a conjunction. For instance, if t i := (ite t i t i+1 false) for i = 1, . . . , n, then t 1 is an AND-ITE equivalent to (and t 1 . . . t n t n+1 ). Naturally, we transform these AND-ITEs into actual conjunctions and we call this process compression. However, we compress an AND-ITE t 1 only if all of its actual AND-ITE subterms t i appear only inside t 1 and only once. If this is the case, then we first replace all occurrences of t 1 in f by a new propositional variable p j and extend our formula f to the formula (and f (= p j t 1 )), where t 1 is the compressed form of t 1 . (If this is not the case, then we simply replace and compress the AND-ITE subterms in t i first.) We do the compression in this way to strengthen the connection of the AND-ITEs that have multiple occurrences in f . The above described technique is called if-then-else compression (turn off with -PPIC 0) and it was first presented by Burch in [START_REF] Burch | Techniques for verifying superscalar microprocessors[END_REF]. However, Burch used it to simplify control circuits and not SMT input problems.

Last but not least, SPASS-SATT also provides a preprocessing technique for pseudo-boolean problems [START_REF] Hammer | Boolean methods in operations research and related areas[END_REF], i.e., linear arithmetic problems where all integer variables x j have bounds 0 ≤ x j ≤ 1, which we call pseudo-boolean variables. To be more precise, SPASS-SATT recognizes clauses that are encoded as linear pseudo-boolean inequalities (i.e., inequalities containing just pseudo-boolean variables) and turns them into actual clauses (turn off with -PPCC 0). (This technique goes back to the NP-hardness proof of 0-1 programming [START_REF] Karp | Reducibility among combinatorial problems[END_REF].) However, SPASS-SATT only transforms inequalities containing at most three variables because it would otherwise fail to solve some of the problems from the pidgeons benchmark family.

The convert benchmark family contains problems that are relatively hard unless SPASS-SATT uses the right combination of techniques. To be more precise, SPASS-SATT solves only two-thirds of the 319 instances if one of the following three techniques is missing: (i) recalculating simplex states during backtracking (Section 2), (ii) decision recommendations (Section 3), or (iii) the small CNF transformation. SPASS-SATT with all three techniques solves all instances in less than 2 seconds (see Figure 3(b)).

The nec smt benchmark family contains problems with many nested if-thenelse and let expressions. SPASS-SATT can handle most of them if we first apply our constant if-then-else simplifications and our conjunctive if-then-else compression. In Figure 4(b), we see that SPASS-SATT without our preprocessing techniques solves only 1422 out of the 2800 benchmark instances and is by far slower on the instances it can solve. SPASS-SATT with our preprocessing techniques solves 2782 out of the 2800 benchmark instances.

The rings benchmark family encodes associative properties on modular arithmetic with the help of if-then-else expressions. With a combination of shared monomial lifting and constant if-then-else bounding, these problems become almost trivial to solve. In fact, SPASS-SATT needs less than one second for each problem instance and needs only techniques for linear rational arithmetic to solve each of them (Figure 4(c)).

The rings preprocessed benchmark family is equivalent to the rings benchmark family except that all if-then-else-operations were eliminated by standard if-then-else elimination [START_REF] Griggio | A practical approach to satisfiability modulo linear integer arithmetic[END_REF]. We can use the same trick as for the rings benchmark family if we first use our if-then-else reconstruction technique that reverses the standard if-then-else elimination.

The pb2010 benchmark family is a set of industrial problems taken from the pseudo-boolean competition 2010. With its pseudo-boolean preprocessing, SPASS-SATT solves 22 additional benchmark instances from the 81 instances in the pb2010 benchmark family (see Figure 4(e)). Moreover, SPASS-SATT solves all of these instances without its branch-and-bound implementation.

Conclusion and Future Work

We have presented SPASS-SATT our complete solver for ground linear arithmetic and have explained which techniques make it so efficient in practice. To summarize, SPASS-SATT is so efficient because (i) we have optimized the data structures in our simplex implementation, (ii) we have combined branch and bound with the unit cube test and the bounding transformation, (iii) we have added decision recommendations to our CDCL(LA) framework, (iv) we have added a small CNF transformation, and (v) we have added specialized preprocessing techniques for if-then-else expressions and pseudo-boolean inequalities.

Almost all of the presented techniques can be applied incrementally, however this is not always useful. For the partial models computed by the SAT solver, we only apply the simplex method, unate propagation, and bound refinements incrementally. If we ever extend SPASS-SATT to handle theory combinations or incremental SMT-LIB problems, then we would also apply branch-and-bound, unit cubes, and bounding transformations incrementally; but only on the models generated during Nelson-Oppen combination or between two (check-sat) calls.

For future research, we plan to extend SPASS-SATT to quantified linear arithmetic. Moreover, we plan to complement SPASS-SATT with several specialized decision procedures. For instance, SPASS-SATT could handle (almost) pseudo-boolean problems (e.g., benchmark families pb2010, miplib2003) much more efficiently if we extended branch and bound with a SAT based arithmetic decision procedure [START_REF] Bromberger | Linear integer arithmetic revisited[END_REF][START_REF] Jovanović | Cutting to the chase solving linear integer arithmetic[END_REF][START_REF] Nieuwenhuis | The IntSat method for integer linear programming[END_REF].

Fig. 1 .

 1 Fig.1(a)

Fig. 3 .

 3 Fig.3(a)

3 FigFig. 4 .

 34 Fig.4(a)

  Assert(L): returns ⊥ if the literal L contradicts another asserted literal. IncLvl(): notifies the theory solver that a new decision level was reached.

			Decide(L,P ): if P is , then L is added
	Solver CDCL(LA) Input: N unsat N	Theory Solver CDCL(LA) Output: C or sat unsat sat	to the model; otherwise, ¬L is added. SelectDecLit(): returns ∅ if the model satisfies all clauses; otherwise, returns a literal L that is undefined in the model. UnitPropagate(N ): returns ∅ if no lit-eral can be unit propagated; otherwise, se-lects a literal L that can be unit propa-gated and adds it to the model; if a clause C ∈ N evaluates to ⊥ under the new model, then returns C; otherwise, returns L. AnalyzeConflict(C): derives a clause C which is the negation of the literals that

led to the conflict in C; ends CDCL(LA) and returns unsat if C is empty; otherwise, returns C . SATBacktrack(C ): adds C to the clause set N and backtracks to the maximum decision level D where C is still satisfiable; retuns D.

These instances belong to the dillig[START_REF] Dillig | Cuts from proofs: A complete and practical technique for solving linear inequalities over integers[END_REF], CAV-2009[START_REF] Dillig | Cuts from proofs: A complete and practical technique for solving linear inequalities over integers[END_REF], slacks[START_REF] Jovanović | Cutting to the chase solving linear integer arithmetic[END_REF], 20180326-Bromberger[START_REF] Bromberger | A reduction from unbounded linear mixed arithmetic problems into bounded problems[END_REF], and prime-cone benchmark families[START_REF] Jovanović | Cutting to the chase solving linear integer arithmetic[END_REF], which together contain more than 1483 instances of absolutely unbounded problems.

These instances belong to the 20180326-Bromberger[START_REF] Bromberger | A reduction from unbounded linear mixed arithmetic problems into bounded problems[END_REF], arctic-matrix[START_REF] Codish | Exotic semi-ring constraints[END_REF], cut lemmas[START_REF] Griggio | A practical approach to satisfiability modulo linear integer arithmetic[END_REF], slacks[START_REF] Jovanović | Cutting to the chase solving linear integer arithmetic[END_REF], and tropical-matrix[START_REF] Codish | Exotic semi-ring constraints[END_REF] benchmark families.

In our own theory solver, GetPhase(L)= if the current theory assignment satisfies L or if Assert(¬L) would return ⊥. Otherwise, GetPhase(L)=⊥. (If unate propagations are enabled, then the case that Assert(¬L) would return ⊥ is impossible because L would have been unate propagated.)

All preprocessing techniques are also contained in CVC4[START_REF] Barrett | CVC4. In CAV[END_REF] with the exception of the small CNF. The implementation is so efficient because we employ a shared term representation and cache all intermediate results.


Elisabeth Longatte
email: elisabeth.lacazedieu@edf.fr

Erwan Liberge

Marie Pomarede

Jean-François Sigrist

Aziz Hamdouni

Parametric study of flow-induced vibrations in cylinder arrays under single-phase fluid cross flows using POD-ROM

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

Introduction

Optimization of safety barrier reliability and Uncertainty Quantification (UQ) of physical models give rise to long-term research programs involving designers and engineers working on maintaining systems under operating conditions. The present article focuses on vibration risk assessment in heat exchangers in the context of lifecycle control and increase of mechanical components in spite of the very constrained conditions they are submitted to. In Pressurized Water Reactors (PWR) steam generators ensure the transition between the primary and secondary loops and are used to convert water into steam from heat produced by the core made of fuel assemblies. Each heat exchanger can measure up to 20 m and is made of regular confined arrangements of several thousands of thin elongated cylinders whose section diameter may be less than 2 cm. These cylinders are conveying the primary fluid and steam is produced, coming from the water on the shell side, so-called the secondary fluid. The heat exchanges take place at this stage in this area. Then the secondary steam is delivered to turbines for electric power generation. A brief dimensionless analysis of this geometrically complex system shows that as a first approximation under current conditions the dynamical behavior of these vibrating cylinders is independent on heat transfers. Therefore it is not necessary to account for thermal effects in the context of dynamical analysis. Due to geometry and to thermohydraulics conditions, most effects responsible for vibrations of cylinder arrays are due to Fluid Structure Interaction (FSI) and Flow-Induced Vibration (FIV) occurring at the external shell of the cylinders. They are due to the hydrodynamic unsteady load exerted by the external fluid on these structures. Moreover the component is such that cylinders are submitted to two-phase flows at the top and to single-phase flows at the bottom. Therefore as a first step, in the present article, one focuses on single-phase FIV at the bottom of the heat exchanger where flows are crossing mostly transversally the array. The purpose is to address the multi-physics multi-scale issue of stability analysis of dynamical response of cylinder arrangements submitted to external single-phase fluid flows in order to be able to establish accurate stability maps applicable for design. Therefore parameter perturbations must be investigated in order to account for possible large variations of parameter values during the lifecycle of these mechanical components.

As far as numerical modeling is concerned, a fully three-dimensional microscopic-scale model of the bottom area of the heat exchanger tube array under single-phase fluid cross flow under real operating conditions characterized by Reynolds numbers of order of 10 5 would lead to solving systems involving more than 10 10 degrees of freedom which is still unreachable in practice even by using most advanced High Performance Computing (HPC) resources. Therefore superposition methods have been introduced in the context of first order asymptotic developments in order to perform linear stability analysis of these systems and exhibit major mechanisms involved in dynamical instability occurrence. Using space domain decompositions dimensionless parameter effects have been studied independently. Nevertheless HPC simulations performed on small-size domains are not yet convenient for design since they still require long computations of several hours to several days according to Reynolds numbers to be considered and to turbulence models to be involved. Therefore Reduced-Order Modeling (ROM) is investigated in order to deal with this problematic of reducing computational time. It is investigated in the present work by considering as a first step a two-parameter configuration where only Reynolds number and reduced velocity can vary while other parameters like mass ratio, Scruton number and pitch ratio of the cylinder arrangement are given. The basic idea is to use few High-Fidelity (HF) computations with a classical method (such as Finite Volumes and Finite Elements) to build low computational time models. Then ROMs have an initial cost, but this cost can be charged off if the models are used for large set of parametric values. The success of such approaches can be quantified by the large field of their applications. For example, ROMs have been used for studying different configurations of thermal management of data centers [START_REF] Samadiani | Reduced order thermal modeling of data centers via proper orthogonal decomposition: A review[END_REF], for controlling wave energy converters [START_REF] Hesam | An Approach to Reduced-Order Modeling and Feedback Control For Wave Energy Converters[END_REF], fluid flow control at high Reynolds numbers [START_REF] Semaan | Reduced-order modelling of the flow around a high-lift configuration with unsteady Coanda blowing[END_REF][START_REF] Balajewicz | Minimal subspace rotation on the Stiefel manifold for stabilization and enhancement of projection-based reduced order models for the compressible Navier-Stokes equations[END_REF] or in aeronautics [START_REF] Kim | Parametric model reduction for aeroelastic systems: Invariant aeroelastic modes[END_REF]. One way also explored is computation on smartphone, for which ROM approaches seem to be the appropriate solution. For example, [START_REF] Modesto | Proper generalized decomposition for parameterized Helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation[END_REF] proposed an application for smartphone for computing disturbances in the Barcelona harbor in quasi real time using a ROM approach. On the same thing, ROM using a few/small resources could be used in a context of volunteer computing [START_REF] Nouman Durrani | Volunteer computing: requirements, challenges, and solutions[END_REF] for a large parametric studies.

Few works exists on ROM development for FIV or FSI context. All works consist in adapting or extending the Proper Orthogonal Decomposition (POD) to this problematic. The POD is the most famous method for ROM in fluid mechanics and it has naturally been extended and adapted to the cases of flow in the presence of moving boundaries. For FSI cases similar to aeroelasticity, i.e where the fluid domain could be considered as fixed, the Navier-Stokes equations are linearized and projected on a POD basis (Barone et al., 2009b;[START_REF] Lieu | Reduced-order fluid/structure modeling of a complete aircraft configuration[END_REF]. For small vibrations at the surface of the structure, [START_REF] Bourguet | Reduced-order modeling of transonic flows around an airfoil submitted to small deformations[END_REF] proposed a Hadamar formulation associated with a ROM-POD. For large fluid domain motion of deformation, [START_REF] Rozza | Reduced basis methods for stokes equations in domains with non-affine parameter dependence[END_REF] used POD on a parametric mesh and next built ROM for stationary problems on different meshes. In case of imposed displacement of the structure, [START_REF] Balajewicz | Reduction of nonlinear embedded boundary models for problems with evolving interfaces[END_REF] proposed an immersed boundary ROM with a modified formulation. The ROM is directly built in an immersed boundary solver and is available only for imposed displacement of the solid domain. Among the previous proposed methods, none allows to solve a FSI problems with a large displacement of the structure, i.e cases where the fluid modifies the behavior of the structure, which interactively also changes the flow. The issue of using POD for moving domains deals with the paradox of computing a POD basis which is a spatial basis on a changing (i.e. a time-evolving) domain. Liberge and Hamdouni (2010) and Liberge et al. (2010) proposed a POD-multiphase formulation available for FSI problems with large displacement of the domain. The method consists in computing a POD basis of an interpolated velocity field on a reference grid, and next in projecting a multiphase formulation on this POD basis. It is recalled in the next section. This method has been successfully used for small and large displacements of a cylinder in a cross fluid flow for known parameters.

This paper explores the behavior of POD-multiphase ROM when the parameter values are different from those used to build the POD basis. In the first part, the principles of the POD, the POD-multiphase and the ROM parameter sensitivity analysis methods are explained. Next, the parameter sensitivity of the cylinder vibrations in a cross fluid flow is studied and results are discussed.

Theory

One considers Ω ⊂ R 3 such that Ω = Ω f (t) ∪ Ω s (t) ∪ Γ i (t) with Ω s the solid domain, Ω f the fluid domain and Γ i the fluid-solid interface. One defines the external boundary of the domain : Γ f = ∂Ω \ Γ i . The normal vector n is oriented as pointing out towards outside the solid domain. One defines a velocity field u over the whole domain Ω × [0, T ] where T corresponds to the last time where the dynamics of the system is considered.

The Proper Orthogonal Decomposition (POD)

The Proper Orthogonal Decomposition method being a well-known method, just a practical approach introduced by [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] and the principal properties are exposed in the following. For more informations, one can refer to [START_REF] Allery | Experimental and numerical POD study of the Coanda effect used to reduce self-sustained tones[END_REF], Liberge and Hamdouni (2010) and Liberge et al. (2010).

Considering a set of M snapshots u (t k ) , k = 1 • • • M of the velocity field, the POD method consists in solving the following eigenvalue problem :

CA k = λ k A k (1) with C ij = ∫ Ω u (x, t i ) u ( x, t j
) dx is a time correlation tensor of the snapshots of the velocity field. Next, the vectors of the spatial basis are obtained with the decomposition :

Φ k (x) = M ∑ m=1 A k m u (x, t m ) (x) (2)
The POD basis has the following properties :

• There is a set of positive eigenvalues (λ i ) i≥1 which decreases to 0,

λ 1 ≥ λ 2 ≥ • • • ≥ λ i ≥ • • • and λ i → 0 (3)
and a set of vectors (Φ i ) i≥1 which is a Hilbertian basis. • The vectors (Φ i ) are orthogonal and can be normalized:

( Φ i , Φ j ) = ∫ Ω Φ i (x) • Φ j (x) dx = δ ij (4)
In case of an incompressible flow, the POD basis fulfills the free divergence, i.e. divΦ i = 0, i ≥ 1. • The values of the temporal coefficients a i at specific times t k , k = 1 • • • M are obtained from the projection of u onto the basis :

a i (t k ) = (u (t k ) , Φ i ) (5)
This procedure is called POD direct projection.

• The eigenvalue λ i is the energy captured by the vector Φ i . For a given N, the POD decomposition is the best energy decomposition that can be obtained.

Usually, in fluid mechanics, the POD is computed for the fluctuating velocity field u ′ and the POD basis is

( Φ ′ i ) .

The POD-multiphase formulation for FSI problems

The classical POD approach leads to a spatial basis, and needs fixed domains. Indeed, if the domain Ω is time dependent, the eigenvalue problem (1) cannot be computed since the snapshots do not have the same domain definition at times t i and t j . To cancel this drawback, Liberge and Hamdouni (2010) and Liberge et al. (2010) have proposed to use a global domain Ω with a fixed grid which recovers all the time variant domains

Ω = Ω f (t) ∪ Ω s (t) ∪ Γ i (t) as displayed in Fig. 2.
The method consists in computing snapshots using a classical Fluid Structure Interaction approach, which can use time variant grid. Next, the snapshots are interpolated from the time variant grid to the reference one. The snapshots can be computed by a classical solver using Arbitrary Lagrange Euler (ALE) method involving moving grid. Then a global velocity field is built on Ω as follows:

u(x, t) = u s (x, t)1 Ωs (x, t) + u f (x, t)1 Ω f (x, t) (6)
with u s (x, t) the solid velocity field defined on space Ω s (t) and u f (x, t) the fluid velocity field defined on space Ω f (t) at time t. Characteristic functions are respectively expressed by 1 Ωs (x, t) and 1

Ω f (x, t) with 1 Ωs (x, t) = 1 -1 Ω f (x, t). Each function 1 Ω i (for i = f or s) is defined by : 1 Ω i (x, t) = { 1 if x ∈ Ω i (t) 0 else
Each variable (density and viscosity) and each field is decomposed following the example of Eq. ( 6).

Examples of time-dependent space distributions of characteristic functions are displayed in Fig. 3 for several grid refinement.

Reduced-order dynamical model

Applying the POD to the velocity field (6) leads to a spatial basis defined on Ω.

We focus in this paper on cases where solid domain can be modeled using a rigid body coupled with springs and dampers. The non fluid constraints are expressed as external forces F f and momenta M s .

Then, the ROM is built using the method proposed by Liberge and Hamdouni (2010). The principle is developed in Liberge and Hamdouni (2010) and it is recalled in the present paper.

The formulation consists in extending the Navier-Stokes equations to the solid domain by adding a rigidity constraint Eq. ( 7) defined on the solid domain and a Lagrange multiplier Λ associated with the constraint as follows:

D(u) = 0 in Ω s (t) (7)
where D (•) is defined by :

D (•) = 1 2 ( ∇ • +∇ T • ) .
This leads to the following variational formulation on the global domain Ω :

H v Γ = { u|u ∈ H 1 (Ω) , u = u Γ (t) on ∂Ω \ Γ I } , H 0 = { u|u ∈ H 1 (Ω) , u = 0 on ∂Ω \ Γ I } , L 2 0 (Ω) = { q ∈ L 2 (Ω) | ∫ Ω qdx = 0 } , ∀u * ∈ H 0 and q ∈ L 2 (Ω) , find u ∈ H v Γ , p ∈ L 2 0 (Ω) , Λ ∈ H 1 (Ω s (t)) such that ∫ Ω ρ ( ∂u ∂t + u∇u ) u * dx - ∫ Ω p∇ • u * dx + ∫ Ω q∇ • vdx + ∫ Ω 2µD (u) : D ( u * ) dx + ∫ Ωs(t) D (Λ) : D ( u * ) dx = ∫ Ω fu * dx (8)
where ρ and µ are defined on the global domain Ω :

ρ = ρ f ( 1 -1 Ωs ) + ρ s 1 Ωs , µ = µ f ( 1 -1 Ωs ) + µ s 1 Ωs . ( 9
)
µ s is the penalization factor associated with the Lagrange multiplier. f denotes the global volumique forces defined as

f = f f ( 1 -1 Ωs ) + f s 1 Ωs (10)
where f f designates the standard external volumic hydrodynamic forces exerted by the fluid and f s is the other part of external volumic forces acting on the rigid body different from fluid forces. f s satisfies the two following equations :

F s = ∫ Ωs f s dx and M s = ∫ Ωs f s ∧ xdx (11)
Next, the POD-ROM is built by choosing POD modes Φ i , i = 1 • • • N for the virtual velocity field and the decomposition of u on the POD basis.

For each n = 1, . . . , N:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N ∑ i=1 da i dt A in + N ∑ i=1 N ∑ j=1 a i a j B ijn + N ∑ i=1 a i C in + D n + E n = 0 D (u) = 0 in Ω s (t) ∂1 Ωs ∂t + u • ∇1 Ωs = 0 (12)
where :

A in = ∫ Ω ρΦ i • Φ n dx D n = ∫ Ω 1 Ωs Tr [ D(Λ)D(Φ n ) ] dx B ijn = ∫ Ω ρ(Φ i • ∇)Φ j • Φ n dx E n = ∫ Ω fΦ n dx C in = 2 ∫ Ω µTr [ D(Φ i )D(Φ n ) ] dx
As explained in Liberge and Hamdouni (2010) the formulation (12) does not need to compute explicitly the fluid forces on the structure because they are intrinsic to the formulation. In addition, the equivalent formulation using the fluctuating velocity field is described in Appendix A. Moreover the method can be extended to the case of deformable structure as described in Appendix B.

Remark. A POD-Galerkin method is involved in the present work. It is well-known that it can lead to instabilities, that is why the use of additional viscosities associated to each POD mode [START_REF] Cazemier | Proper orthogonal decomposition and law dimensional models for driven cavity[END_REF] or Least-Square approach [START_REF] Tallet | A minimum residual projection to build coupled velocity-pressure POD-ROM for incompressible Navier-Stokes equations[END_REF][START_REF] Carlberg | Galerkin v. Least-squares Petrov-Galerkin projection in nonlinear model reduction[END_REF] is often preferred in the context of fluid mechanics applications. However this is not necessary in the present study since as explained in Liberge and Hamdouni (2010), the solid viscosity used to penalize the rigid body plays the role of the viscosity of [START_REF] Cazemier | Proper orthogonal decomposition and law dimensional models for driven cavity[END_REF] and avoids instabilities.

Application to FIV of cylinder array

Configuration

In what follows, the model reduction method is applied to a configuration involving a periodic cylinder array made of one central moving cylinder and 8 truncated non-moving cylinders submitted to an external cross flow. The problem is two-dimensional thanks to the choice of parameters since the Reynolds number does not exceed 1000. Reduced velocity varies between 1 and 4. Other parameters are fixed. Pitch ratio is 1.44, Scruton number is taken to 0 and mass ratio equals 7.8. (See Fig. 4).

These parameters are defined as follows: pitch ratio P r = P/D with D the cylinder diameter and P the distance between two centers of in-line neighbor cylinders in the square array; Scruton number or mass-damping parameter

Sc = (2πξ m) / ( ρ f D 2 )
with ξ the solid reduced damping in presence of flow, m the solid mass and ρ f the fluid density; The expression of f s defined in Eq. ( 10) is :

f s = k |Ω s | (x -x 0 ) (13)
with |Ω s | the volume of the solid, and x 0 the solid position at rest. Thus,

F s (t) = ∫ Ωs(t) f s dx = k (x G (t) -x 0 ) (14)
with x G the coordinate of the gravity center of the rigid body. For the studied case, f f = 0.

Using the fluid, solid and geometric characteristics, the dimensionless formulation of system (8) leads to the definition of two significant nondimensional parameters, the Reynolds number and the reduced velocity. This leads to a nondimensional form for the ROM as follows:

for each n = 1, . . . , N, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N ∑ i=1 da i d t Ãin + N ∑ i=1 N ∑ j=1 a i a j Bijn + N ∑ i=1 a i Cin + Dn + Ẽn = 0 D ( ũ) = 0 in Ωs ( t) ∂1 Ωs ∂ t + ũ • ∇1 Ωs = 0 (15)
where:

Ãin = ∫ Ω h ( x, t) Φi • Φn dx Dn = ∫ Ω 1 Ωs Tr [ D(Λ) D( Φn ) ] dx Bijn = ∫ Ω h ( x, t) ( Φi • ∇) Φj • Φn dx Ẽn = 1 M 4π 2 U 2 R ∫ Ω ( x -xG ) Φn dx Cin = 2 Re ∫ Ω g ( x, t) Tr [ D( Φi ) D( Φn ) ]
dx φ denotes the dimensionless form of the corresponding variable or operator φ.

POD basis generation

For POD basis generation HF simulations are performed on this configuration for several values of Reynolds number and reduced velocity. With the previous approach, each set of parameter values corresponds to an optimal POD-ROM. In the fluid domain, the HF computation involves a FV solver with colocalized pressure and velocity fields and a fractional time step approach as described in [START_REF] Archambeau | Code Saturne: a finite volume code for the computation of turbulent incompressible flows[END_REF]. An ALE formulation based on an elliptic equation for grid deformation is involved to account for solid wall motion and the solid dynamics is modeled by using an oscillator equation solved by using a Newmark time scheme. A relaxation is introduced in order to ensure the stability of the explicit solver used to describe the dynamics and kinematics boundary conditions at the time-dependent fluid structure interface. An example of flow velocity field is reproduced in Fig. 5. This field was evaluated on the original reference grid. Then before ROM generation, it is interpolated on a cartesian coarser grid as displayed in Fig. 1.

ROM accuracy and control

Results provided by the reduced-order models of FIV in cylinder arrangements are displayed below. An example of mean flow velocity u and first mode shapes Φ ′ i (for i = 1 to 3) is plotted on Fig. 6 in a two-dimensional configuration involving a periodic cylinder arrangement under cross flow. The mean flow velocity does not allow to capture the velocity of the rigid body. This velocity is captured using the following modes which are non zero on the area traveled by the cylinder. These modes are also needed to have the velocity fluctuations. For example, the mode 3 has to be used to capture the fluid flow and the solid velocity for cylinder position near the maximum amplitude.

a ′ i = (u ′ , Φ ′ i ) (16)
Good adequacy has been found between the POD-ROM solution and the reference one. Finally the solution is reconstructed over the whole domain using the decomposition:

u(x, t) = u(x) + N ∑ n=1 a ′ n (t)Φ ′ n (x)
The reconstruction of the solution is compared to the solution of the reference system over the interpolation domain. The flow velocity field as well as the solid displacement are displayed in Fig. 8. With a truncate threshold fixed at 99%, the reduced order model solution reproduces reference values with a good agreement. The evolution of the relative error on the solution according to the mode number of the reduced model is given in Fig. 9. 10 modes are sufficient to stabilize the error. The size of the ROM is very small compared with the HF model. The sensitivity to grid refinement is also investigated. When increasing the number of cells of the cartesian grid, the slope is improved and the accuracy of the model for a given mode number is better. For the present study,refine the grid more than a 50 × 50 grid does not improve the quality of the ROM.

The Fig. 3 shows a good description of the geometry with this grid. As displayed in Fig. 9, this grid is sufficient to capture the physics of the phenomena.

ROM stability

The sensitivity to penalization factor is plotted on Fig. 10. The stabilization compensates the effect of basis truncation which ensures a preservation of most energetic mode development but in the same time suppresses modes responsible for dissipation [START_REF] Cazemier | Proper Orthogonal Decomposition and Low Dimensional Models for Turbulent Flows[END_REF][START_REF] Cazemier | Proper orthogonal decomposition and law dimensional models for driven cavity[END_REF][START_REF] Rempfer | Evolution of three-dimensional coherent structures in a transitional flat plate boundary layer[END_REF]. In the elementary configurations to be considered, an optimum is obtained for a solid viscosity close to the fluid value.

Performances and CPU time reduction

The predictive capability of ROM over a long period of time is pointed out in several configurations for a given Reynolds number and several values of reduced velocity. The phase portrait of the dynamical response is displayed in Fig. 11 in three cases: for reduced velocity below, close to and above the critical value of the system dynamical stability limit. In each case the solution predicted by the ROM is in good agreement with reference [START_REF] Bourguet | Reduced-order modeling of transonic flows around an airfoil submitted to small deformations[END_REF] and Barone et al. (2009a).

In terms of computational time the gain is displayed on Fig. 12. Comparisons of CPU times are performed by solving the ROM and High Fidelity (HF) model. In the present study, the HF model involves an ALE formulation of Navier-Stokes equations on a moving grid. The interpolation step from the moving grid to the reference grid as well as the POD basis computation are excluded. These two steps are considered as offline steps and only the online steps are compared. The gain in CPU time is expressed in % of gain per iteration, processor and cell. As shown in Fig. 12, the gain of CPU time is about of a factor 100 in the present configurations.

Small ROM parameter perturbations

Principle of sensitivity analysis

Here it is proposed to study the behavior of the POD-ROM when one considers different values of parameters than those used to compute the POD basis. Considering the problem (P λ ):

Find u in space H for parameter λ such that : When a POD basis obtained for a parameter value λ 0 is used for a different value λ, Akkari et al. (2014b, a) proposed to bound the error by the number of POD modes involved in the POD-ROM : where u λ is the approximation of the problem (P λ ) solution and ũ is given by :

⎧ ⎨ ⎩ ∂u ∂t = f (u, λ) u(0) = u 0 (17)
|| u λ -ũ ||≤ g(N)+ | λ -λ 0 | γ h(N) (18)
ũ = N ∑ n=1 a λ,λ 0 n Φ λ 0 n or ũ = u λ,λ 0 + N ∑ n=1 a ′ n λ,λ 0 Φ ′ n λ 0 (19) a λ,λ 0 n
is obtained by solving the POD-ROM resulting from the projection of the problem (P λ ) Eq. ( 17) on

( Φ λ 0 )
.

It is expected that for parameter values not too far from λ 0 , the ROM provides an accurate solution for a number of POD vectors N not too big.

Numerical results

In what follows, pitch ratio, Scruton number and mass ratio are given and the system is characterized by two parameters : Reynolds number and reduced velocity. The purpose is to evaluate the influence of small variations of these both parameters on the accuracy of the reduced model. Values of Reynolds number and reduced velocity are changed and reduced order model associated to each case is used for other cases.

In the present study, the Reynolds number is modified by changing the inlet velocity, and the Reduced velocity by changing the inlet velocity and the spring stiffness. Therefore, the initial values of the a λ,λ 0 n (or a ′ n λ,λ 0 ) time coefficients of the ROM have to be modified. We note Γ the part of Γ f with non zero Dirichlet boundary condition, and x Γ a point of this boundary. The Reynolds number and the reduced velocity can be written :

Re = ρ f u (x Γ ) µ f U R = u (x Γ ) f s D (20)
At the initial step, one gets :

u λ (x Γ , 0) = u λ,λ 0 (x Γ ) + N ∑ n=1 a ′ n λ,λ 0 (0) Φ ′ n λ 0 (x Γ ) (21)
It is proposed to use a scaling coefficient α defined by :

α(λ, λ 0 ) = u λ (x Γ , 0) u λ 0 (x Γ , 0) . (22)
Then

u λ,λ 0 = α(λ, λ 0 )u λ 0 (23)
and

a ′ n λ,λ 0 (0) = α(λ, λ 0 )a ′ n λ 0 (0) (24)
Several configurations for stable and unstable dynamical behaviors are considered. Values of parameters of the configurations to be considered are described in Table 1. Three values of the Reynolds number are considered : 480,612 and 818. In that range of perturbation the flow regime is not changing. Three values are considered for reduced velocity too : a dynamical stable case (U r = 1), a second dynamical stable case close to the critical stability (U r = 2.44) and an unstable case (U r = 3.27). Sensitivity to Reynolds number deviation is displayed in Figs. 13 and 14 respectively in dynamically unstable and stable cases. Reduced velocity deviation is investigated in Fig. 15.

POD basis associated with stable case is used for prediction of unstable solution and conversely. Finally both parameters Reynolds number and reduced velocity are changed simultaneously in Fig. 16. For each case, good behavior of the rigid body has been obtained by the ROM. It means that the ROM built for a parameter couple can be used for another near values of parameters. In the case treated, the ROM is enough robust to model another behavior, i.e from stable to unstable or conversely. Change only one parameter or two in the same time does not affect the results. Increasing number of modes improves accuracy of results.

Conclusion

The multiphase-POD method has been applied for fluid structure interaction study in case of a 2D model of tube array. The proposed reduction method is efficient for velocity field prediction in a domain involving flows in interaction with solid dynamics. It relies on a standard POD-Galerkin projection method combined with an interpolation over a non-moving cartesian grid and a penalization method for dealing with heterogeneity of medium in the whole domain with possible unsteadiness. Both fluid and solid may feature a dynamical behavior. The robustness of the multiphase-POD methods is investigated in order to make possible real-scale real-time parameter sensitivity analysis and preserve the optimality of the generated ROM in terms of error control and computational time reduction. Good agreement has been found between the ROM solution and HF model, even if the ROM is used for another set of parameters than these used to build the POD basis. An extension of the method is now under development to account for turbulence and thermal effects.

The same procedure as that presented in Section 1.3 leads to the following ROM:

for each n = 1, . . . , N, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ N ∑ i=1 da ′ i dt A ′ in = N ∑ i=1 N ∑ j=1 B ′ ijn a ′ i a ′ j + N ∑ i=1 C ′ in a ′ i + D ′ n + E ′ n D(u) = 0 in Ω s ∂1 Ωs ∂t + u • ∇1 Ωs = 0 (27)
with :

A ′ in = ∫ Ω ρΦ ′ i • Φ ′ n dx B ′ ijn = - ∫ Ω ρ(Φ ′ i • ∇)Φ ′ j • Φ ′ n dx C ′ in = - ∫ Ω ρ [ (u • ∇)Φ ′ i + ((Φ ′ i ) • ∇)u ] • Φ ′ n dx -2 ∫ Ω µTr [ D(Φ ′ i )D(Φ ′ n ) ] dx -2 ∫ Γ f µD(Φ ′ i )Φ ′ n ndγ D ′ n = - ∫ Ω ρ(u • ∇)u • Φ ′ n dx -2 ∫ Ω µTr [ D(u)D(Φ ′ n ) ] dx - ∫ Ω 1 Ωs Tr [ D(Λ)D(Φ ′ n ) ] dx E ′ n = ∫ Ω 1 Ωs f s Φ ′ n dx where ρ = ρ f ( 1 -1 Ωs ) + ρ s 1 Ωs , µ = µ f ( 1 -1 Ωs ) + µ s 1 Ωs ,
and Λ is the Lagrange multiplier.

Appendix B. Multiphase formulation for deformable structure with small deformations

The case of a domain containing a fluid and a solid (with a linear elastic behavior) is considered. The aim is to propose a monolithic formulation of the system. Small deformations of the solid are superposed by large displacements of a rigid solid. The extension to large deformations and large displacements of the solid are considered in [START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF].

The solid is first considered as a structure with viscoelastic characteristics. As [START_REF] Hachem | Stabilized finite element method for incompressible flows with high Reynolds number[END_REF] and El Feghali et al.

⎧ ⎨ ⎩ ρ f ∂u f ∂t + ρ f (u f • ∇u f ) -∇ • σ f = f f in Ω f ∇ • u f = 0 in Ω f (28)
with associated boundary conditions, where ρ f is the fluid density, p f the fluid pressure term and σ f the stress tensor in the fluid part, defined by:

σ f = -p f I d + 2µ f D(u f ) (29)
µ f is the fluid dynamical viscosity and I d represents the identity tensor and D(u f ) is the deformation rate tensor, itself defined by

D(u f ) = 1 2 ( ∇u f + ∇ T u f )
.

The solved equation in the solid for u s with u s (t) = ∂d ∂t , where d is the solid displacement, is:

ρ s ∂u s ∂t + ρ s (u s • ∇u s ) -∇ • σ s = f s in Ω s (30)
ρ s is the solid density. The solid stress tensor σ s is defined by:

σ s = 2αD(d) + βTr(D(d))I d (31)
where α and β are the Lamé constants. This leads to its time derivate form σs :

σs = 2αD(u s ) + βTr(D(u s ))I d (32)
If we consider the hypothesis of small deformations, we can define the linearized Lagrangian deformation rate tensor as:

D(u s ) = 1 2 (∇u s + ∇ T u s ) (33)
To Eqs. ( 28) and ( 30), we have to add continuity equations for velocity and stress fields (n the external normal unity vector):

{ u s| Γ i = u f | Γ i σ s| Γ i • n = -σ f | Γ i • n (34)
We can choose to rewrite Eq. ( 32) as:

σs = 2αD(u s ) -ṗs I d (35)
with ṗs = -βTr(D(u s ))I d and we note τ = 2αD(u s ). By coupling Eqs. ( 28) and ( 30) and with the stress tensors definition, we deduce the following total strong formulation:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ρ ∂u ∂t + ρ(u • ∇u) -∇ • (-pI d + 2µ f 1 Ω f D(u) + 1 Ωs τ ) = f in Ω ∇ • u + 1 κ ṗs = ∇ • u + 1 κ [ ∂p s ∂t + u∇p s ] = 0 in Ω s ( 36
)
where p = p f 1 Ω f + p s 1 Ωs and κ = 1 Ωs β, κ represents a compressibility coefficient. Taking into account Eq. ( 36), a total weak form of Navier-Stokes equations in Ω reads: Find u ∈ H(Ω), p ∈ L 2 (Ω) and τ ∈ L 2 (Ω) N×N such as: ∀u * ∈ {v|v ∈ H(Ω), u = 0 on ∂Ω \ Γ i } and for q ∈ L 2 0 (Ω) = {r ∈ where Φ n , Ψ m are elements of POD basis respectively for the velocity and the pressure fields for each n and each m, and a n (t), b m (t) are time coefficients. The final low-order dynamical system is the following, for each Φ i , i = 1, . . . , N and each Ψ i , i = 1, . . . with the system coefficients:

A ni = ∫ Ω ρΦ n • Φ i dx C u nli = ∫ Ω (Φ n • ∇Φ l ) • Φ i dx B p 1,mi = ∫ Ω Ψ m ∇ • Φ i dx B u 1,ni = ∫ Ω µD(Φ n ) : D(Φ i )dx B u 2,ni = ∫ Ω α1 Ωs D(Φ n ) : D(Φ i )dx E i = ∫ Ω f Φ i dx B p 2,ni = ∫ Ω (∇ • Φ n ) • Ψ i dx B κ mi = ∫ Ω 1 κ Ψ m • Ψ i dx C p nmi = ∫ Ω 1 κ (Φ n • ∇Ψ m ) • Ψ i dx
Taking into account the viscoelastic characteristic of the solid, this formulation has the advantage to offer a quasi-classical Navier-Stokes dynamical system to solve, with only a penalization term for the solid viscosity, and the introduction of τ leads to an integro-differential equation.

Fig. 1 .

 1 Fig. 1. Fluid structure interaction scheme.

Fig. 2 .

 2 Fig. 2. Examples of grids for moving (left) and reference (right) domains to be introduced in order to account for time dependency of fluid-solid interface in a configuration involving a periodic cylinder arrangement with a moving central cylinder.

Fig. 3 .

 3 Fig. 3. Space-and time-dependent characteristic functions in the whole domain equal to 1 in the fluid and 0 in the solid domains for fluid and conversely for solid for coarse with 25 × 25 cells (left), refined with 50 × 50 cells (middle) and fine with 200 × 200 cells (right) cartesian grids for a two-dimensional configuration involving a periodic cylinder array.

Fig. 4 .

 4 Fig. 4. Application to tube array.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Velocity field in m.s -1 for Reynolds number R E = 612 and reduced velocity V R = 1. Pitch ratio is 1.44, Scruton number is 0 and mass ratio 7.8.

Fig. 7 .

 7 Fig. 7. Expected values deduced from reference solution versus reduced-order model values of time coefficients a i (for i = 1 to 3) deduced from computation of the dynamical system for case of Figs. 5 and 6.

Fig. 8 .

 8 Fig. 8. Predicted flow velocity field on the cartesian grid using the POD-ROM (left). Solid displacement (right): reference solution (dots) versus predicted solution (lines) for case of Figs. 5 and 6.

Fig. 9 .

 9 Fig. 9. Norm of the relative error on the solution according to number of modes of truncation from 1 to 30 (left) and to grid refinement for grids involving 25 × 25 to 200 × 200 cells (right) for case of Figs. 5 and 6.

Fig. 10 .

 10 Fig. 10. Norm of relative error on the solution according to solid viscosity for case of Figs. 5 and 6.

Fig. 11 .

 11 Fig. 11. Phase portrait of reference (lines) and ROM (dots) solutions for Reynolds number R e = 612 and reduced velocity U r = 1 (top left), U r = 3.273 (top right) and U r = 2.448 (bottom).

Fig. 12 .

 12 Fig. 12. Linear evolution of the gain in computational time compared to the reference simulation according to mode number of the reduced order model for case 3 of Fig. 11.

Fig. 13 .

 13 Fig. 13. Small perturbation of Reynolds number: case 1 (left), case 2 (right).

Fig. 14 .

 14 Fig. 14. Small perturbation of Reynolds number: case 3 (left), case 4 (right).

Fig. 15 .

 15 Fig. 15. Small perturbation of reduced velocity: case 5 (left), case 6 (right).

Fig. 16 .

 16 Fig. 16. Small perturbation of both Reynolds number and reduced velocity: case 7 with 10 (left) and 15 modes (right).

(

  2010) propose, one can write a monolithic formulation of Navier-Stokes equations in Ω. Hachem et al. (2010) and El Feghali et al. are interested by this monolithic formulation in the framework of finite-element models. If one considers the equation governing the fluid domain, one gets the classical incompressible Navier-Stokes equations:

  al. (2010) choose to introduce a viscosity µ s in the solid as a penalization term. Thus, this equation becomes: µ = µ s 1 Ωs + µ f (1 -1 Ωs ). In the case of the POD, trial functions u * are the POD modes. Global velocity and pressure fields are decomposed along space and time as: t)Φ n (x) p(x, t) = M ∑ m=1 b m (t)Ψ m (x) (39)

Table 1

 1 Parameter values used for ROM interpolation with λ 1 and λ 2 respectively the Reynolds number and the reduced velocity of the problem.

	Case	Values of	POD basis	ROM
		parameters for	generation	generation
	1	λ 1	612	818
		λ 2	3.27	3.27
	2	λ 1	818	612
		λ 2	3.27	3.27
	3	λ 1	480	612
		λ 2	1.00	1.00
	4	λ 1	612	480
		λ 2	1.00	1.00
	5	λ 1	612	612
		λ 2	2.44	3.27
	6	λ 1	612	612
		λ 2	3.27	2.44
	7	λ 1	818	612
		λ 2	3.27	2.44

Acknowledgments

This work was performed thanks to the support of DCNS Research, EDF R&D and CNRS. A special thank is addressed to Gaetan Mangeon at EDF R&D for its helpful contribution.

Appendix A. POD-Multiphase ROM for fluctuating velocity field

It has been observed that in fluid mechanics, the first POD mode is similar to the average velocity field and capture the most important part of the energy. Then, a practical approach consists in writing the POD-ROM for the fluctuating velocity field u ′ such as:

where u designates the mean flow velocity field and :

The Navier-Stokes equations for a domain Ω i , i = f , s become :


M J Cánovas
email: canovas@umh.es

A Y Kruger
email: a.kruger@federation.edu.au

H X Phu
email: hxphu@math.ac.vn

Michel Théra
email: michel.thera@unilim.fr

Marco A López

• Michel

Marco Antonio

López Cerdá

August

Segundo Gutiérrez

Cabria He

of Continuous Optimization in Spain

Keywords: assignment problems, convex optimization, graphs, scheduling, game theory, and variational analysis

research documents, whether they are published or not. The documents may come

He has been the leader of many national and international research projects, both theoretical (basic research) and applied (for the university community services, the city and province of Alicante, and bank branches networks).

Marco has written several books. His book Linear Semi-Infinite Optimization [START_REF] Goberna | Linear semi-infinite optimization[END_REF], co-authored with Miguel Angel Goberna and published by John Wiley & Sons in 1998, has become a classical monograph in the field.

Marco's research activity has been increasing with time: out of his 146 publications appearing in MathSciNet more than 60 papers have been published in the last decade. These are mostly papers in top journals like, for instance, Mathematical Programming, Mathematics of Operations Research, SIAM Journal on Optimization, Journal of Optimization Theory and Applications, European Journal of Operational Research, Optimization. Marco's collaborative ability is remarkable. The list of his co-authors includes 60 names.

Pioneer among Spanish mathematicians of international caliber in the field of operations research, Marco jointly with his PhD student at that time Enriqueta Vercher were the first Spanish authors to publish in Mathematical Programming in 1983 (see [START_REF] López | Optimality conditions for nondifferentiable convex semi-infinite programming[END_REF]).

We want to point out Marco's service to the OR community. He served as coeditor-in-chief of TOP (the OR journal of SEIO, the Spanish Society of Statistics and Operations Research) from 2000 to 2007. He serves as Associate Editor of several international journals, including International Abstracts of Operations Research, Journal of Industrial and Management Optimization, Journal of Optimization Theory and Applications, Set-Valued and Variational Analysis, RACSAM, etc. Marco was elected vice-president of SEIO in 1985 and 1998, elected president of SEIO in 1986, and served as chair of the Managing Board of EUROPT, the EURO Working Group on Continuous Optimization, since 2008 to 2010. He has been organizer or co-organizer of several important national and international conferences.

From July 2008 to April 2012 Marco López was the coordinator of i-MATH Consolider, a huge project of the Spanish Ministry of Education for all areas of mathematics, gathering more than 300 research groups, and with a budget of 7.5 million euros.

Marco has enormous general knowledge and is a very versatile person. He is a strong advocate for improving the quality of education and mathematics research at Spanish universities, and promoting the advantages of research for the society in general.

We would like to conclude this part by saying that besides every great man, there is an even greater woman, and Marco (as he said several times) is the first to recognize the amazing contributions of his wife, Maria Pilar, in enabling a fulfilled academic life and family life with 2 children (Pili and Marco) and four grandchildren!

Some key qualities of Marco as observed by his colleagues

The following paragraphs quote verbatim the words of some colleagues. They were written -and gathered in an unpublished book of 2010 -with the aim of paying tribute to Marco on the occasion of his 60th birthday. Now, we take up these words to illustrate a selection of the key personal and professional attributes of Marco.

-Knowledge of mathematics "... At this occasion, I further realized, but this time at a more personal level, the deepness of Marco's knowledge in mathematics which covers so many areas such as convex analysis, functional analysis, optimization from theoretical formalism to numerical methods, as well as applications. .. In those days, research was not very popular within the Spanish OR community but since the very beginning my generation realized that things were about to change and we fixed our scientific benchmark in a few relevant researches that had gained respect for their doubtless trajectory and excellence... There are a few names that would deserve to be mentioned but in this occasion I shall only mention professor López-Cerda, our friend Marco... Marco has been one of the main responsible persons for the current international character of TOP, the Spanish OR journal." Justo Puerto, University of Seville, Spain. -Sense of humor "... Moreover, Marco is an outstanding speaker, not only when he speaks of mathematics. It is impossible to be bored when he is around, He has a contagious humor sense and has a special ability to mix it with work, with people, with live." Ignacio García Jurado, Coruna University, Spain.

Research contributions

Semi-infinite programming

The main research interests and contributions of Marco López have been in the area of semi-infinite programming (SIP). The first series of papers co-authored with Miguel Ángel Goberna, Jesus Pastor and Enriqueta Vercher, in the period 1980-2000, dealt with the geometry of semi-infinite linear inequality systems, as well as with the the optimization theory of linear semi-infinite optimization problems, for the primal and the dual problem (in the sense of Haar). A special attention was paid to deriving semi-infinite versions of theorems of alternative (by Farkas, Gale, Gordan, Motzkin, etc.); and topics like consistency, redundancy, discretizability and reducibility, which are essential for the numerical approach to solving linear SIP problems, were extensively studied. The results obtained in these papers, most of which are co-authored with Miguel Ángel Goberna, were included in the monograph Linear Semi-Infinite Optimization [START_REF] Goberna | Linear semi-infinite optimization[END_REF]. Among these contributions, we higlight below some papers which seem to us to be the most relevant. The choice is very subjective.

-In [START_REF] Goberna | Optimal value function in semi-infinite programming[END_REF], by means of an appropriate use of the subdifferential mapping, the connections between duality theory, solvability properties and convergence properties of sequences of discretized linear semi-infinite programs are made explicit.

-In [START_REF] Goberna | A theory of linear inequality systems[END_REF], the behavior of consistent semi-infinite linear inequality systems is analyzed with respect to many features of the solution set such as boundedness, dimension, description of its boundary, redundant inequalities, minimality, finite reduction, etc. -Duality in linear SIP, and its relationship with discretizability and reducibility, was deeply studied in [START_REF] Goberna | On duality in semi-infinite programming and existence theorems for linear inequalities[END_REF], where it was proved that the Farkas-Minkowski property is crucial for the so-called uniform duality.

-In [START_REF] Fajardo | Locally Farkas-Minkowski systems in convex semi-infinite programming[END_REF], the powerful tool of the Farkas-Minkowski constraint qualification is exploited in the framework of convex SIP, exploring its relationship with stronger qualification conditions. -Infinite-dimensional versions of the Farkas-Minkowski constraint qualification are established, and their role in the optimality conditions in convex optimization is clarified in [START_REF] Dinh | New Farkas-type constraint qualifications in convex infinite programming[END_REF] and [START_REF] Goberna | Necessary and sufficient constraint qualifications for solvability of systems of infinite convex inequalities[END_REF].

Algorithms

Marco López has made significant contributions to the algorithmic solution of SIP problems. We comment on the most relevant ones.

-In [START_REF] Anderson | A: Simplex-like trajectories on quasi-polyhedral convex sets[END_REF] two simplex-like methods are proposed for solving linear SIP's whose feasible regions are quasi-polyhedral sets, that is, sets whose intersection with an arbitrary polytope is a polytope. -In [START_REF] Zhang | A new exchange method for convex semi-infinite programming[END_REF], an exchange-type method for convex SIP is proposed with the feature that no global optimization needs to be carried out at each iteration to detect the (almost) most violated constraints.

-In [START_REF] Anderson | A: Simplex-like trajectories on quasi-polyhedral convex sets[END_REF], the authors formulate the strong separation of two sets in a normed space as an infinite-dimensional optimization problem, and give a method based on an approximation scheme which solves a SIP problem in each iteration without requiring any global optimization.

-In [START_REF] Auslender | Penalty and smoothing methods for convex semi-infinite programming[END_REF], penalty and smoothing methods (of Remez-type) for solving min-max convex SIP problems are proposed and their convergence is proved, whereas in [START_REF] Auslender | Comparative study of RPSALG algorithms for convex semi-infinite programming[END_REF] some computational estimates with these methods, for comparison purposes, is provided.

-In [START_REF] Aragón-Artacho | Local convergence of quasi-Newton methods under metric regularity[END_REF], some quasi-Newton methods to solve generalized equations in Banach spaces are examined under metric regularity.

Stability in SIP

Starting in 1995, in main collaboration with Maria Josefa Cánovas, Miguel Ángel Goberna, Juan Parra, Maxim Todorov, Francisco Javier Toledo and Virginia Vera de Serio, Marco López has published a series of papers on qualitative stability in SIP (linear and convex), whose general aim is to provide conditions under which sufficiently small perturbations of the nominal problem provoke small changes in the optimal value, the optimal set and the feasible set. The stability notions involved in this analysis are lower and upper semi-continuity (in Berge and Hausdorff senses), ill-posedness, etc.

-In [START_REF] Goberna | Stability theory for linear inequality systems[END_REF], the Berge lower-semicontinuity of the feasible set mapping in linear SIP is completely characterized. -In [START_REF] Goberna | Stability theory for linear inequality systems II: upper semicontinuity of the solution set mapping[END_REF] the upper-semicontinuity of the feasible set mapping is studied; its complete characterization is given in [START_REF] Cánovas | Upper semicontinuity of the feasible set mapping for linear inequality systems[END_REF] through the notion of reinforced system. -In [START_REF] Cánovas | Stability and well-posedness in linear semi-infinite programming[END_REF], the lower/upper semicontinuity of the optimal set mapping and the wellposedness property (in the sense of Hadamard) are characterized for linear SIP's, whereas these results are extended to the convex setting in [START_REF] López | Stability of the feasible set mapping in convex semi-infinite programming[END_REF].

-In [START_REF] Cánovas | Ill-posedness with respect to the solvability in linear optimization[END_REF] and [START_REF] Cánovas | Sufficient conditions for total ill-posedness in linear semi-infinite optimization[END_REF], the ill-posedness in linear SIP with respect to various properties such as consistency, solvability, etc. (i.e., the boundary separating problems enjoying such a property and problems for which the property fails) is analyzed in detail.

-Qualitative stability when only some of the coefficients are perturbed and when the mappings are restricted to their domains, is the main subject of [START_REF] Daniilidis | Lower semicontinuity of the feasible set mapping of linear systems relative to their domains[END_REF] and [START_REF] Daniilidis | Stability in linear optimization under perturbations of the left-hand side coefficients[END_REF].

Infinite-dimensional counterparts of this qualitative analysis are performed in [START_REF] Dinh | On the stability of the feasible set in optimization problems[END_REF] and [START_REF] Dinh | On the stability of the optimal value and the optimal set in optimization problems[END_REF].

Lipschitz-type stability properties

Also known as quantitative stability, as these properties aim to measure the level of stability of the mappings of interest (feasible set, optimal set, and optimal value mappings). These properties have important consequences in the overall stability of the SIP, and affects also the numerical complexity of algorithms. The properties studied by Marco López and his co-authors are, mainly, pseudo-Lipschitzness, calmness, and distance to ill-posedness. Special attention is paid to characterizing, or estimating, the corresponding moduli as they quantify the level to which such a stability property holds.

-In [START_REF] Cánovas | Metric regularity of semi-infinite constraint systems[END_REF][START_REF] Cánovas | Stability of indices in KKT conditions and metric regularity in convex semi-infinite optimization[END_REF][START_REF] Cánovas | Metric regularity in convex semi-infinite optimization under canonical perturbations[END_REF][START_REF] Cánovas | Variational analysis in semi-infinite and infinite programming, I: Stability of linear inequality systems of fesible solutions[END_REF] ( co-authored with M.J. Cánovas, A. Dontchev, A. Hantoute, D.

Klatte, B. Mordukhovich, and J. Parra), the pseudo-Lipschitz property of the feasible and the optimal set mappings (or metric regularity of the inverse mappings), both in the linear and in convex settings, are approached; the last one is developed in infinite-dimensional spaces exploiting the properties of coderivatives, as it is also done in [START_REF] López | On coderivatives and Lipschitzian properties of the dual pair in optimization[END_REF]. -In the field of calmness, the main contributions of Marco López ( jointly with M.J. Cánovas, R. Henrion, A. Kruger, J. Parra and M. Théra) are as follows: an expression for the calmness modulus of the mapping in linear programming under canonical perturbations (objective function and right-hand side of the constraints), involving limits of subdifferentials [START_REF] Cánovas | Calmness modulus of linear semi-infinite programs[END_REF]; left-hand-side perturbations of the constraints system added into the analysis in [START_REF] Cánovas | Calmness of the feasible set mapping for linear inequality systems[END_REF]; and the outer limits of subdifferentials of max-functions and calmness moduli for feasible and optimal set mappings dealt with in [START_REF] Cánovas | Outer limit of subdifferentials and calmness moduli in linear and nonlinear programming[END_REF]. -In [START_REF] Cánovas | Distance to ill-posedness and the consistency value of linear semi-infinite inequality systems[END_REF][START_REF] Cánovas | Distance to solvability / unsolvability in linear optimization[END_REF][START_REF] Cánovas | Distance to ill-posedness in linear optimization via the Fenchel-Legendre conjugate[END_REF], the distance to ill-posedness (in terms of the distance to infeasibility and to unsolvability) is studied.

Subdifferential calculus

The subdifferential of the supremum of convex functions is a key tool in convex analysis. In fact, Fritz-John and KKT-type optimality conditions for (semi-infinite) convex optimization can be approached by appealing to the subdifferential of the supremum of the constraint functions. In [START_REF] Hantoute | Subdifferential calculus rules in convex analysis: a unifying approach via pointwise supremum functions[END_REF] ( co-authored with Hantoute and Zȃlinescu) the first general formula for the subdifferential of the supremum function was established. Later on, in a series of papers co-authored with R. Correa and A. Hantoute [START_REF] Correa | Towards supremum-sum subdifferential calculus free of qualification conditions[END_REF][START_REF] Correa | Valadier-like formulas for the supremum function I[END_REF][START_REF] Correa | Valadier-like formulas for the supremum function II: The compactly indexed case[END_REF][START_REF] Correa | Moreau-Rockafellar-type formulas for the subdifferential of the supremum function[END_REF], various formulas for the subdifferential of the supremum were established, extending and improving some classical results (due to Brøndsted, Valadier, Volle, etc.) as some standard continuity assumptions are dropped or weakened. In the last two quoted papers, Fritz-John and KKT optimality conditions for convex SIP are derived under very weak assumptions by applying the obtained subdifferential rules.

Game theory

The (unpublished) doctoral thesis of Marco López was on cooperative game theory, and it focused on Markov processes of cooperation and the valuation of a game. During most of his career, Marco has been fond of the game theory and delivered many courses and seminars on this topic. We quote some of his contributions.

-In [START_REF] López | Convex semi-infinite games[END_REF], a generalization of the two-person zero-sum semi-infinite game introduced by A. L. Soyster [Management Sci. 21 (1974/75), no. 7, 806-812] is provided. The pure strategies for player I are picked from an infinite family of convex functions, whereas the set of mixed strategies available to player II is a closed convex set. -The paper [START_REF] López | Dynamic strategic models for markets with fixed prices. (Spanish) Trabajos Estadíst[END_REF] formulates a generalization of the model of M. Balch [Differential games and related topics, Proc. Internat. Summer School, Varenna, 1970, North-Holland, Amsterdam (1971)], and analyzes the oligopoly (and the Nash-Cournot and Pareto equilibria) from the point of view of differential games theory. -In [START_REF] Sánchez-Soriano | On the core of transportation games[END_REF] the benefit sharing in holding situations is modeled as a big boss game, and its most relevant properties are examined. -Finally, in [START_REF] López | Convex semi-infinite games[END_REF] the authors prove that the transportation games have nonempty cores, and an LP-procedure describing this solution is shown.

Miscellaneous

Marco López has contributed to many other areas which include: robustness, infinite-dimensional duality theory, stationarity and regularity, error bounds, and Hölder stability.

Marco has published several surveys on SIP. The most cited one is [START_REF] López | Semi-infinite programming[END_REF] (coauthored with G. Still), whereas the most recent survey [START_REF] Goberna | Recent contributions to linear semi-infinite optimization: an update[END_REF] reviews the state of the art in the theory of deterministic and uncertain linear semi-infinite programs, presents some numerical approaches to this type of problems, and describes a selection of recent applications in a variety of fields. Recently Marco López has published another two books: monograph Post-Optimal Analysis in Linear Semi-Infinite Optimization [START_REF] Goberna | Post-optimal analysis in linear semi-infinite optimization[END_REF] with M.A. Goberna, and Nonlinear Optimization [START_REF] Aragón-Artacho | Nonlinear optimization[END_REF] with F.J. Aragón-Artacho, M.A. Goberna and M.M.L. Rodríguez. The last one is a great companion for upper level undergraduate students in optimization.

  From the performance of Marco in the meeting we held in the Executive Council of i-Math 1 , I highlight his good character, his calmness, his ability to negotiate and reach agreements, and his dedication to the job." Eduardo Casas, University of Cantabria, Spain.-Thoroughness in writing and reviewing "The first time that I have been in contact with the mathematical opus of Marco was probably several years ago. I was one of the referees of a paper he submitted to a journal. Then, I understood that I was facing one of the best specialists in the field of linear or convex semi-infinite programming and well posedness problems ... Also, I have asked Marco to act as a referee for some papers that I was handling, the job was done in due time, carefully, with appreciated helpful comments for the authors. I consider Marco as a very good and deep mathematician." Jean-Pierre Crouzeix, Blaise Pascal University, France. -Scientific benchmark "... I met Marco in the mid nineties. We were at a conference... I was curious about the reasons a senior (senior, not old) man could have to be one of the few persons in the audience. Not only that: since he discussed with the speakers after their talks, I realized he was even paying attention! ... Enthusiasm for learning and teaching maths, scientific thoroughness, high estimation of the scientific work of others, and unbounded working capacities are qualities I mentally linked to his face during this conference. Time has given me the chance of meeting Marco

	.. I truly appreciated his
	remarkable abilities, his enthusiasm for investigating new fields, and his sense of
	perfection, all qualities that are very stimulating for his collaborators." Alfred
	Auslender, Institut Camille Jordan, France.
	-Research capability
	"... We all know that Marco is an outstanding mathematician, an international
	leader in semi-infinite optimization and related topics, the founder of the strong
	research Spanish group in these areas. But only those who actually worked with
	him will share my feelings on how deeply he approaches mathematical problems
	trying to understand the very essence of the issue and how strongly he contributes
	in joint projects, from the point of view of both ideas and techniques." Boris S.
	Mordukhovich, Wayne State University, USA.
	-Undergraduate teaching skills
	"... He was a genuine lord. Etched on my memory ... his clear handwriting and
	organization on the blackboard. When writing on it, he sometimes got further to
	see in perspective his writing, checking if what he wrote was clear; if something,
	even a minimal thing like the tail of the 'a' was not completely perfect, then he
	deleted it and wrote it again.

.. I have never had another professor who knew so much from so many things... I have to say that he has been the best teacher I have ever had." Isabel Molina, Universidad Carlos III de Madrid Spain.

-Hospitality "... Of course, I benefited much from our scientific discussions and from your profound knowledge of many areas of mathematics. But today I like to emphasize some other side of your personality: your warm, perfect hospitality. I was several times your guest, and I was impressed and overwhelmed each time. To be your guest means to have a big cultural experience because of many, many impressions of your country and the interesting information on Spanish and local history and traditions. It also means a culinary highlight..." Diethard Klatte, University of Zurich, Switzerland. -Training of young researchers "... Another contribution which is not immediately seen in my point of view, at least by Internet, is the big work of Marco in forming scientific groups and the establishment of recognized centers... I had the opportunity (and the chance) to be one of his postdoctoral students, and I realized that together with the famous scientific personality, he is a very fine person with many wonderful qualities. A student not only learns theorems and corollaries from him but also the love of doing them. He is a very patient person and very gracious at the same time as a good listener..." Abderrahim Hantoute, Center for Mathematical Modeling, University of Chile. -Leadership abilities "... much closer... Saying that he is my scientific benchmark in OR in Spain is an accurate first-order approximation to what I do feel." Emilio Carrizosa, University of Seville, Spain. -Internationalization of the Spanish research "... I still remember my first SEIO meeting in 1988.

It refers to the period when Marco was the coordinator of the huge Spanish project i-Math Consolider .