L Cassayre

G Plascencia

T Marin

S Fan

T Utigard

GAS EVOLUTION ON GRAPHITE AND OXYGEN-EVOLVING ANODES DURING ALUMINIUM ELECTROLYSIS

Keywords: aluminium electrolysis, gas evolution, slotted anode, inert anode

Anode gas evolution, growth and flow behaviour during aluminium electrolysis have been investigated using various experimental techniques, including water modelling, X-ray visualization and direct observation. Video recordings of oxygenevolving anodes (SnO 2 , Cu, Cu-Ni) and carbon anodes were performed in laboratory electrolysis cells of various scales. The water model also investigated the effects of slotted anodes on the gas escape from beneath large anodes. Beneath large horizontal anodes, individual bubbles form and are then subsequently swept away by large sweeping bubbles flowing rapidly beneath the bottom surface. The gas behaviour on oxygen evolving anodes was very different from that with carbon, with the formation of very small bubbles that released very quickly from the anode due to improved wetting of the anode by the electrolyte.

Introduction

In the Hall-Héroult process, aluminium is produced by electrolysis of alumina at 960 °C. The cell contains normally 15 to 20 cm of liquid aluminium which also acts as the cathode. Carbon anodes are suspended in the cryolite based electrolyte and separated from the aluminium pool by a distance of 4 to 5 cm. Modern cells with cell currents above 250 kA are capable of producing aluminium at 13-14 kWh/kg of aluminium and with a current efficiency of 95-97%. These cells represent major improvements in technology over the last two decades. These advancements are based on extensive theoretical and mathematical studies, laboratory and plant investigations, and plant operating experience. However, most of these investigations are based on indirect measurements since very little of what actually takes place within the cell can be observed. To "look" into the process, we have developed various techniques to visualize the gas evolution, including both wetting characteristics and flow behaviour beneath the anodes [START_REF] Cassayre | Visualisation of Gas Evolution on Graphite and Oxygen-Evolving Anodes During Aluminum Electrolysis[END_REF][START_REF] Utigard | Radiographic observations of the Hall-Heroult electrolysis in bench scale cells[END_REF].

The two main objectives of this work were to investigate the effect of slots in the anodes on the gas escape behaviour and to study the gas behaviour using inert anodes. With the use of large anodes, the gas escape from beneath such large horizontal surfaces may lead to excessive electrolyte stirring as well as increased anode voltages due to the large surface coverage by gas bubbles. Slotted anodes have recently been introduced to allow easier gas escape. The use of inert anodes is likely to affect gas evolution and flow behaviour. One aim of the current work was to compare the gas behaviour between oxygen evolving and graphite anodes, and to characterize the bubble size and the bubble layer depth in the inter-electrode area.

While various materials have already been patented for inert anode application, none is yet in commercial use [START_REF] Pawlek | Inert anodes: an update[END_REF]. This is mainly because the severe operating environment and the need of good electronic conductivity limit the number of suitable materials. Most recent research has concentrated on metals, oxide ceramics and cermets. In the present study, copper, copper-nickel alloy (wt% 75/25) and tin oxide materials were tested as oxygenevolving anodes, even if none of them is stable enough to totally resist to corrosion.

Video recordings of electrolysis at various current densities (CD) were performed in laboratory cells, using direct observation from above, see-through cell and radiography techniques [START_REF] Popelar | Anode Effect[END_REF][START_REF] Utigard | Direct observation of the anode effect by radiography[END_REF]. The influence of the alumina content in the electrolyte was studied for graphite and tin-oxide. A comparison of the wettability of the materials by the electrolyte was also carried out to try to explain the very different gas behaviour observed.

Experimental Techniques

The experimental techniques used include: Direct observations from above, Quartz see through cell, X-ray radiography, Hotthermocouple microscope, Hot-stage microscope and Watermodelling.

Direct Observations From Above

Experiments were carried out in graphite crucibles with 1.5 kg of cryolite and placed inside an inconel container and heated inside a furnace (Fig. 1A). The voltage measurements were made directly from the stainless steel buss bars holding the anode and the cathode (crucible). Anode vertical positioning was achieved via a screw-driven linear movement device. The current and voltage signals were transmitted to isolated inputs of the data acquisition. The graphite (Union Carbide AGSR grade) used to make the anodes had a density of approx. 1.7 g/cm 3 . Quartz See Through Cell (Fig. 1B) Electrolysis was also carried out in transparent quartz tubes (∅=54 mm) placed in an electric furnace allowing lateral observations through a hole in the front panel, as described by Haupin et al. [START_REF] Haupin | See-through Hall-Heroult cell[END_REF]. In this see-through cell, the cathode was a graphite cylinder (∅=43 mm, h=15 mm) located at the bottom of the tube, and the anode was a graphite or tin-oxide cylinder (∅=9 mm, h=50 mm), vertically positioned 2 cm above the cathode. A new quartz tube was used for each experiment.

Two different electrolyte compositions were used, based on a mixture of cryolite with AlF 3 (11 wt%), CaF 2 (5 wt%) and Al 2 O 3 . Weight fraction of alumina was either 2.5 wt% (electrolyte A) or 9 wt% (electrolyte B), electrolyte B being saturated in alumina. The mix of powders was pre-melted at 1000°C in graphite crucibles, then cast, crushed and stored dry.

X-ray Radiography (Fig. 1C)

The experiments were carried out in graphite crucibles with inner diameters from 17 to 23 mm with 2 to 9 mm graphite or inert anodes. Images of the aluminium drop and the graphite anode immersed in the electrolyte were formed on a photographic film by a horizontal beam of X-rays. Gas bubbles appear white because they allow X-rays to pass through. More dense materials such as aluminium, appears darker. The experiments were recorded using a TV camera-monitor and a video recorder.

Hot-thermocouple Microscope

A hot-thermocouple method was employed in conjunction with a microscope and a video camera. A very small sample (0.5 g) of electrolyte was attached to and heated by passing current through a Pt/Pt-Rh thermocouple. A pair of platinum electrodes (0.3 mm dia.) was inserted in the sample after melting and electrolysis was initiated and observed using a microscope.

Hot-stage Microscope

A 'hot-stage microscope' was developed in order to allow for bigger samples than what could be achieved by the hotthermocouple method. A 1-3 g sample of cryolite was placed on a boat and heated by a platinum foil and contained in a gas-tight Pyrex container. A pair of graphite electrodes (1-1.3 mm dia.) and a thermocouple could be placed within the sample. A microscope was used to monitor the various phenomena taking place.

Water-modelling

The water model was constructed from acrylic plastic. The purpose of this technique was to study the gas behaviour beneath large horizontal, sloped and slotted anodes. The model represents half of an anode which is suspended inside a tank and can be moved relative to the cathode. Three independent chambers in the anode can be pressurized with gas to force bubbles to form on the face of the brass sheet. The entire model is situated on top of two scissor jacks which allows the angle of the cathode and anode to be adjusted.

Bubbling at Graphite Anodes

At current densities below 0.2 A/cm 2 , bubble birth occurrs on specific nucleation spots on the graphite anode. Increasing CD leads to a large number of nucleation sites. Round-shaped bubbles grow and coalesce until reaching 2-3 mm in diameter, then detach and escape vertically. As shown on Figure 2, the surface coverage of the anode is very high and the small bubbles appears stuck to the anode. Tests in the quartz see-through cell were limited to less than 7 minutes and CD up to 0.5 A/cm 2 , because streamers rapidly arose in the cell and lead to an opaque electrolyte. It was observed that the anode gas was released periodically. At small CDs (up to around 0.5 A/cm 2 ), several bubbles form on the bottom surface of the anode, then grow and coalesce until a single large bubble covers the whole area and is finally released, leaving a clean surface free of gas. Figure 3 shows a bubble trapped under the anode just before release, during see-through cell electrolysis at 0.2 A/cm 2 .

A B C

Figure 1: Schematics of the experimental conditions A: Observation from above B: See-through cell C: X-Ray furnace Generally, the bath tended to wet the graphite anode fairly poorly as compared to the cathode [START_REF] Utigard | Radiographic observations of the Hall-Heroult electrolysis in bench scale cells[END_REF]. Gas was generated over the entire anode/electrolyte interface during electrolysis. The gas bubbles on the side of the anode varied in size up to about 3 mm and the gas bubble under the anode penetrated a maximum of 5 mm into the electrolyte. With increasing CD, the bubbles became smaller and were formed at a higher frequency. This may be due to the increased amount of gas formed, creating more agitation in the melt.

At higher CDs, coalescence didn't happen as much as at low CDs, and the bubbles escaped before covering the whole anode, as illustrated on Figure 4. The average bubble diameter before release, calculated as the mean diameter of 10 bubbles for each CD, is presented on Figure 5. The diameter decreases from around 18 mm at 0.2 A/cm 2 to 7 mm at 1.6 A/cm 2 . The alumina content in the electrolyte doesn't have a major influence on the bubble size. The release frequency of the CO 2 bubbles is strongly related with the current density. Using the X-ray technique it was observed that simultaneous with the occurrence of the anode effect, there was a sudden formation of a gas layer under the anode. This sudden non-wetting occurred at the same time as the voltage increased. It was even possible to promote the anode effect using the crucible as the anode. When the anode effect occurred, the melt became separated from the crucible, and a gap of up to 1 mm was observed between the crucible bottom and the electrolyte. Using the hot-stage microscope, electrolysis was carried out with a 1 mm diameter graphite anode. With 10 wt% alumina very fine bubbles form along the anode surface. When only 1 wt% alumina was used, slightly larger gas bubbles formed, and in pure technical cryolite, large bubbles form at the anode indicating poor wetting of the anode. After the anode effect was initiated, the anode became non-wetted and it was difficult to immerse the anode into the electrolyte. It was like pushing a finger into a blown-up rubber balloon with the electrolyte skin stretching around the tip of the anode.

The water modeling study was used to investigate the gas bubble behaviour on a scale similar to that in industry and also to investigate slotted anodes. By the use of a water model, the effects of anode/cathode distance, cathode inclination, and anode gas evolution rate on the mixing rates of electrolyte in the cell could be investigated. The experiments were conducted by filling the tank with water so that the anode was immersed 1 to 15 cm. The gas flowrate, anode angle, and anode-cathode distance were set and blue dye was added at one end of the tank. Five optical sensors were used to measure the light intensity passing through the water from the back to the front of the tank. When the anode was slightly tilted(a fraction of a degree), large bubbles formed and flowed beneath the anode with a leading edge penetrating up to 2 cm into the electrolyte and a long tail as seen in Figure 6. This was also observed by Fortin et al [7] who explained this behaviour by variations in the hydrostatic pressure due to the motion of the bubble. The results are in general agreement with previous researchers [START_REF] Haupin | See-through Hall-Heroult cell[END_REF][7][START_REF] Welch | Gas driven electrolyte flow[END_REF].

When a vertical slot was made in the anode, it was observed that for slots as narrow as 0.5 cm, the moving gas front entered the slot nearly 100%. The gas would flow along the horizontal anode and then up through the slot and leave the anode at its upper side. It was found that even when the slot was fully submersed in the water, as would be the case for anodes just before they have to be replaced, the gas would flow up into the slot, and then out along the side. For penetration into narrow slots, water and cryolite should behave very similar since the surface tension to density ratio is nearly the same for the two. Industrially, it is expected that with slotted anodes, gas splashing of cryolite will partly fill and block the upper section of the slots. However, this should not cause any problems since the gas will still be able to escape as seen in this study. Another issue is that the heat up of cold anodes will be more uneven with the center section heated the fastest and initially draw more current than the outer and colder section. On the other hand, for a fully heated anode the current profile through the anode will be more vertical, potentially improving the magnetics of the pot. The reduced gas coverage beneath the anode will lower the cell resistance slightly, but this may be counteracted by a slight increase in the total anode resistance.

Bubbling at Inert Anodes

Gas evolution is pretty much the same for tin-oxide, copper and copper-nickel anodes, but was found to be very different as compared to what was observed on the graphite anode. Observations of a tin-oxide anode from above showed tiny bubbles escaping continuously without much coalescence, forming a foam around the electrode, as illustrated in Figure 7. The bubble diameter, measured under the microscope was estimated to be around 0.1 mm, which is 10 to 30 times smaller than on graphite anodes. Some tests were run in see-through cells, but observations were difficult because the froth of bubbles hides the electrode even at small CDs (0.1 A/cm 2 ). This resulted in poor contrast between the electrolyte and the gas. Radiographic observation proved to be more adapted to the study of the gas evolution on oxygen evolving anodes. The anode appears in black and the bath in dark grey, while gas takes the form of a light grey halo around and under the electrode (Figure 8). As almost no coalescence occurs between the tiny bubbles, gas release looks continuous, and the gas foam under the electrode keeps roughly a constant shape at fixed CD. The depth of the gas layer under the anode is plotted in Figure 9 as a function of CD and material. The evolution is almost identical for the three oxygen-evolving electrodes: the bubble layer thickness grows regularly with CD until reaching 2-2.5 mm at 1 A/cm 2 . Increasing CD over 1 A/cm 2 doesn't increase the penetration of the gas under the anode, maybe because the volume fraction of the gas in the bubble layer increases. However, the contrast of the X-ray images wasn't good enough to observe noticeable changes in the brightness of the gas froth, as it could be expected with a higher gas density. Another possibility is that the release velocity of the bubbles increases, leading to a shorter residence time of the bubbles under the anode. Wetting experiments on graphite versus metallic and tin oxide substrate materials showed that on graphite the electrolyte gave contact angles of about 120-130 o , indicating poor wetting. However, metallic substrates and tin oxide were completely wetted by the electrolytes forming contact angles of near 0 o . This dramatic difference most likely is the reason for the dramatic difference in gas evolution behaviour on graphite versus metallic "inert" anodes. The bubbles formed on tin-oxide, copper and copper-nickel, all very well wetted by the electrolyte, will detach easily from their nucleation site, while bubbles formed on carbon are stuck and grow before detachment.

mm

Conclusions

Video film showing various aspects of the Hall-Héroult process has been developed using X-ray radiography, hot-stage microscopy, direct observations from above the bath and water modeling techniques. During anode effects, the anode is nonwetted by the electrolyte with a thin film surrounding the anode.

The characterization of the bubble layer under the anode showed many differences between graphite and metallic anodes. The study of three different oxygen-evolving anodes showed that a good wetting by the electrolyte leads to the formation of tiny bubbles. Contact angle measurements of the system electrolyteanode-oxygen might then allow predicting the gas evolution on other candidate as inert anodes materials without having recourse to heavy equipment like X-ray radiography.

Water modelling of slotted anodes showed that they are effective in releasing gas bubbles from beneath large anodes. Even narrow slots of about 0.5 cm, allowed for rapid removal of gas as it swept across the surface of the anode.

Figure 2 :

 2 Figure 2: Observation from above -side bubbles on graphite anode [∅ anode =5mm, CD=0.2 A.cm -2 , T=960°C, 9%Al 2 O 3 ]

Figure 3 :

 3 Figure 3: Graphite anode during electrolysis in see-through cell [∅ anode =9mm, CD=0.2 A/cm 2 , T=960°C, 9%Al 2 O 3 , ACD=2 cm]

Figure 4 :Figure 5 :

 45 Figure 4: Graphite anode during electrolysis in X-Ray furnace [∅ anode =18mm, T=960°C, 2.5%Al 2 O 3 , ACD=2 cm]

Figure 6 :

 6 Figure 6: Gas bubble flow behaviour under carbon anode

Figure 7 :

 7 Figure 7: Observation from above -side bubbles on tin-oxide anode [∅ anode =5mm, CD=0.2 A/cm 2 , T=960°C, 9%Al 2 O 3 , ACD=2 cm]

Figure 8 :

 8 Figure 8: Tin-oxide anode during electrolysis in X-Ray furnace [∅ anode =18mm, T=960°C, 9%Al 2 O 3 , ACD= 2 cm]

  Test with electrolytes containing either 2.5% or 9% Al 2 O 3 on tinoxide anode didn't reveal any difference in the gas behaviour. The bubble layer under the graphite anode is slightly decreasing with CD from 5 mm to 4.2 mm, confirming both the values and the trend reported by Aaberg et al.[START_REF] Aaberg | The gas under anodes in aluminium smelting cells Part II: Gas volume and bubble layer characteristics[END_REF]. At 1 A/cm 2 , gas penetration under the anode is roughly twice as deep with graphite as with oxygen-evolving anodes, which is important with regards to the anode-cathode optimum distance.

		6.0
		5.5
		5.0
	Bubble layer depth (mm)	2.0 2.5 3.0 3.5 4.0 4.5 1.5									Cu-Ni Cu
		1.0									SnO 2
		0.5									Graphite
		0.0
		0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8
					Current density (A.cm	-2 )
	Figure 9: Influence of CD on bubble layer thickness under various
		anode materials [∅ anode =18mm, T=960°C, ACD=2 cm]

Table I

 I

		. Gas evolution summary
	Anodic Material	Bubbling	Bubble	Bubble layer
			diameter	depth
	Graphite	Periodic release	2-10	~ 5 mm
		of large bubbles	mm
	Oxygen-evolving	Froth of tiny	0.1 mm	~ 2 mm
		bubbles

Acknowledgement

The financial support from the National Science and Engineering Research Counsel of Canada and Pechiney is greatly appreciated.


Simon Bagy
email: simon.bagy@isl.eu

Michel Libsig
email: michel.libsig@isl.eu

Bastien Martinez
email: bastien.martinez@isl.eu

DESIGN AND OPTIMIZATION OF LONG-RANGE GUIDED AMMUNITION

This paper describes the use of optimization approaches to increase the range of near-future howitzer ammunition. The performance of a gliding ammunition concept, developed at the French-German Research Institute of Saint-Louis, is assessed using an aeroballistic workflow comprising aerodynamic characterization and flight trajectory computation. First, a single-objective optimization is run to find the maximal attainable range for this type of projectile. Results show significant range improvement along the optimization process. Then, a multi-objective formulation of the problem is proposed to consider the compromise between range and time of flight.

INTRODUCTION

On the battlefield, howitzer batteries are used for artillery fires on long-distance objectives. As a consequence, their range of fire is of strategic interest and several existing technologies are developed to extend it, such as base-bleed [START_REF] Sahu | Navier-stokes computations of projectile base flow with and without mass injection[END_REF][START_REF] Nietubicz | Navierstokes computations for a reacting, m864 base bleed projectile[END_REF][START_REF] Lyle D Kayser | Flight testing for a 155mm base burn projectile[END_REF] or propulsion using rocket or air-breathing engines [START_REF] Andersson | Solid fuel ramjets for projectile propulsion-summing up of a joint FOA-TNO research project[END_REF]. Still, these devices come at high price, which is often limiting for their use in operational situations. In this context, the development of affordable extended-range ammunitions for current and future artillery systems appears particularly relevant. Several studies have been conducted on this matter at the French-German Research Institute of Saint-Louis (ISL) [START_REF] Decrocq | Aerodynamic study for capability increase of a spin-stabilized projectile fitted with fins[END_REF][START_REF] Decrocq | Aerodynamic prediction of a projectile fitted with fins[END_REF][START_REF] Decrocq | Aerodynamic Study for Capability enhancement of Artillery Ammunition Fitted with Fins[END_REF]. The present work is part of a project focusing on the design of a long-range guided projectile without base-bleed or rocket propulsion [START_REF] Libsig | CS 2019 -long range guided projectile (LRGP)[END_REF][START_REF] Libsig | CS 2020 -long range guided projectile (LRGP)[END_REF]. In particular, this study aims at introducing an optimization approach in the design process.

In recent decades, optimization methods have know a great development thanks to computer science and the constant increase of computational resources. Today, they regroup various types of algorithms, that can be applied on most engineering processes, including aerodynamic design [START_REF] Mohammadi | Applied shape optimization for fluids[END_REF]. In the research field on projectiles, some studies have applied such approaches for ammunition design. For example, Vasile et al. [START_REF] Vasile | A Multi-Disciplinary Approach to Design Long Range Guided Projectiles[END_REF] improved the geometry of a long-range guided projectile using a Particle-Swarm Optimizer (PSO). In particular, the performance of the configurations is evaluated with a Semi-Empirical Aerodynamic Prediction (SEAP) code. Another application can be found in the work of Arnoult et al. [START_REF] Arnoult | Control surface geometry surrogate-based optimization for spin-stabilized projectile course correction[END_REF][START_REF] Arnoult | Modélisation de la trajectoire d'un projectile gyrostabilisé muni d'un dispositif de contrôle[END_REF], that performed optimizations of a trajectory control device for spin-stabilized ammunitions. The optimal geometrical parameters of this device, consisting in an aerodynamic spoiler deployed during the flight, are researched with a Surrogate-Based Optimization (SBO) approach. Moreover, this study uses a Multi-Layer Perceptron (MLP) neural network as surrogate model to predict the aerodynamic effect of the spoiler.

In the following, Section 2 introduces the developed long-range projectile and its specific flight scenario. Then, the optimization approach is detailed in Section 3 with its formulation and workflow. Section 4 provides a deeper explanation of aerodynamic characterization and the use of surrogate modelling. Finally, results obtained with single-objective and multi-objective optimizations are discussed in Sections 5 and 6.

STUDY CASE

The studied configuration is based on a Secant Ogive Cylinder Boat Tail (SOCBT) geometry. The full-scale projectile with a calibre of 155 mm is considered. For long-range applications, the projectile is non-rotating and comprises four tail-fins, unfolded immediately after the muzzle exit. In addition, two canards are deployed during the flight, therefore resulting in the two configurations depicted in Fig. 1.

The projectile is shot from a standard howitzer gun with an initial spin rate supposed to be null. The flight strategy associated to such ammunition concept is depicted in Fig. 2 and can be decomposed in two phases, namely a ballistic phase with tail-fins only, and a gliding phase after the canards deployment. The canards are then used as lifting surfaces but also as flight control steering devices during the gliding flight.

OPTIMIZATION APPROACH

Formulation

The main objective of this study is to improve the range of the gliding ammunition concept, denoted as L f light . This In addition, attention is given to the total flight time of the projectile, which can be of equal interest as the range from the operational perspective. Thus, it is considered as an additional function of interest, and minimizing t f light is the second objective defined for optimizations.

The decision parameters considered in this work consist in both geometrical and flight parameters. The geometrical design parameter is the sweep angle of the canards after deployment, represented by the angle between the SOCBT revolution axis and the canards and named A Can (as depicted on Fig. 3). Among flight parameters, the gun elevation angle, the maximum angle of attack of the projectile and the time of canards deployment can be modified by the optimizer. A summary of the decision parameters and objective functions is given in Tab. 1 and 2 respectively.

Workflow and tools

The functions of interest (summarized in Tab. 2) imply that the optimization workflow has to be aeroballistic, i.e. include aerodynamic projectile characterization and fullflight trajectory computations. The process defined for this purpose is depicted in Fig. 4.

First, the aerodynamic characterization consists in computing the projectile aerodynamic force and moment coefficients with respect to flight conditions and design parameters. Due to the flight scenario, a wide range of Mach numbers is expected. In this situation, Computational Fluid Dynamics (CFD) simulations appear wellsuited, because they allow to study subsonic, transonic and supersonic flows. Still, a large number of simulations is required to characterize one geometry. This step has to be repeated when decision parameters are modified, which results in a significant overall computational cost. For this reason, the CFD simulations cannot be performed in-line during the optimization process. As After aerodynamic characterization, the trajectory simulation computes the two-phase flight path. First, the BALCO trajectory computation tool [START_REF] Wey | Balco 6/7dof trajectory model[END_REF] is used to simulate the ballistic flight phase and integrate the aerodynamic coefficients of the projectile with folded canards. After canards deployment, the flight is assumed to be planar and a Three-Degrees-Of-Freedom (3DoF) step-bystep computation of the gliding trajectory is performed. In addition, a basic guidance law is implemented in the simulation during the gliding phase. The canards operate as control surfaces and their steering is updated at each iteration in the code. The complete trajectory computation process is managed with MATLAB [START_REF]MATLAB. 9.9[END_REF].

The last element in the workflow is the optimizer, which uses algorithms to find the best set of decision parameters with respect to the formulation described in Section 3.1. Two types of optimization problems are consid- Having a surrogate model instead of CFD computations inside the analysis workflow higly reduces the cost of the function evaluations. Therefore, it allows the use of algorithms that require thousands of evaluations to reach convergence. Genetic algorithms are efficient when it comes to finding an optimum in an unknown complex design space [START_REF] Gen | Genetic algorithms and engineering optimization[END_REF]. For this study, the choice is made to use a Differential Evolution (DE) algorithm [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF], and more particularly its implementation proposed in the SciPy Python library [START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. This algorithm has the advantage of managing boundaries on decision parameters, as specified in Tab. 1. However, considering range maximization as a single-objective problem does not consider the induced flight time increase. In order to take it into account, a multi-objective optimization approach is introduced in a second time. Instead of looking for maximal range only, this method focuses on finding the best compromises between flight distance and time. The Non-Sorted Genetic Algorithm II (NSGA-II) [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] appears well suited in this case and is used as multi-objective optimizer. The implementation provided by the Platypus [START_REF]Platypus -multiobjective optimization in python[END_REF] Python library is chosen to perform the optimization. A summary of the algorithms and the associated libraries is presented in Tab. 3.

AERODYNAMIC CHARACTERIZA-TION AND SURROGATE MODEL

This section is dedicated to providing more explanation on the computation of aerodynamic coefficients and the use of surrogate modelling. Depending on the flight phase, the presence or absence of canards influences the behaviour of the projectile. The aerodynamic characterization as well as the trajectory computation method change accordingly.

Ballistic configuration

During the ballistic phase, the BALCO trajectory simulation uses a table of aerodynamic coefficients for the projectile without canards. Considering that this tool performs a Six-Degrees-Of-Freedom (6DoF) computation, the required aerodynamic coefficients consist in the axial, lateral and normal force coefficients denoted as (C A ,C Y ,C N ) as well as the associated moments computed at the center of gravity (C l ,C m ,C n ). The coefficient table

Figure 5: Sketch of the ballistic projectile geometry with its vertical and horizontal symmetry planes consists in their evaluation at various flight conditions, i.e. for discrete values of the Mach number M, the angle of attack α and the angle of sideslip β . Under some assumptions, the set of CFD simulations required to generate this table can be simplified.nd First, the ballistic projectile geometry comprises several symmetry planes, as illustrated in Fig. 5. Thanks to these symmetries, the aerodynamic behaviour can be only considered in the vertical plane, and the remaining coefficients necessary for 6DoF can be deduced from it. Similarly, the dependency on the sideslip angle is obtained from the angle of attack. This implies that only (C A ,C N ,C m )(M, α) have to be considered in the following. In addition, the symmetries imply that the normal force and pitching moment are null at zero angle of attack (Eqs. 1 and 2).

C N (α = 0°) = 0 (1) C m (α = 0°) = 0 (2)
Then, the projectile stability has been assessed previously in [START_REF] Libsig | CS 2020 -long range guided projectile (LRGP)[END_REF] and showed that the ammunition remains stable during the complete ballistic flight. Consequently, the angles of attack and of sideslip are supposed to vary in a limited interval around zero, and a small-angle approximation can be introduced. This hypothesis translates to a constant axial force coefficient in Eq. 3. Moreover, the lift an pitching moment coefficients are assumed to follow a linear evolution as written in Eqs. 4 and5.

C A (α) ≈ C A (α = 1°) (3) C N (α) ≈ C N (α = 0°) + α × ∂C N ∂ α (α = 0°) (4) C m (α) ≈ C m (α = 0°) + α × ∂C m ∂ α (α = 0°) (5)
By approximating the coefficient derivatives using their values at α = 0°and α = 1°, then using Eqs. 1 and 2, the final formulation is obtained in Eq. 6.

  C A C N C m   (α, M) ≈   C A (α = 1°) α ×C N (α = 1°) α ×C m (α = 1°)   (M) (6)
To summarize, under the preceding assumptions, only one value of the angle of attack is required at each Mach number to evaluate the aerodynamic coefficients of the projectile. Simulations are performed at α = 1°, and the total cost solely depends on the investigated Mach numbers, ranging from M = 0.3 to M = 3.0.

Finally, the decision parameters detailed in Tab.1 do not affect the projectile aerodynamics before canards opening. Therefore, the coefficient table is computed once and used in the ballistic flight phase for every trajectory simulation.

In order to perform the aerodynamic computations of the projectile, a mesh comprising 5 × 10 6 cells is generated. Reynolds Averaged Navier-Stokes (RANS) simulations are performed with the implicit density-based solver implemented in ANSYS Fluent. The kω Shear Stress Transport (SST) turbulence model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] is chosen, as well as second-order spatial discretization.

Each CFD simulation needs between 1 and 2 hours to converge on 224 cores of the supercomputer at ISL. Consequently, the computational cost associated to the ballistic projectile characterization is acceptable, with a duration of several days at most.

Gliding configuration

After canards deployment, the projectile trajectory is supposed to be contained in a vertical plane, allowing a 3DoF simulation. In this case, only 3 aerodynamic coefficients are required as inputs for the trajectory computation, namely (C A ,C N ,C m ). Moreover, the sweep angle of the deployed canards, A Can , has an impact on the aerodynamic behaviour of the projectile. Each time this design parameter is modified, the coefficient tables have to be updated accordingly. If the aerodynamic characterization is done using in-line CFD simulations in the optimization workflow, the set of computations has to be performed at each iteration.

The aerodynamics of the projectile with deployed canards not only depend on M and α, but also on the deflection angle of the canards δ . Assuming that N α , N M , N δ are the numbers of discrete values taken by each parameter for the characterization, usually around a dozen, then the total cost can be expressed as N α × N M × N δ ∼ 10 3 . Consequently, the full CFD characterization of a gliding projectile geometry represents a significantly increased cost in comparison with the ballistic configuration. This one-time computational cost can hardly be afforded on the supercomputer at ISL, and as this is to be repeated at each optimizer iteration, this solution is not feasible.

Superposition hypothesis

In order to mitigate computational cost and enable the characterization with CFD nevertheless, a superposition hypothesis is proposed. This assumption considers that the aerodynamic effect of canards deflection can be superposed with the effect of the angle of attack. The resulting equations are developed in the following.

First, the effect of canards deflection on the overall aerodynamic coefficients is assumed to be obtained from a set of simulations at zero angle of attack, as expressed in Eq. 7.

  δC A δC N δC m   (δ , M) =   C A C N C m   α=0°( δ , M) -   C A C N C m   α=0°( δ = 0, M) (7)
Then, this contribution is summed with coefficients from simulations at δ = 0°, in order to approximate the total aerodynamic coefficients for any flight condition. This leads to the expression given in Eq. 8.

  C A C N C m   (α, δ , M) ≈   C A C N C m   δ =0°( α, M) +   δC A δC N δC m   (δ , M) (8)
The superposition hypothesis allows in fact a decomposition of the coefficients dependency on α and δ . Thanks to this formulation, the aerodynamic characterization of a projectile geometry can be estimated with two sets of CFD simulations. Considering that each variable is discretized over a dozen of values for the characterization, the resulting total number of computations is

(N α × N M ) + (N δ × N M ) ∼ 10 2 .
Although the overall computational cost appears now acceptable for a onetime characterization, it remains too significant to be operated inside the optimization process. For that reason, the cost mitigation brought by the superposition hypothesis appears interesting but not sufficient.

Surrogate models

To further reduce the cost of aerodynamic characterization during the gliding flight, a surrogate modelling approach is introduced. The role of the models is to replace in the workflow the time-consuming sets of CFD simulations required to evaluate projectile aerodynamic coefficients. In particular, it should be able to predict the coefficient table depending on the value of the decision parameter A Can .

Under the superposition assumption, two CFD coefficient tables at α = 0°and δ = 0°are required to characterize the ammunition. Thus, two extensive databases of CFD simulations have been performed to evaluate (C A ,C N ,C m ) δ =0°( α, M, A Can ) and (C A ,C N ,C m ) α=0°( δ , M, A Can ). A mesh comprising 12 × 10 6 cells is generated for the gliding geometry, and the CFD parameters detailed in Section 4.1 are re-used. Each database comprises hundreds of CFD simulations and the overall computational time was approximately of two weeks on the supercomputer at ISL.

The computations are post-processed according to Eqs. 7 and 8 to generate training datasets. An approach using Radial-Basis Functions (RBF) interpolation method [START_REF] Martin D Buhmann | Radial basis functions: theory and implementations[END_REF] is chosen for the surrogate models. More precisely, the RBF process is operated using Gaussian basis functions, as provided in the Surrogate Modelling Toolbox (SMT) [START_REF] Bouhlel | A python surrogate modeling framework with derivatives[END_REF] Python framework. The two surrogate models are then trained on the training databases to generate

(C A ,C N ,C m ) α=0°( α, M) and (δC A , δC N , δC m )(δ , M) for a given value of A Can .
Finally, the quality of the surrogate models is assessed before using the models for optimization. They appear to have a satisfying accuracy over the complete range of flight conditions required for trajectory simulation and parameters for optimization.

SINGLE-OBJECTIVE OPTIMIZA-TION

The first optimization is conducted with a singleobjective formulation. Therefore, the genetic algorithm aims at maximizing the range of the projectile. The differential evolution process explores the design space and stops after 47 generations of 40 samples, totalizing 1880 evaluations.

A representation of the function evaluations can be observed in Fig. 6. In particular, it shows how the first generations of samples are widespread over the design space, while the last generations are evolving towards sets of parameters resulting in maximal range values. The algorithm evolution, in terms of the functions of interest, is presented in Fig. 7. The global trend of the graph depicts the strong correlation between the distance attained and the time of flight. Despite this expected dependency, it can also be observed that the decision variables allow to reach a given distance in a significant time interval, leaving room for future optimizations. Between the first and last sample generations, the mean range has been improved by 57%. As a consequence, the time of flight is increased by 52%. The optimal set of decision parameters is summarized in Tab. [START_REF] Arnoult | Control surface geometry surrogate-based optimization for spin-stabilized projectile course correction[END_REF]. The optimizer has found the best 39.9 s compromise between canards-induced lift and drag over the entire gliding path by setting a low value of A Can . The maximum angle of attack is limited at 10, 7°, probably in order to avoid the interaction of the canards wake with the tail-fins. A significant elevation angle is retained, in order to provide a sufficient altitude before starting the glide. Moreover, the canards opening happens shortly after the apogee of the ballistic flight in order to take advantage of this altitude.

In conclusion, this first application of optimization for the design of long-range guided projectiles achieves a significant improvement of the flight range. However, a single-optimization formulation, only focusing on maximizing range without time constraint, is not representative of operational reality. Thus, the next step consists in considering multi-objective optimization problems.

MULTI-OBJECTIVE OPTIMIZA-TION

Multi-objective optimizations consider multiple functions of interest, which can be of great benefit for the long-range guided ammunition design. Such approach is implemented in the following, with a research of the best The research process is launched over 100 generations of 50 samples, therefore leading to a total of 5000 analysis evaluations. All evaluated points can be observed on the visualization matrix displayed in Fig. 8. The evolution of range and time of flight over the successive samples generations is plotted in Fig. 9, as well as the resulting Pareto front. The successive generations appear to be progressively converging towards the final front. On the Pareto border, the evolution of optimal flight time with respect to the range appears to be non-linear. Due to the projectile deceleration during the gliding phase, increasing the range at long distance is more time-demanding than at shorter distances. / In practice, the Pareto front can be of great interest for operational use and tactical decision making. It indicates the maximal reachable distance for a given flight time, or on the contrary the trajectory of minimal flight time when considering an objective at a certain range.

Finally, a deeper analysis of the set of points belonging to the Pareto front is proposed. The evolution of the decision parameters along the front is plotted in Fig. 10. Some tendencies appear clearly on this plot and give a better understanding of the points constituting the frontier. First, it seems that increasing the projectile angle of attack and the gun elevation angle are mandatory conditions to reach long-distance targets. Then, the optimal value of canards sweep angle is barely modified along the optimal line; thus, this parameter could be fixed in future studies. In addition, the optimal moments for canards deployment seem to have a less critical influence on

CONCLUSIONS

In this paper, an alternative is presented to usual technical devices extending the fire range of howitzers. The proposed solution consists in a guided ammunition concept with tail-fins and unfoldable canards. Accordingly, the flight scenario of such projectile is decomposed in two phases, with a ballistic phase followed by canards deployment and a steered gliding phase.

Then, an optimization approach is introduced, in order to find the most-suited design and flight parameters. The workflow developed for this purpose comprises an aero- To conclude, the proposed solution has the advantage of using additional lifting surfaces and guidance, instead of requiring base-bleed or propulsion device. Hence, this long-range guided concept is expected to come at reduced cost. As a perspective, a hybrid ammunition configuration combining propulsion, additional lifting surfaces and guidance could be considered to reach very-long-range targets.

1

 1

Figure 1 :Figure 2 :

 12 Figure 1: Top-views schematics of the projectile before (a) and after canards deployment (b)

Figure 3 :

 3 Figure 3: Definition of the canards design parameter

Figure 4 :

 4 Figure 4: Optimization workflow for long range ammunition

Figure 6 :

 6 Figure 6: Evaluations performed by the genetic algorithm (coloured by flight range)

Figure 7 :

 7 Figure 7: Evolution of flight range and time throughout genetic algorithm iterations

Figure 8 :Figure 9 :

 89 Figure 8: Evaluations performed by the NSGA-II algorithm (coloured by flight range)

Figure 10 :

 10 Figure 10: Evolution of normalized decision parameters along the Pareto front, depending on the associated flight time

Table 1 :

 1 Recapitulative table of optimization parameters

	Decision parameters Lower -upper boundaries
	A Can	30,0 °-90,0	°Pro
	jectile AOA max	0,0 °-16,0	°Gun
	elevation angle	30,0 °-80,0	°Opening
	time canards	20,0 s -80,0 s

function is to be maximized by the optimization process.

Table 2 :

 2 Recapitulative table of functions of interest

	Functions of interest	Objective
	L f light	Maximization
	t f light	Minimization

Table 3 :

 3 Recapitulative table of optimization algorithms

	Formulation	Algorithm	Python library
	Single objective	Differential evolution	SciPy.optimize
	Multi-objective	NSGA II	Platypus
	ered in this work. First, a single-objective optimization
	is implemented, in order to maximize the flight range.

Table 4 :

 4 Resulting optimal configuration, obtained with single-objective formulation

	Decision parameter Optimal value
	A Can	36.2	°Pro
	jectile AOA max	10.7	°Gun
	elevation angle	61.8	°Opening
	time canards


David Monniaux

Completeness in Static Analysis by Abstract Interpretation, a Personal Point of View

Static analysis by abstract interpretation is generally designed to be "sound", that is, it should not claim to establish properties that do not hold-in other words, not provide "false negatives" about possible bugs. A rarer requirement is that it should be "complete", meaning that it should be able to infer certain properties if they hold. This paper describes a number of practical issues and questions related to completeness that I have come across over the years.

Introduction

The concept of completeness, under several definitions, permeates mathematical logic. A proof system is deemed complete with respect to a semantics if it can prove all properties that are true in this semantics. For instance, Gödel's completeness theorem states that first-order logic (that is, any reasonable proof system for it) can prove any property true in all models. A decision procedure for a logic (that is, a procedure that answers "true" or "false" for any formula in the logic F ) is expected to be sound (it does not declare to be true formulas that are not true) and complete (it declares to be true all formulas that are true).

The concepts of completeness (a) for a proof system with respect to a class of properties (b) for a procedure that searches for proofs of such properties within such a proof system are related, but distinct. Obviously, if a proof system is incomplete, then so is any procedure that searches for proofs within that system. However, it is possible to have a complete proof system and an incomplete procedure for searching within it, for instance for practical efficiency reasons.

In abstract interpretation, completeness, at least at the global level (proving properties of whole programs), is often forgone straight from the start when designing the abstraction (for instance, interval arithmetic is used even though it is obvious that it cannot prove safety in general, even if the property to prove is itself an interval), and there thus may be little reluctance to adding further incompleteness in the solving procedure by using approaches such as widening operators [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF][START_REF] Cousot | Abstract Interpretation Frameworks[END_REF]. Yet, it is interesting to control this further incompleteness, both from a theoretical and a practical points of view. Completeness in the solving method ensures some predictability in the analysis outcome, while the brittleness of some incomplete approaches (the approach succeeds or fails depending on seemingly inconsequential aspects of the program and the property) surprises end users: for instance, providing more precise information on the precondition of a program may prevent the analysis procedure from proving a property that it could prove with a less precise precondition.

In this paper, I discuss various instances of completeness and incompleteness that I have came across, without any pretense of thorough theoretical discussion and, ironically, no pretense of completeness. For more theoretical discussion, Giacobazzi et al. [START_REF] Giacobazzi | Making abstract interpretations complete[END_REF]Sec. 7] surveyed the literature extensively. For a study of completeness with respect to the relationship between the semantics and the abstract domain, and constructions to make the domain complete, see [START_REF] Giacobazzi | Completeness in Abstract Interpretation: A Domain Perspective[END_REF].

Completeness of the Abstraction: the Case of LRU Caches

It is well-known that it is equivalent to (a) evaluate an integer multivariate polynomial over integer inputs, then take the final result modulo N (b) do the same computation but reducing modulo N at every step. 1 Due to a strong mathematical property (a ring morphism), it is possible to replace an expensive computation, possibly involving large numbers, by a simpler one that abstracts the expensive one by operating on small abstractions of the data.

This example is ideal, and one could doubt the existence of complete abstractions that significantly simplify computations for nontrivial program analysis problems; yet we came across one such example, where we perform an exact analysis defined by a composition of exact abstractions and an exact solving procedure.

All current processors use some form of cache memory: frequently accessed code or data is retained in fast memory close to the processor core, avoiding much slower accesses to main memory. Cache analysis consists in determining, given a program and a description of a processor cache subsystem, which memory accesses are always "cache hits" (data already in cache), which are always "cache misses" (data not in cache, thus access to external memory), or more complex properties ("the first time this statement is executed, the access is a cache miss, then subsequent accesses are hits"). Such an analysis is for instance used as a prerequisite for worst-case execution time analysis.

A possible approach to solving the cache analysis problem would be to apply a model-checker to the combination of the program under analysis (or some suitably simplified model thereof) and of the cache system. A first and obvious abstraction is, for the cache system, to retain only the addresses of the memory blocks currently stored within the system, but not their contents. This abstraction is exact: the cache system finds data, or decides to evict data to make room for new contents, according to addresses only, never according to contents. The resulting model of the cache system, for a fixed processor model, is finite: the cache consists of 1 For N = 9, this is the basis for the method of "casting out nines", named so because reduction modulo N can be implemented over numbers in decimal form by summing the digits of the number, with 9 being replaced by 0 (repeat the process until the result consists on one single digit). Schoolchildren used to apply this method to check their hand computations: the final result of an error-prone integer computation, taken modulo 9, should be identical to the result obtained with reduction at every step. This is of course not a sound check: if results differ, one is sure that there has been some error somewhere (no false positives), but they can coincide by chance.

a fixed number of cache lines, buffers, etc., labeled with addresses taken within a fixed range. However, such a model is in practice intractable on all but the simplest examples.

The naive vision of cache memory is that of a fully associative cache: any cache line may be used to store any block from the main memory, regardless of its address. In reality, a cache is usually split into several cache sets, each of which able to retain blocks present only at a given class of addresses in the main memory. The classes of addresses for the different cache sets form a partition of the possible addresses in main memory. In almost all cases2 these cache sets operate completely separately from each other. It is thus possible to perform cache analysis separately for each cache set; this is again an exact abstraction. However, this still results in intractable models.

Each cache set is composed of a fixed number of cache lines; the number N of lines is also known as the number of ways or the associativity of the cache. Each cache line may be empty or nonempty; a nonempty cache line stores a block from main memory: its address and its contents (we have already seen that the contents are not relevant for analysis). For simplicity, we shall consider here the case of single-level, read-only caches. A read from a memory location not currently in the cache triggers a load from main memory into the cache. Except in the rare case when the relevant cache set still has empty lines, this involves evicting a block from the cache set. The selection of the block to be evicted is made according to a cache replacement policy. The most obvious replacement policy is to evict the least recently used (LRU) block. 3 A cache set is therefore modeled as a sequence of at most N block addresses, in increasing age: the youngest block is the one accessed most recently, the oldest is the least recently accessed. The positions of blocks within the set change as data is accessed. Let us take an example with N = 4. Assume a, b, c, . . . denote distinct addresses in memory, and assume the cache set initially contains abcd, meaning that a is the youngest block and d the oldest. If the processor requires block d, then it is rejuvenated and the cache set then contains dabc. If then the processor requires block e, then the oldest block (c) is evicted and the cache set contains edab.

Let us now see our approach to analyzing programs over LRU caches [START_REF] Touzeau | Fast and Exact Analysis for LRU Caches[END_REF]. Consider the case where we want to know whether a certain fixed block a is within the cache at certain program locations. Then, the blocks older than a in the cache are unimportant; so is the ordering of the blocks younger than a. The idea is that if a block is younger than a, then accessing it again will just change the ordering of blocks younger than a, but not the age of a; and if it is older, then accessing it will increment the age of a.

A cache set may thus be described, with respect to a block a, by an abstract configuration: either "a is absent" or "a is here and the set of blocks younger than a is {. . . }". Again, this is an exact abstraction.

Cache analysis using that abstraction specialized for accesses to a thus amounts to collecting the sets (for all program locations) of reachable abstract configurations in the cache set where a is to be stored. This is still very costly.

A collection of abstract configurations for that cache set consists of a Boolean "is it possible that a is absent" and a set of sets of size less than N of blocks distinct from a. A crucial observation is that, in this set of sets, only the minimal and maximal elements (with respect to the inclusion ordering) matter: if it is possible that a is preceded by {b}, {b, c}, {b, c, d}, then, for checking if it is possible that a (at some later point) is out of the cache, it is only necessary to retain {b, c, d}, which subsumes the other two cases; and for checking if it is possible that a (at some later point) is in the cache, it is only necessary to retain {b}, which again subsumes the other two cases. The reason for this subsumption is that if a sequence of steps leads from a state where a is preceded in the cache by a set P of blocks to the next access to a, and that access is a miss, then the same sequence of steps of steps, starting with any cache state where the set of blocks younger than a is a superset of P , also leads to a miss (mutatis mutandis for subsets and hits). Retaining only the extremal elements is thus, again, an exact abstraction from the point of view of cache analysis.

To summarize, we first gradually simplified the cache state by abstracting away parts and aspects that are irrelevant to the analysis under way, then we simplified the collection of abstract cache states by using a subsumption relation. The resulting abstraction is exact with respect to the properties under analysis.

A this point, what remains is how to implement this abstraction. What is needed is an efficient way to represent sets of incomparable sets of blocks (antichains); we opted for zero-suppressed binary decision diagrams (ZDD). For better efficiency, this analysis, still somewhat costly, is applied to memory accesses only after some cheap approximate analyses [START_REF] Touzeau | Ascertaining Uncertainty for Efficient Exact Cache Analysis[END_REF] have failed to classify these accesses.

An advantage of an exact approach is that, since the result is uniquely defined, it is possible to test implementations of various exact abstractions against each other; the results should be identical. This way we discovered a subtle ZDD bug by comparison with a (less scalable) model-checking approach.

Granted, this approach is "exact" or "complete" only because a simplified model of the program (a control-flow graph, not taking into account the data being processed) is used. This model introduces approximation: for instance, in the following program if ( flag ) access ( a ); else access ( b ); if ( flag ) access ( a ); else access ( b ); the only feasible sequences of accesses are a, a and b, b, (and thus in either case the second access is a hit), but an analysis based on control-flow only will not track the value of the flag and consider that sequences a, b and b, a are also feasible, and conclude that the second access may be a miss. However, tracking the value of data precisely makes the model undecidable. An intermediate solution that retains decidability is to consider control-flow with precise modeling of procedure calls [START_REF] Monniaux | The complexity gap in the static analysis of cache accesses grows if procedure calls are added[END_REF].

On a final note, on many examples [START_REF] Touzeau | Fast and Exact Analysis for LRU Caches[END_REF], worst-case execution time (WCET) analysis ended up being overall faster with the help of this apparently expensive analysis. The reason is that the WCET analysis in use involved enumerating all reachable pipeline states. If the cache analysis is imprecise and cannot conclude about some access, where the exact analysis would conclude to a "always miss", this does not normally change the bound on WCET, but this increases the number of pipeline states to consider. In short:

• incomplete analyses are first applied in order to answer most subproblems exactly; only the subproblems for which the answer is definitely unknown are passed to the more expensive complete analysis;

• completeness is ensured by a composition of abstractions that do not lose any precision with respect to the properties we are interested in, but which greatly simplify the model to analyze;

• the resulting model is analyzed exactly;

• complete analyses have an exactly defined result, thus testing for bugs is easy by comparing results;

• for the same reason, it is possible to study the complexity of the problem [START_REF] Monniaux | On the Complexity of Cache Analysis for Different Replacement Policies[END_REF][START_REF] Monniaux | The complexity gap in the static analysis of cache accesses grows if procedure calls are added[END_REF];

• the extra cost of the complete analysis may be offset by savings in client analyses, because a more precise result means fewer case analyses down the road.

Completeness or Incompleteness of the Analysis Method

A classical question of completeness in abstract interpretation is whether the abstract domain in use is sufficient to prove the properties being sought. Yet, even if the abstract domain is sufficient, it may happen that the abstract interpretation method being used cannot find the necessary invariants within the domain.

Widening Operators

Let us see an example: some simplified version of a dataflow program where i is an index into some circular array of length 42.

i = 0; while ( true ) { if ( trigger ()) { i = i +1; if ( i > 42) i = 0; } assert ( i < 1000); }
Textbook forward static analysis by intervals [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF][START_REF] Cousot | Méthodes itératives de construction et d'approximation de points fixes d'opérateurs monotones sur un treillis, analyse sémantique de programmes[END_REF] will compute, for variable i, a sequence of intervals [0, 0], [0, 1], [0, 2], . . . , and widen it to [0, ∞). If narrowing iterations are used, in their simplest form by starting from the [0, ∞) invariant, running the loop body once more and adding the initialization states, one obtains [0, 999]. The assertion cannot be proved using that invariant. Yet the assertion would have been proved if the analysis had guessed the least possible interval, [0, 42]. Clearly, the problem is with the inference method (iterations with widening), not the abstract domain (intervals). Iterations with widening are an incomplete method with respect to the abstract domain: they may fail to discover suitable inductive invariants when some exist in the domain.

Note also a surprising characteristic of this incomplete approach: if instead of the precondition i = 0 we had analyzed the loop with the precondition 0 ≤ i ≤ 42, we would be able to prove the assertion correct. This non monotonic behavior (a less precise precondition, or a coarser abstraction of some program semantics leads to more precise analysis results) may surprise end users.

Another surprising characteristic of widening operators for relational domains is that they do not commute, in general, with projection, meaning that adding some extra variables may change the result of the analysis on the original variables even if the extra variables have no influence on the original ones (for instance, if these variables are mere "observers"). Consider for instance the program:

i =0; j =0; while ( true ) { if (*) i =0; else { i =1; j =0; } if ( i ==0) j = j +1; }
In this program, j observes the number of last iterations during which i was 0. The set of reachable states of this program for (i, j) is (1, 0) ∪ {(0, n)|n ∈ N}; in particular the set of reachable values for i is {0, 1}, and it can be computed exactly with interval analysis by postponing widening by one step, or by applying one step of narrowing. Now consider the sequence of polyhedra, that is, polygons, computed by polyhedral analysis with widening [START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF]: at step n, the polygon is defined by vertices (0, 0), (1, 0), (0, n). The lines passing through vertices (0, 0) and (1, 0) (constraint j ≥ 0) and the one passing through vertices (0, 0) and (0, n) (constraint i ≥ 0) are stable; the line passing through vertices (1, 0) and (0, n) is unstable and will be discarded by widening. The resulting system of constraints is i ≥ 0 ∧ j ≥ 0; we have lost the i ≤ 1 constraint obtained by interval analysis. In some cases, a stratified approach, where successive analyses take increasing subsets of variables in consideration, may be able to recoup precision [START_REF] Monniaux | Stratified Static Analysis Based on Variable Dependencies[END_REF].

Exact Solving

A safety property is established using forward abstract interpretation by inferring some inductive invariants in the abstract domain and checking that these invariants imply the property. A complete abstract interpretation method is one that would compute always such invariants if they exist.

Clearly, any method that computes the least inductive invariants in the abstract domain is complete: if there exist invariants in the domain that can prove the safety property, then, a fortiori, the least inductive invariants in the domain can prove that property.

One such method is ascending iterations without widening within a domain that satisfies the ascending chain condition (no infinite strictly ascending sequences), which ensures termination. This includes domains of finite height (there is a bound on the length of strictly ascending chains), and in particular finite domains, such as powersets of finite sets. If f : L → L is a monotone function, and ⊥ is the infimum of L, then the sequence f n (⊥) is ascending and, in all the above cases, becomes stationary. When f n (⊥) = f n+1 (⊥), then f n (⊥) is the least solution of x = f (x). The algorithms used in practice are algorithmic improvements over this idea. 4 An example of an infinite domain of finite height is that of solution sets of systems of linear equations over a fixed number of variables [START_REF] Karr | Affine Relationships Among Variables of a Program[END_REF]. When a solution set S is strictly included in a solution set S ′ , the dimension of S ′ is strictly greater than that of S; this dimension cannot be more than n, which thus bounds the height of the lattice.

The domain of intervals has none of these characteristics; but it is however possible to compute the least inductive invariant within that domain of programs such as the one above, by specifying this least inductive invariant as a least fixed point, writing a system of numeric equations that this fixed point should satisfy, and solving this system for the least solution. Indeed, for the above example, let us write down the equations that h should satisfy for [0, h] to be a fixed point for the loop:

h = min(max(min(42, h + 1), h), 999) (1)
Let us solve this equation. In any solution, the outermost "minimum" operator must be equal to one of its arguments. Assume it is the second argument, 999. One can indeed verify that 999 is a solution to this equation. Now consider the case where it evaluates to its left argument. The equation then becomes h = max(min(42, h + 1), h). Similarly, the "maximum" operator evaluates to one of its operands. Assume it evaluates to its right argument. The system then simplifies to h = h; but this is under the condition that h ≤ 999 and h ≥ min(42, h + 1). When h + 1 > 42, that is, h ≥ 42, this minimum is equal to 42; thus all 42 ≤ h ≤ 999 yield valid solutions. When h + 1 ≤ 42, this minimum is equal to h + 1; but then the condition becomes h ≥ h + 1, which would imply h = ∞, which contradicts h ≤ 999. Now assume the "maximum" operator evaluates to its left argument. The system then simplifies to h = min(42, h + 1), which yields only h = 42 as solution. We conclude that h = 42 is the least solution.

The above reasoning by case analysis over the minimum and maximum operators can be automated through exhaustive case analysis or, most subtly, by a SMT (satisfiability modulo theory) solver. Many such solvers can also optimize within the solution space, and thus directly produce the least fixpoint. Alternatively, for proving safety properties, it is sufficient to query the solver for any solution sufficient for proving the safety property.

A more principled approach to the case analysis for the maximum operators is ascending policy iteration [START_REF] Gawlitza | Precise Fixpoint Computation Through Strategy Iteration[END_REF], which considers an ascending sequence of fixpoints of systems obtained by picking arguments to the "maximum" operators. All the above approaches generalize to domains other than intervals, for instance to zones [START_REF] Martin | Precise Program Analysis through Strategy Iteration and Optimization[END_REF], at the expense of some algorithmic complications.

Let us go further. Assume the abstract domain consists of the sets defined by a parameterized predicate I(p, x) where p is the vector of parameters, x is the program state. Note that this encompasses intervals, octagons, template polyhedra domains, or even polyhedra with a fixed number of faces. Also note that I may contain disjunctions and express non-convex relations. 4 In many cases, L = L m b where L b is a base lattice. An element x of L is then decomposed into x 1 , . . . , xm, and f is decomposed into its m projections f 1 , . . . , fm. The problem is then to find a solution of a system of equations x 1 = f 1 (x 1 , . . . , xm), . . . , xm = fm(x 1 , . . . , xm), typically using a system of working set of variables being updated, with some judiciously chosen ordering for choosing the next update to process.

The inductiveness condition for a transition τi,j from program location i to program location j may thus be written as the Horn clause ∀x, x ′ I(pi, x) ∧ τi,j(x, x ′ ) ⇒ I(pj, x)

Safety conditions may be expressed as ∀x Ii(pi, x) ⇒ Ci(x), and program startup may be specified as ∀x Si(x) ⇒ Ii(pi, x). Inductive invariants within the abstract domaine suitable for proving the safety properties are thus expressed as the solutions, in the parameters pi, of a system of Horn clauses.

If I, the transitions τi,j , the safety conditions Ci and the start conditions Si are all expressed within a theory admitting algorithmic quantifier elimination, such as Presburger arithmetic or the theory of real closed fields, then the existence of such invariants is decided by quantifier elimination. For instance, if we consider a program operating over real numbers, and invariants of the form A(p)x ≤ B(p) where A(p) and B(p) have a fixed number r of rows (a polyhedron with at most r faces), then the existence of a suitable invariant is established by quantifier elimination in the theory of real closed fields. It is even possible to specify that one wants the least inductive invariant and obtain it by quantifier elimination [START_REF] Monniaux | Automatic Modular Abstractions for Linear Constraints[END_REF].

Imprecise Abstract Transfer Functions

Widening operators are the best known source of non-monotonicity and incompleteness in finding suitable invariants in abstract interpretation, but they are not the only one. There are, and this may be surprising, cases of non-monotonicity and incompleteness due to transfer functions, particularly in combinations of abstractions.

It is well-known that even if each individual step of abstract interpretation computes the least possible post-condition within the abstract domain compatible with the precondition, this property does not extend to composition of steps; in other words, optimality does not compose. Let us see a simple example:

y = x ; z =x -y ;
If x is known to lie in [0, 1], then so does y, and this is optimalnon-optimal but sound answers would be intervals containing [0, 1]. Then interval arithmetic deduces that z is in [-1, 1], which is sound, but not optimal: the optimal interval is [0, 0]. However, in order to reach that conclusion, one must know that x = y at the intermediate control point, that is, a relation between x and y, which is not possible in a non-relational abstraction such as intervals.

One workaround to this weakness is to propagate a system of relations along with the interval analysis. For instance, we we had propagated y = x, then by substituting y → x into xy and simplifying, we would have obtained z = 0.

An approach to implement this workaround [START_REF] Cousot | Combination of Abstractions in the ASTR ÉE Static Analyzer[END_REF], for instance applied in the Astrée static analyzer, is to propagate a terminating rewriting system (for instance, populated in chronological order) along with the interval information, perform interval analysis on both the original expressions and the expressions obtained by rewriting and simplification, and then take the intersection. Here, the rewriting system would contain y → x, and the rewritten and simplified expression would yield z = 0. This approach however creates non-monotonicity. 5 Consider this example [8, ex. 4]:

Code

Symbolic computation Less precise symbolic int i, j=i+1;

j → i + 1 nothing int k=j+1; j → i + 1, k → j + 1 → i + 2 k → j + 1 if (j > 0) l=k; j → i + 1, k → i + 2, l → i + 2 k → j + 1, l → j + 1
Assume we don't know anything about the values of the variables initially, and consider the interval obtained for l at the end of the "then" branch of the test. Simple interval propagation yields no information about l. The rewriting system yields l → i + 2 and since no information is known about i, no information is known about l. Yet, if we forget that j → i + 1, and apply the resulting rewriting system, we obtain l → j + 1 and thus l > 1.

Granted, some improvements could be made by more clever propagation. 6 For instance, since j → i + 1, the guard could be rewritten into i + 1 > 0, and this could itself be rewritten into i > -1. This amounts to "inverting" j → i + 1 to i → j -1. More generally, one would need to consider the rewriting system not as a directed system, but as a system of equations. The combination of a system of linear equations and intervals forms the basis of the simplex algorithm for linear programming, and obtaining optimal interval bounds from a combination of equalities and interval bounds amounts to linear programming. In fact, any convex polyhedron can be expressed as the projection of a set defined by the conjunction of linear equalities and interval bounds. 7 The appropriate domain for dealing with such constraints precisely would thus be that of convex polyhedra [START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF][START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF].

Undecidability of an Abstraction

The argument from Section 3.2 for establishing decidability of the existence of an inductive invariant within the domain of convex polyhedra with at most r faces does not extend to the domain convex polyhedra with any number of faces. In fact, there is no known algorithm for deciding the existence of inductive invariants within that domain for any nontrivial class of programs. Invariant inference approaches for that domain, starting from the early proposals from Cousot and Halbwachs [START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF][START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF], are typically based on iterations with widening. An intriguing question is whether it is inevitable to resort to heuristics.

Polyhedral Abstraction

I have attempted proving that there is no algorithm deciding the existence of polyhedral invariants for linear transition systems (real or integer), and have so far failed. Because of my efforts in raising attention to Figure 1: When the process is non terminating, any polyhedron containing the reachable points inevitably goes below the parabola this issue, some colleagues named this the "Monniaux problem". Several distinguished researchers 8 have also reported trying to solve it and failing.

Undecidability with Quadratic Guards

Progress has however been made on variants of the so-called "Monniaux problem" [START_REF] Fijalkow | On the Monniaux Problem in Abstract Interpretation[END_REF]. I was able to prove that this problem becomes undecidable if one is allowed to use a quadratic transition guard: it is then possible to encode a deterministic counter machine reachability problem into the invariant inference problem. Let us recall how (proof sketch) [START_REF] Monniaux | On the decidability of the existence of polyhedral invariants in transition systems[END_REF].

Let the counter machine M operate over variables z1, . . . , zn, to which we add two variables x and y, initialized to 0. The transitions of the counter machine are synchronously combined with steps (x, y) → (x + 1, y + x). The combined machine thus simulates the counter machine on variables zi together with a parabola n, 1 2 n(n -1) on variables (x, y), where n is the number of steps taken so far. Now modify the resulting machine by conjoining to all transitions the guard y = 1 2 x(x -1); clearly this does not modify the behavior of the machine. Add a special control state σ b , meant to be unreachable, and transitions from any state y < 1 2 x(x -1) to that state. Call the resulting machine M ′ .

Assume M terminates in N steps, then so does M ′ . Consider the convex hull H of the finite family of points such that x = k, y = 1 2 k(k -1), z1, . . . , zn is the state reached after k steps in the execution of M , and 0 ≤ k ≤ N . The only points in H that lie on the y = 1 2 x(x -1) parabola project to the reachable states of M . H, in general, contains extra states (integer points) such that y > 1 2 x(x -1), but due to the nonlinear guards, these states cannot transition to other states. Thus H is an inductive invariant for M ′ that proves the inaccessibility of σ b .

Assume M does not terminate, and thus so does M ′ . Any inductive invariant for M ′ must contain the infinite family of points such that x = k, y = 1 2 k(k -1), z1, . . . , zn is the state reached after k steps, for k ≥ 0. Any convex polyhedron that contains this infinite family of points must have a point y < 1 2 x(x -1): assume there is none and consider the vertex V * with maximal x * ; this vertex must lie on the parabola; there is for 8 Including some at the 2022 CSV workshop in Venice, whence this volume comes.

x ≥ x * at least one infinite face of the form y ≥ L(x, z1, . . . , zn); but then inevitably there will be points of that face lying below the parabola (Fig. 1). However, any convex polyhedron that contains a point where y < 1 2 x(x -1) is not an inductive invariant capable of proving that σ b is unreachable. It follows that there is no convex polyhedron that is an inductive invariant for M ′ capable of showing the unreachability of σ b .

Thus, M ′ has an inductive polyhedral invariant suitable for proving the unreachability of σ b if and only if M terminates, and this for arbitrary M in a class of machines with undecidable termination.

The same reasoning applies whether the state space is considered over integers or over the reals.

Perspectives and Dead Ends

The role of the parabola y = 1 2 x(x -1) in the above proof could be played by other sets, such as a circle x 2 + y 2 = 1: the process that enumerated points on the parabola can be replaced by one that enumerates points on the circle. Initialize (x, y) = (1, 0), and, as next step, apply a rotation matrix defined by a Pythagorean triple ((a, b, c) integers such that a 2 + b 2 = c 2 , for instance [START_REF] Cook | Soundness and Completeness of an Axiom System for Program Verification[END_REF][START_REF] Cousot | Méthodes itératives de construction et d'approximation de points fixes d'opérateurs monotones sur un treillis, analyse sémantique de programmes[END_REF][START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF]):

a c b c b c -a c
. This however leads to more complex proof arguments, and does not gain anything: we still need a nonlinear guard.

Essentially, what we used in both cases is a strictly convex set S (y ≥ 1 2 x(x -1)) with boundary definable by a guard (y = 1 2 x(x -1)), with exterior (y < 1 2 x(x-1)) definable by a guard, and with a way to enumerate points in the boundary while still remaining in the class of deterministic transitions under consideration. Strict convexity is used so that taking the convex hull of points in S, as when using polyhedral invariants, does not add "parasitic" points on the boundary. This ensures that the guard that keeps only points on the boundary removes "parasitic" points.

Could we achieve similar goals using only linear constraints? Suppose we can express such a guard with a formula, possibly containing disjunctions. Then, by distributivity, this formula expresses a finite union of polyhedra P1 ∪ • • • ∪ Pn. Then there is a i such that the enumeration process eventually picks two points (in fact, an infinity of them) in Pi; and thus the guard cannot exclude "parasitic" points between these two, at least over the reals.

A possible course of action could be to set the problem over integers and take advantage of the implicit guard that non-integer points are discarded. This however was also a dead end.

Richer Domains

If the abstract domain (class of invariants) under consideration is rich enough to express exactly the set of reachable states of the programs to be analyzed, then obviously the problem is undecidable: a suitable inductive invariant (the set of reachable states) exists in the domain if and only if the safety property is true.

This happens for instance if the domain includes Σ 0 1 formulas of Peano arithmetic, that is, formulas of the form ∃v1 . . . ∃vn F where F only contains bounded quantifiers and the usual arithmetic operations and comparisons. By a form of Gödel's encoding, for any integer program with m variables built using the usual arithmetic operators, it is possible to build a Σ 0 1 formula F (k, s1, . . . , sm) satified if and only if (s1, . . . , sm) is the state of the program after k steps of [START_REF] Winskel | The formal semantics of programming languages -an introduction[END_REF]Ch. 7][3].

Such domains are obviously too rich. In fact, with the domain of Σ 0 1 formulas described above, it is not even possible in general to check for inductiveness, because the formulas belong to an indecidable class. However, even semilinear sets (sets defined by Boolean combinations of linear inequalities with algebraic coefficients, thus disjunctions of convex polyhedra) are sufficient to obtain undecidability even for a very restricted class of programs (one single loop control state, nondeterministic choice between two linear transformations, no guards) [10, Th. 9].

Perspectives and Conclusion

The question of the completeness of the domain with respect to the properties to prove, and the class of programs under consideration, is distinct from the problem of algorithmically finding suitable invariants within that domain-or the related problem of computing the least inductive invariant within the domain.

The traditional approach for abstract interpretation in domains with infinite ascending chains is iterations with widening. This approach has known weaknesses: the analysis may miss the invariants necessary for the proof and the analysis behavior is non-monotone (increasing knowledge about preconditions or transitions may decrease the precision of the analysis), which is a source of "brittleness" (a minor change in the program causes the analysis to fail to prove the property). Many tricks are thus used to improve the invariants obtained with widening: narrowing iterations [START_REF] Cousot | Abstract Interpretation Frameworks[END_REF][START_REF] Cousot | Méthodes itératives de construction et d'approximation de points fixes d'opérateurs monotones sur un treillis, analyse sémantique de programmes[END_REF], lookahead widening [START_REF] Gopan | Lookahead Widening[END_REF], guided static analysis [START_REF] Gopan | Guided Static Analysis[END_REF]. Though they help in practice, none of them guarantees that analysis will succeed.

Over time, completely different approaches were designed to compute suitable arguments for proving safety properties. While widening is a form of extrapolation (look at sets of reachable states in 0, 1, 2 . . . steps and try to extrapolate for arbitrary number of steps), these methods are based on interpolation: given some under-approximation of the set of reachable states, and the property to prove, find a simple separation between the two (a Craig interpolant) and hope that the interpolant, or components thereof, or formulas constructed from components of successive interpolants, becomes inductive. One such approach is propertydirected reachability [START_REF] Bradley | SAT-Based Model Checking without Unrolling[END_REF][START_REF] Eén | Efficient implementation of property directed reachability[END_REF], later enriched with a number of extensions and improvements (dealing with arithmetic theories, dealing with non-linear Horn clauses and not just transition systems). Variants of this approach are in particular implemented in the popular Z3 solver. 9 Unfortunately, this approach, too, is brittle [START_REF] Braine | Verifying Programs with Arrays and Lists[END_REF]: minor differences in the problem to be solved (names of variables, etc.) may result in dramatic changes in outcome (finding invariants or timing out).

Approaches that are guaranteed to find the least inductive invariant in the chosen abstract domain (or, more generally, an inductive invariant suitable for proving a given property, if it exists) are more resilient; I have listed several of them in Section 3.2. Some of these methods however do not scale up well (those based on quantifier elimination); policy iteration seems to be the most scalable.

In our view, the case for continued use of widening operators, despite their known weaknesses, or other non monotonic and/or brittle meth-ods, would be strengthened by a proof that the existence of a suitable invariant in the domain is undecidable, or at least has high complexity. Unfortunately, such undecidability proofs are hard.

In contrast, complete methods have many desirable properties. We have illustrated this with the example of LRU cache analysis (Section 2). Since the end result has a unique definition, there is no brittleness, and several algorithms or implementations can be used to compute the same result, which allows performance comparisons (the meaning of performance differences between methods producing non comparable results is unclear) and validation by testing that results are identical.

The exception is the "pseudo round robin" cache replacement policy, which involves a counter global to all cache sets. Also, dependencies between cache sets arise if one considers an integrated model of the processor pipeline, pre-fetching units, and caches, because whether or not some data is in some cache set influences indirectly whether some other data will be fetched.

Some processors implement this LRU policy. It is however more common to implement so-called pseudo-LRU policies, meant to provide the same practical performance at a fraction of the hardware cost[START_REF] Heckmann | The influence of processor architecture on the design and the results of WCET tools[END_REF][START_REF] Gille | Study of Different Cache Line Replacement Algorithms in Embedded Systems[END_REF]. Unfortunately, these pseudo-LRU policies are very different from the point of view of cache analysis[START_REF] Heckmann | The influence of processor architecture on the design and the results of WCET tools[END_REF][START_REF] Monniaux | The complexity gap in the static analysis of cache accesses grows if procedure calls are added[END_REF][START_REF] Monniaux | On the Complexity of Cache Analysis for Different Replacement Policies[END_REF].

This non-monotonicity resulted in some "false alarms", that is, warnings about nonexistent problems, one some industrial programs.

[START_REF] Cousot | Abstract Interpretation Frameworks[END_REF] More generally, advanced tools such as Astrée implement complex combinations of domains and iteration strategies, which make many examples of non-monotone behavior and false alarms on simple analyses actually

succeed.[START_REF] Cousot | Automatic Discovery of Linear Restraints among Variables of a Program[END_REF] Consider the representation of the polyhedron as a system of inequalities l i (x 1 , . . . , xn) ≤ B i , create a new variable y i for any linear combination l i for which there does not already exist a variable, keep this inequality as y i ≤ B i and y i = l i (x 1 , . . . , xn), then the set defined by the original system is the projection of the set defined by the transformed system on the x variables

https://github.com/Z3Prover/z3