Francis Adzima

Py Y Manach

Tudor Balan
email: tudor.balan@ensam.eu

Laurent Tabourot

Sébastien Toutain

Jean-Luc L Diot

Simulation of ultra-thin sheet metal forming using phenomenological and crystal plasticity models

come   L'archive ouverte pluridisciplinaire

Introduction

The ongoing trend on device miniaturization has increased the demand for microparts. Micromanufacturing technologies are being developed and optimized continuously to meet the needs of industry. Particularly, micro-forming processes can achieve high productivity and low costs while delivering good parts. However, the knowledge and the simulation tools that are successfully operated in macro-scaled forming cannot be applied to accelerate the design phase of micro-forming processes. Indeed, the employed ultra-thin sheet metals exhibit peculiarities in their deformation behavior. In the present study, crystal plasticity and phenomenological material models are implemented for the numerical simulation of industrial micro-forming processes. Their predictive capabilities are compared and the influencing parameters for accurate predictions are investigated.

Materials and demonstrator parts

The current study focuses on two demonstrator parts manufactured from metal sheets of CuBe2 and CuFe2P copper alloys. The gull-wing microelectronic package leads (see Fig. 1 (a)) are made of 127 µm thick CuFe2P sheets while 100 µm thick CuBe2 sheets are processed to produce the connector displayed on Fig. 1 (b). Thus, due to the part geometries, small (chip lead) and large (connector) R/T (radius/thickness) ratios situations are under investigation.

Material models and parameters

Phenomenological modeling as well as crystal plasticity based modeling are considered since there is an average of 20 grains through sheets thickness. The crystal plasticity model is based on the works of [START_REF] Peirce | An analysis of nonuniform and localized deformation in ductile single crystals[END_REF] and [START_REF] Raphanel | Three-dimensional rate-dependent crystal plasticity based on Runge-Kutta algorithms for update and consistent linearization[END_REF]. The slip rate

. () s
 on a slip system (s) is given by a visco-plastic law as a function of the resolved shear stress on the system () s  and its resistance to motion ()

s c  , reading (1)
In the previous equation, . () 0 s  is the reference strain rate and n is the strain rate sensitivity coefficient. The hardening at slip system level is based on the so-called PAN (Peirce-Asaro-Needleman) model [START_REF] Peirce | An analysis of nonuniform and localized deformation in ductile single crystals[END_REF].

The phenomenological plasticity model is elastic-visco-plastic and accounts for material anisotropy through the yield function of Bron and Besson [START_REF] Bron | A yield function for anisotropic materials Application to aluminum alloys[END_REF]. Mixed hardening is assumed with an isotropic hardening of Hockett-Sherby and a combined Armstrong-Frederick & Prager kinematic hardening. The phenomenological and crystal plasticity approaches were implemented in ABAQUS/Standard.

The sheets behavior was characterized through mechanical tests and a large experimental database was established. The latter consisted of tensile tests in 0, 45 and 90° to the rolling direction (R.D.), equi-biaxial loading via hydraulic bulge test, monotonic shear tests and reversed shear tests after different pre-strain levels. The crystal plasticity model was calibrated on the tensile test in the R.D. on a Representative Volume Element which grains orientations were sampled from EBSD measurements. As illustrated on Fig. 2, the CPFEM (crystal plasticity finite element method) model proved capable of providing good predictions of all experimental tests although it was calibrated on a single tensile test. On the other hand, the phenomenological plasticity model required a calibration procedure, implemented in SiDoLo [START_REF] Cailletaud | Identification and Inverse Problems: A Modular Approach[END_REF], on the full database in order to describe accurately all the experimental strain paths.

Springback predictions and influencing parameters

The two modeling approaches (CPFEM and phenomenological) were applied to the prediction of springback in bending processes. A numerical study was conducted on several parameters of the CPFEM model. Both demonstrator parts were employed in the study.

Crystalline elastic anisotropy was found to significantly impact the CPFEM predictions of springback shapes. Indeed, assuming whether isotropic or anisotropic elasticity yielded very little difference in the forming step as anisotropic discrete slip flow drives the mechanical response in both cases. However, springback is dominated by the elastic behavior and thus affected by the type of elasticity considered. The differences that arise are illustrated on Fig. 3 which represents the profiles of a portion of the connector part. Furthermore, the ratio T/D (sheet thickness over grain diameter) which can be approximated to the number of grains in the thickness was linked to the amount of springback predicted. When this ratio and consequently the number of grains of grains in the part were increased, the predicted mechanical response converged towards the prediction of the phenomenological model. The CPFEM springback angle also grew closer to the phenomenological value and the experimental measurements (see Fig. 4). Fig. 4: Predictions and measurement of the springback angle of the chip leads A similar trend on the evolution of springback with the T/D ratio was reported by [START_REF] Jiang | Grain Size Effect on the Springback Behavior of the Microtube in the Press Bending Process[END_REF] in their experimental study.

Summary and conclusions

The deformation behavior of copper alloys ultra-thin sheets was investigated with phenomenological plasticity and CPFEM models in the framework of micro-forming. The CPFEM model was more robust with respect to parameter identification. The springback predictions of the CPFEM model were highly sensitive to the assumption made on elasticity (isotropic or anisotropic). Moreover, increasing the number of grains in the sheet thickness led to more springback amplitude and to converging predictions with the phenomenological models in terms of stress fields and springback angle.

Fig. 1 :

 1 Fig. 1: Demonstrator parts: (a) CuFe2P made LQFP package (small R/T ratio) and (b) CuBe2 made connector (large R/T ratio)

Fig. 2 :

 2 Fig. 2: Predictions by the CPFEM model of (a) tensile tests and (b) shear and bulge tests

Fig. 3 :

 3 Fig. 3: Influence of elastic anisotropy on the profile (b) of a portion of the formed connector (a)


Romain Canu

Stefano Puggelli

Mohamed Essadki

Benjamin Duret
email: duret@coria.fr

Thibaut Menard

Marc Massot

Julien Reveillon

F X Demoulin

Where does the droplet size distribution come from?

Keywords: interface, curvature, DNS, two-phase flows

This study employs DNS of two-phase flows to enhance primary atomization understanding and modeling to be used in numerical simulation in RANS or LES framework. In particular, the work has been aimed at improving the information on the liquid-gas interface evolution for modeling approaches, such as the Eulerian-Lagrangian Spray Atomization (ELSA) framework. Even though this approach has been already successfully employed to describe the complete liquid atomization process from the primary region to the dilute spray, improvements are still expected on the derivation of the drop size distribution (DSD).

The main aim of the present work is the introduction of a new framework to achieve a continuous description of the DSD formation during the atomization process. The attention is here focused on the extraction from DNS data of the behavior of geometrical variable of the liquid-gas interface, such as the mean (H) and Gauss (G) surface curvatures. The use of a Surface Curvature Distribution is also proposed and studied.

A Rayleigh-Plateau instability along a column of liquid and a droplet collision case are first of all considered to analyze and to verify the capabilities of the code to correctly predicting the curvature distributions. A statistical analysis

requires the possibility to count the number of droplets. Thus, it is necessary to separate the liquid phase in a set of discrete elements. However, usual atomization process starts with a continuous liquid flow (for instance a liquid jet) and during the atomization a splitting of the continuous liquid phase occurs. Once it is created, for fixed external conditions and generally considering the whole spray, the NDD may evolve towards an asymptotic state, for which numerous theoretical and experimental works are reported in the literature (Babinsky and Sojka, 2002). To address more complex situations or to determine its function in space and time, a transport equation for the DSD, like the Williams Boltzmann Equation (WBE) (Williams, 1958), is required.

Although the flow inside the injector as well as surface instabilities and the methods to solve the WBE have been the subject of several research works in technical literature, few studies deal with the droplet generation step to estimate the DSD. One reason is the lack of data to justify any proposal as well as the absence of a definition of the NDD during the primary breakup. In-fact, considering a liquid jet before the detachment of a liquid parcel, the notion of a diameter cannot be reasonably introduced. Nonetheless, considering an initial instability over a liquid surface, such unsteadiness should evolve in time and create wrinkles at a certain length scale that is somehow related on a theoretical point of view to the diameter of the droplet that will be then generated. It is clear that a link between these two subsequent steps is missing. Regarding the data in this transition zone, lately there have been great progresses both in experimental techniques [START_REF] Déjean | Experimental study on the influence of liquid and air boundary conditions on a planar air-blasted liquid sheet, Part I: Liquid and air thicknesses[END_REF]Warncke et al., 2017) and in numerical simulations [START_REF] Navarro-Martinez | Large eddy simulation of spray atomization with a probability density function method[END_REF][START_REF] Dumouchel | Multi-scale analysis of simulated capillary instability[END_REF][START_REF] Fox | Large-Eddy-Simulation Tools for Multiphase Flows[END_REF]. However, an extended definition of the DSD should be first introduced to be able to deal with the dense spray region as well as with primary atomization.

It should be pointed out that the same problem appears in defining a mean diameter: it is a moment of the NDD and therefore it can be defined only once the continuous liquid jet has been broken in a set of liquid droplets.

The main goal of the present study is therefore the proposition of an innovative extended definition of the drop size distribution to be used all along the atomization process, leading to the so called Surface Curvature Distribution (SCD). In the following, the definition of the SCD will be firstly introduced and related to some characteristic geometrical properties of the liquid/gas interface, such as the mean (H) and Gauss (G) curvatures. Then, numerical simulations will be used to extract the evolution of such geometrical properties on fundamental test cases. The well-known Archer code [START_REF] Tanguy | Application of a level set method for simulation of droplet collisions[END_REF][START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF], where a combined VOF-Level set approach is used to capture the interface and a ghost-fluid method is applied to represent accurately the jump of variable across the liquid-gas interface, has been used to this end. A database of DNS data focused on interfacial geometrical quantities will be presented to provide valuable information to improve atomization modeling in particular for the Eulerian Lagrangian Spray Atomization (ELSA) model [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF][START_REF] Chesnel | Large eddy simulation of liquid jet atomization[END_REF] in the region where the transition to a Lagrangian frame for the liquid phase is realized.

The paper is structured as follows: the next section is devoted to the description of the SCD and to the introduction of a relation between the SCD and the NDD for a spherical droplet based on surface curvatures. In the second part of the work, a Rayleigh-Plateau and a droplet collision test case will be presented to characterize the two curvatures and to show how the breakup/coalescence process from a liquid jet to the final droplet can be directly described employing such surface topological properties. Finally, an Homogeneous Isotropic Turbulence test case will be considered and the evolution of H and G addressed for several values of liquid volume fraction in order to mimic the atomization process.

An extended definition of the drop size distribution

The purpose of this section is the proposition of an extended definition of the drop size distribution that may applicable all along the atomization process including the primary breakup zone. From the pioneering work of [START_REF] Vallet | Modélisation eulerienne de l'atomisation d'un jet liquide[END_REF]; Vallet et al. (2001), it is known that the definition of the mean diameter can be replace beneficially by the more general mean surface density Σ, that is the surface area per unit of volume. This quantity can be defined everywhere whatever is the liquid phase topology and, combined with the liquid volume fraction α, it gives the mean Sauter diameter once the spray is formed: Essadki et al. (2018a) used high order fractional moments of the DSD for disperse phase, where the size is given by the surface area of droplet, to recover some interface geometrical quantities already used in describing the gas-liquid interface in their previous work [START_REF] Essadki | Adaptive Mesh Refinement and High Order Geometrical Moment Method for the Simulation of Polydisperse Evaporating Sprays[END_REF].

D 32 = 6α Σ . Recently,
These quantities are the volume fraction, the mean surface density and the two averaged Gauss G = κ 1 * κ 2 and mean H = κ1+κ2 2 curvatures, where κ 1 and κ 2 are the two principal curvatures of the surface.

To overcome the limits of the DSD in modeling a gas-liquid interface during the whole atomization process, the first idea here proposed is the exploitation of some surface characteristics to look at the DSD since they carry the information of the size distribution which is hidden all along the wrinkling process. The proposal is to look at the curvature (κ) distribution. In this case the distribution of curvature is F S (κ) and the F S (κ )dκ is the area of surface with a curvature κ ∈ [κ , κ + dκ[. This leads to the definition of a surface curvature distribution (SCD). However, an accurate description of the interface cannot be restricted to one geometrical variable as it is often used in the DSD for disperse phase.

Therefore, the Gauss and mean curvatures are proposed here as phase-space variables, leading to the definition of the novel bi-variate SCD: F S (H, G). H and G are much more meaningful to be used in this context with respect to κ 1 and κ 2 . Indeed, the curvatures κ 1 and κ 2 need to be distinguished by considering an order relation κ 1 ≤ κ 2 , while the curvatures H and G are mathematically well defined and provide some interesting topological properties as proposed in [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF].

The definition of F S (H, G) and its general link with the standard number density function is based on the introduction of the Gauss-Bonnet theorem ( See [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF] for details). Considering M as a 3D object, delimited by a surface Σ(M ), the following relation stands:

x∈Σ(M ) G(x)dS(x) = 2πχ(M ) (1)
where χ(M ) is the Euler characteristic that is equal for any homeomorph object to M. Such quantity is particularly meaningful since for instance it is equal to 2 for objects homeomorph to a sphere, considering a set of droplets homeomorph to a sphere, it can be then used to count the number of particles. In the case of a sphere κ 1 = κ 2 = 2 D and the Gauss curvature can be simply expressed as G = 4π S , where S = πD 2 is the droplet surface area. Hence, the two functions GF S (H= 2 √ G,G) 4π and F D (D) are equivalent, since they count the number of droplets per volume and per radius or Gauss curvature. The demonstration of this equivalence between the two distributions in a more general case can be found in [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF]. In Figure 1, the link between the DSD and for instance the probability distribution of H is shown for a cloud of droplets with a known normal distribution (i.e. D mean =5e-05 m and σ=1e-05). A continuous link between the two distributions can be introduced and the DSD can be easily recast starting from the corresponding curvature values. Nonetheless, the SCD is defined all along the atomization process and therefore, once the spray is formed, the DSD is a particular case of F S (H, G). We should mention that the SCD is a distribution that measures the probable surface area having some given interfacial properties. As a matter of fact, it can not cope with the volume occupied by the liquid or the gas phase. Another approach to characterized the spray consists in using a volumetric distribution.

For a spray, the DSD becomes F V (D) in such way that F V (D )dD is the liquid volume fraction of droplets with a diameter ∈ [D , D +dD[. This volumetric distribution can be also extended to any liquid topology, leading to the volumetric distance function distribution (V ϕD), that can be a second proposal to describe the transport of curvatures. It is defined by F V (ϕ) in such way that F V (ϕ )dϕ is the volume fraction of flow that is at a distance ϕ from the interface, with ϕ ∈ [ϕ , ϕ + dϕ[, ϕ is positive in the liquid phase and negative in the gas phase.

With this definition the following relation stands: The gas volume fraction is The V ϕD gives the distribution of the flow with respect to the distance to the interface. The link with the DSD is not straightforward. However, this distribution can be clearly identified on a given geometric object and in particular for a spherical droplet.

α g = 0 -∞ F V (ϕ)dϕ, the liquid volume fraction is α l = +∞ 0 F V (ϕ)dϕ,
A similar idea has been also explored by [START_REF] Dumouchel | Towards an interpretation of the scale diffusivity in liquid atomization process: An experimental approach[END_REF][START_REF] Dumouchel | Multi-scale analysis of simulated capillary instability[END_REF] with the concept of scale distribution, E 3 (d). The definition of this function for a given object (for instance a droplet) is based on the total object volume V 0 and the volume V (d) defined by all points at a distance d or greater to the object surface. The scale distribution is given by E

3 (d) = V0-V (d)

V0

. The link between the scale distribution and the V ϕD can be express readily by

E 3 (d) = F V (ϕ=d) α l
. [START_REF] Dumouchel | Towards an interpretation of the scale diffusivity in liquid atomization process: An experimental approach[END_REF] have worked on experimental data in 2D and more recently they have used numerical simulations to extract 3D results [START_REF] Dumouchel | Multi-scale analysis of simulated capillary instability[END_REF]. A connection with curvatures is surely present but further investigations are required on this point.

Following the idea of [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF], the emphasis of this work is to analyze the surface curvature distribution and the curvature behavior (Gauss, mean and principal curvatures) in multiple two-phase flows configurations. These quantities give an accurate geometric description of the interface behavior. This article aims to propose a new framework to study the atomization process, based on the curvature and SCD time evolution.

Governing equations and numerical methods

Navier-Stokes equations and interface capturing method

Here, a joint level set/VOF method are coupled with a projection method to carry out the direct numerical simulation of incompressible Navier-Stokes equations:

∂V ∂t + (V.∇)V = - ∇p ρ + 1 ρ ∇.(2µD) + f + 1 ρ σHδ(ϕ)n (2)
where p is the fluid pressure, V is the velocity vector, and D is the viscous deformation tensor. At the interface, the surface tension force can be considered based on the Dirac function δ(ϕ). Concerning the interface capturing method, a CLSVOF interface capturing method is used to ensure mass conservation and a accurate description of the interface topology (See [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] for details).

This work benefits for the latest improvements of the ARCHER code : the convective term is solved using the improved Rudman technique presented in Vaudor et al. ( 2017), allowing a better accuracy and description of high density ratio flows. This method is based on the computation of mass fluxes from the VOF method, which can then be used in a conservative formalism of the convective term. The diffusive term is now computed with the approach of Sussman et al. ( 2007), the main advantage of this method is the implicit jump condition of the viscous tensor. For further details, a review of the numerical methods dedicated to the computation of the viscous term in two-phase incompressible flows is available in [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid Method[END_REF].

Fluid dynamics equations are solved in the context of a low Mach number approach, based on a projection method for the direct numerical simulation of incompressible Navier-Stokes equations (detailed in Vaudor et al. ( 2017)).

The viscosity depend on the sign of the level set function according to each phase (liquid and gas). To finalize the description of the two-phase flow, jump conditions across the interface are taken into account with the ghost fluid (GF) method. In the GF method, ghost cells are defined on each side of the interface [START_REF] Kang | Boundary conditions capturing method for multiphase incompressible flow[END_REF]; [START_REF] Liu | A Boundary Condition Capturing Method for Poisson's Equation on Irregular Domains[END_REF]). This prolongs each phase to allow smooth derivatives in the vicinity of the interface. As defined previously, the interface is characterized through the distance function, and jump conditions are extrapolated on a few nodes on each side of the interface.

Curvature computation

Two methods have been used to extract the curvature distribution. The first one is based on the distance function (ϕ) that is part of the numerical procedure in the resolution of the two-phase flow Navier-Stokes equations. The normal to the surface n = -∇ϕ | ∇ϕ| is first calculated. Then, the matrix P = I -n n T with I the identity matrix, and the Hessian matrix H m = ∂ 2 ϕ ∂xi∂xj are calculated in order to have the matrix G m = -P HmP | ∇ϕ| . Finally, the two principal curvatures κ 1 and κ 2 are obtained as

κ 1 = T + 2 √ 2F 2 -T 2 2 and κ 2 = T -2 √ 2F 2 -T 2 2
with T , the trace of the matrix G m and F , the Frobenius norm of the matrix G m . The reader interested in a more detailed explanation of such method is addressed to [START_REF] Kindlmann | Curvature-based transfer functions for direct volume rendering: Methods and applications[END_REF]. The second one has been instead developed in [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF], where the gas-liquid interface is discretized with a 2D triangulated mesh using the Marching Cube algorithm. The curvature is computed at each vertex of the generated mesh through a spatial-averaging process that preserves some topological feature of the gas-liquid surface such as the Gauss-Bonnet formula [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF]. The topological features of the interface in terms of computing the droplet number from the bi-variate SCD, while filtering the noise, are in this way guaranteed. On the other hand, results for the V ϕD have been obtained directly through the distance function.

Results and discussion

This section is dedicated to the analysis of the curvature behavior and statistics in two different configurations. First, a Rayleigh-Plateau instability and a droplet collision configurations are investigated to illustrate our analysis on simple validation cases. Then, a more complex configuration closer to the atomization regime (high Weber and Reynolds number) presenting liquid structures of various shape and size is studied. Namely, the two-phase flows Homogeneous Isotropic Turbulence configuration.

Rayleigh-Plateau instability

Numerical configuration

As a first configuration, the Rayleigh-Plateau instability, where a set of droplets is created from an initial column of liquid, is reported. A quarter of a cylinder has been here computed and symmetric boundary conditions have been employed in a computational box of 1, 5.10 -4 × 1, 0.10 -4 × 1, 0.10 -4 m length. A mesh sizing of 96 × 64 × 64 elements as in [START_REF] Ménard | Développement d'une méthode Level Set pour le suivi d'interface, application de la rupture de jet liquide[END_REF] has been chosen. This is a standard test case for the ARCHER code that is able to correctly reproduce the dispersion analysis diagram on this configuration [START_REF] Ménard | Développement d'une méthode Level Set pour le suivi d'interface, application de la rupture de jet liquide[END_REF]. The liquid and gas properties are reported in Table 2. The cylinder has a radius R = 3.34.10 -5 m and the initial perturbation, essential to observe the instability, has an amplitude of 10% of the radius of the cylinder and a wavelength of R = 3.0.10 -4 m (twice the length of the box). The initial conditions have been chosen in order to have a wave number satisfying k w R = 0.7 which correspond to the fastest growth rate.

ρ l µ l ρ g µ g σ kg.m -3 kg.m -1 .s -1 kg.m -3 kg.m -1 .s -1 kg.s -2 1000 1, 0.10 -3 1.0 1, 879.10 -3 0, 072

Comparison between both curvature evaluation method

Both method aforementioned to extract the curvature distribution are tested on one time frame of the Rayleigh-Plateau test case, but no significant differences are observed as shown in Figure 2. The mean and Gauss curvatures are provided as dimensionless using the radius of the initial cylinder as reference length, whereas for the Gauss curvature the radius squared has been used to this end. Considering that the two approaches behave in a completely consistent manner, in the following all results have been obtained using Method 1. The parabola represents the spherical shape (i.e. k 1 =k 2 ) and the dotted black line is the cylinder (G = 0).

Curvature evolution and statistics

To have a deep insight in the evolution of the breakup event, four time

represents R 1 R 2 with R 1 = 1
κ1 and R 2 = 1 κ2 , the two curvature radii. The final graph represents the (κ 1 -κ 2 , κ 1 + κ 2 ) space which is the (κ 1 , κ 2 ) graph with a 45˚rotation. The dashed lines correspond to a cylindrical topology ; it means that one of the two main curvatures is equal to zero and therefore the Gauss curvature G. The solid lines correspond to a spherical topology so κ 1 = κ 2 . On the (H, G) and (S, P ) graphs, it forms a parabola and no points can be inside this parabola. By the way it is calculated, κ 1 is always greater than κ 2 so κ 1 -κ 2 is always positive. On the (κ 1 -κ 2 , κ 1 + κ 2 ) graph, the dash-dotted line separates "bubbles" on the left from "droplets" on the right.

Besides, on the (H, G) and (S, P ) graphes, the limits of the mesh resolution are also represented by dotted lines. In the following, only the (H, G) graph and the (κ 1 -κ 2 , κ 1 + κ 2 ) graph are conserved, mainly due to their mathematical properties (see Section 2). On Figure 9 and 10, the distance function ϕ is also reported to better understand the obtained results and it is dimensionless using dx as reference and the curves are shown between -l x and l x .

As long as the cylinder is weakly deformed (Figure 9), in terms of PDF profile for ϕ, points are very close to the theoretical shape of a cylinder with a linear curve inside the liquid part. In the figure, line without markers represents the initial cumulative profile for a cylinder, whereas the dashed line separates liquid (positive values) and gas (negative values). For the cumulative of total curvature (i.e. 2H) in Figure 8-left and Gauss curvature in Besides, on this graph, the points, which represent the big part of the cylinder, get closer to the parabola representing the sphere.

When the break-up happens (t c ), such arc breaks up and two clouds of points appear (see Figure 5-bottom-left) with one set located on the black curve.

These points correspond to the big droplet. Indeed, the mean curvature is approximately 0.5 and the Gauss curvature about 0.25 which correspond to a radius twice as big as the initial jet radius. The other points, located below the parabola, correspond to the small ellipsoid. Here, the mean curvature is quite high because of the small thickness and the Gauss curvature is lower because the curvature is close to zero along the length. The big droplet appears also as the step on the total curvature cumulative shown in Figure 8-right and the small ellipsoid by all the high positive values. Such break-up process is mainly Finally, when the two droplets are generated (t d ), for the mean and Gauss curvatures graph (see Figure 5-bottom-right) most of the points are on the spherical curve. The points corresponding to the big droplet are located approximately at the same place than the previous time. Conversely, points corresponding to the small droplet are those having an higher curvature. On Figure 10, PDF profile becomes now parabolic in the liquid part, as expected for spheres, with higher values of the distance function because the generated droplet is bigger than the initial cylinder. On the (κ 1 -κ 2 , κ 1 + κ 2 ) graph (Figure 6), same comments than on the (H, G) graph can be done ; the arc when the cylinder is deformed, the two clouds of points when the break-up happens and the location of the points with respect to the spherical and cylindrical shapes.

Hence, the analysis of the evolution of Gauss and mean curvatures can be reliably employed to describe a breakup process from a liquid jet to the final droplet. Distributions of H and G follow the physical process here investigated suggesting that their evolution all along the atomization process can be characterized.

Curvature evolution and statistics

In Figure 11 the resulting evolution of the liquid surface on four time frames is shown, where it should be pointed out the initial configuration (t a ), the moment of coalescence (t b ), the generation of one deformed droplet (t c ) and 325 its subsequent breakup with the generation of a torus structure (t d ). In Figure 12 the corresponding evolution on the (H, G) plane is reported, where it can be pointed out that, similarly to the Rayleigh-Plateau test case, the curvatures evolve consistently with the investigated physical process.

Considering this preliminary assessment on fundamental test cases, in the 330 next section, the same surface geometrical properties will be used to analyze a Homogeneous Isotropic Turbulence (HIT) box where breakup and coalescence events appear together due to the turbulence forcing.

Homogeneous Isotropic Turbulence configuration

Numerical configuration

335

A further assessment and exploitation of the procedure regarding curvature extraction has been performed using a three dimensional cubical domain with periodic boundaries similar to the investigation of [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF].

This configuration represents a extension of single phase mixing studies to two phase flow and it has been already analyzed to improve the understanding of 340 primary atomization (see [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]). The mesh resolution used in the present investigation is based on previous studies on the same configuration and on single-phase flow simulations to capture the Kolmogorov length scale [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]. In [START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF], two mesh sizing are compared on the same configuration and in the present work the fine computational domain counting 256 3 mesh elements has been employed.

The validation of this set-up for single-phase simulation is not here shown for the sake of clarity. The interested reader is addressed to [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF], where this topic is discussed in detail.

In order to define a configuration with realistic interactions between the liquid surface and the turbulent gas flow-field, the following dimensionless parameters σ the surface tension, ν the kinematic viscosity in liquid or gas phase, ρ the density and L is the length of the box. In Table 3 the corresponding set of physical properties is reported.

Several values of the liquid volume fraction α l inside the box have been studied

in the present work for curvature extraction (i.e. ranging from 1% to 90%) in order to mimic the whole atomization process as shown in Figure 13 using the data obtained from DNS calculation on a diesel jet. The main aim is the extraction of curvatures evolution all along the jet evolution in order to point out firstly the link between the DSD and the F S (H, G) and also how curvatures evolve from the dense region up to the dispersed spray (i.e. reducing α l ).

Curvature evolution and statistics

Results obtained for the HIT varying α l are here presented starting from the 1% test case, which can be reliably considered as representative of the dilute spray region. A snapshot of the liquid-gas interface distribution together with the evolution in time of Σ and H, which represent the volume average for the whole domain of liquid-gas interface and mean curvature, as well as their time averaged values is reported in Figure 14. In the following, t * is the time culation the liquid-gas interface associated for each class of H. Such procedure has been gathered for nearly 15 τ * in order to reduce the statistical noise and to obtain a time-averaged distribution. Black dotted lines represent the mesh sizing, which can be considered as a numerical limit in curvature calculation. A similar evolution has been obtained also for the PDF(G), which is not shown here for the sake of clarity. The peak of PDF(H) is located at positive values since a dilute spray is here considered, whereas the probability to have gas bubbles (i.e. H < 0) is negligible. Starting from such distribution, much more information about the evolution of the liquid-gas interface can be now extracted.

For instance, the spray NDD can be now derived, associating to each curvature class (H i ) the corresponding number of spherical droplets (N d.i ) characterized by the class amount of surface (A i ) collected during the calculation:

N d,i = A i 4π(H i ) -2 (3)
This represents clearly an abrupt approximation since each curvature class is associated to the corresponding spherical shape, but it highlights a direct link

P DF tot = V 1 P DF 1 + V 2 P DF 2 where V 1 and
V 2 represent the liquid volume associated to each class, has been obtained (See Figure 16). It is worth noting that V 1 represents nearly the 5% of the total liquid volume, showing that, even for such low value of liquid volume fraction, the contribution of collisions is not at all negligible. The resulting log-normal P DF tot , both in terms of NDD or VDD, reliably agrees with DNS data.

Hence, starting from the H-PDF it has been possible to derive a NDD that is consistent with a standard drop-size distribution in the dilute regime like the log-normal one. Furthermore, it is worth pointing out that the mean value calculated for P DF tot is coherent with the characteristic Sauter Mean Diameter (SMD) that can be calculated as SM D = 6α l /Σ. Such evidence enforces the idea that curvatures PDF can be reliably used to characterized the spray distribution both in terms of mean value and dispersion.

This comparison clearly shows that, in a test article nearly representative for a dilute spray (i.e. α l = 1%), a direct link between the spray NDD and the curvature evolution can be determined. However, beyond the characterization of the spray distribution, such DNS analysis can bring also further information about the curvature evolution that can be useful to better understand the topology of the liquid-gas interface under investigation and to give some hints about the curvature modeling. In particular, in Figure 17 the time averaged joint PDF(H, G) is shown together with the time averaged joint PDF of the (κ 1 -κ 2 , κ 1 + κ 2 ) space. On the joint PDF(H, G), the parabola represents the spherical shape (i.e. H 2 = G) whereas the dashed line specifies a cylindrical topology (i.e. G = 0). On the joint PDF(κ 1 -κ 2 , κ 1 + κ 2 ), the dashed lines represent a cylindrical shape and the spherical shape is shown by the horizontal solid line. The vertical dash-dotted line separates bubbles topology from droplets topology. It is possible to extract several interesting information about the evolution of the liquid-gas interface looking at such joint PDF distributions of surface curvature:

4 3 2 1 0 1 2 3 4 G[ -] 3 2 1 0 1 2 3 H[ -] Cylindrical Shape Spherical Shape 4 2 0 2 4 6 8 κ1 + κ2[ -] 1 0 1 2 3 4 κ1 -κ2[ -] C y li n d r ic a l S h a p e Spherical Shape
• the PDF(H, G) shows that, in terms of time-averaged distribution, the topology of the liquid surface is mainly stretched along two directions, i.e. several points are collapsed on the parabola representative of the spherical shape and in a similar manner also the cylindrical asymptotic limit appears. Such bi-modal topological distribution derives from the presence of big cylindrical ligaments that, subjected on the action of the turbulence field, collide and generate a cloud of spherical droplets.

• the PDF(κ 1 -κ 2 , κ 1 + κ 2 ) gives similar information in terms of timeaveraged properties with respect to the PDF(H, G), but, looking now at the contour iso-lines, it should be pointed out that the most probable value of the two principal curvatures is not associated to a spherical shape (i.e. it not collapsed on the κ 1 -κ 2 = 0 line).

Hence, the joint PDFs, which can be extracted from surface curvatures, may

give an idea about the topological evolution of the two phase flow together with some hints about the generated spray population. Considering these results, together with data already presented on the Rayleigh-Plateau test case, the capabilities of surface curvatures in describing the liquid breakup as well as its direct link with the NDD have been clarified.

Therefore, higher values of liquid volume fraction have been then investigated in order to understand how the curvatures evolve and if it is possible to determine a general behavior of the curvatures PDFs all along the atomization process.

In Figure 18 the liquid-gas interface distribution together with the time averaged PDF(H) is shown for α l = 5%. Clearly, much more ligaments are generated now due to the augmentation of the liquid volume inside the box with respect to the previous case: elongated and stretched cylindrical elements appear and from their turbulent breakup some spherical droplets are produced. The PDF(H)

shows here an evolution in line with Figure 16 with a mode value shifted towards negative values of H. Now, the probability of finding H < 0 is no more zero and this is linked to the presence of regions where gas bubbles are surrounded by the liquid phase leading to an opposite curvature sign. In Figure 19 the joint PDF(H, G) and PDF(κ 1 -κ 2 , κ 1 +κ 2 ) are reported for the 5% case. Considering the PDF(H, G), with respect to Figure 17 less points tends now to collapse along the parabola and the probability peak is strongly stretched along the line at G = 0. This is again representative of the topological variation in the liquid shape from a sphere to a deformed and elongated cylinder while increasing the liquid volume fraction. In a similar fashion, looking at the PDF(κ 1 -κ 2 , κ 1 + κ 2 ) the distribution is no more just collapsed on the spherical

4 3 2 1 0 1 2 3 4 G[ -] 3 2 1 0 1 2 3 H[ -] Cylindrical Shape Spherical Shape 4 2 0 2 4 6 8 κ1 + κ2[ -] 1 0 1 2 3 4 κ1 -κ2[ -]
C y li n d r ic a l S h a p e

Spherical Shape

Figure 19: Joint PDF of the surface curvatures for α l = 5%. Left : joint PDF(H, G), right :

joint PDF(κ 1 -κ 2 , κ 1 + κ 2 ).
topology, but tends to be straightened also along the line corresponding to a cylindrical shape (κ 2 = 0), with the PDF peak moved away from the spherical asymptote with respect to the previous test point. This clearly shows that for α l = 5%, the interface has a wider range of possible curvatures than the previous α l = 1% configuration, and is subject to more deformation and wrinkles.

Increasing further the liquid volume fraction, as shown in Figure 20 from the instantaneous snapshots from α l = 10% up to 90%, an important number of continuous liquid structures appear and a set of bubbles is at the end generated for the very dense case. The time averaged PDF(H) as well as the joint and 22. Several considerations can be derived from these plots. For α l = 10%, some pockets of gas are closed inside the continuous liquid structures and indeed, with respect to previous cases, now the probability of having H<0 is no more negligible. , bottom : α l = 90%.

The probability distribution of the mean curvature has a similar shape with respect to Figure 16 with the mode value that is shifted towards lower H-values.

The topology of the the interface is now strongly modified and the PDF(H, G)

points out a leading cylindrical shape together with two clouds of points collapsed on the parabola. They represents the generation of a set of bubbles and droplets nearly spherical. In terms of principal curvatures, the peak is moved from the κ 1 = κ 2 limit (κ 1 -κ 2 = 0), even if a significant stretch along the spherical line can be still pointed out.

As soon as the liquid volume fraction is further augmented, points characterized by a positive mean curvature disappear since the probability of having some iso-lated droplets is now negligible. The imposed level of turbulent kinetic energy is not high enough to disrupt the ligaments and the liquid-gas interface preserves its continuity. The gas phase is trapped inside these structures and on the (H, G) plot a huge amount of bubbles collapses on the spherical topology.

Coherently, the principal curvatures change their sign. Finally, at α l = 90%

the characteristic values of H and G are higher (i.e. bubbles with a smaller radius) and the PDF(κ 1 -κ 2 , κ 1 + κ 2 ) changes completely its orientation and the system is stretched on the negative principal curvatures space. The HIT now behaves in terms of topology in a complementary manner with respect to the 10% case as shown also by the time averaged PDF(H), where the peak is moved towards negative values because of the generated smaller bubbles.

By comparing all the results shown so far, it is interesting to point out that the curvature evolution is completely consistent with the observed physical phenomena. Both the PDF(H) and the joint curvature probability distributions evolve all along the atomization process maintaining the same shape and moving the peak of the distribution from positive to negative values. Moreover, as far as the spray becomes dilute a link between the curvature evolution and a standard DSD is directly available. This surely represents a novel and general framework to study the atomization process and the obtained mathematical shape of the PDF(H) is particularly promising having in mind the introduction of such characteristics of the liquid-gas interface in the modeling context (for instance in the ELSA formalism).

Conclusions

This study shows how DNS of two-phase flows can be employed to improve understanding and modeling of primary atomization to be then used in RANS or LES framework. In particular, the study is aimed at enhancing the information available in ELSA framework through the introduction of liquid/gas interface curvatures. Firstly, on a mathematical point of view, an extended definition of the drop size distribution to be used all along the atomization pro-cess is proposed, leading to the so called surface curvature distribution. Hence, the extraction of surface curvatures is shown for simple configurations such a Rayleigh-Plateau instability and a droplet collision. Two different methods to extract the curvatures have been tested and validated. The evolution of Gauss and mean curvatures is analyzed in detail in order to show how the breakup process from a liquid jet instability to the final droplet can be easily described in terms of such surface geometrical properties. Then, the same procedure has been applied on the HIT, where several cases, varying the liquid volume fraction, have been studied. The analysis has been firstly focused on a dilute case (i.e.

α 1 = 1%) where the link between the DSD and SCD has been clarified. Indeed, starting from the PDF(H) it has been possible to extract a NDD and VDD that follow a standard log-normal distribution. By investigating higher values of the liquid volume fraction (up to 90%), it has been shown that the SCD is defined even if droplets/bubbles or large liquid/gas structures are present. Joint PDFs of H and G, as well as joint PDF of κ 1 -κ 2 and κ 1 +κ 2 , have been investigated in this configuration and give a good insight on the behavior of the two-phase flows and the evolution of liquid (or gas) structures shapes. It has been found that curvature distributions are able to describe the transition from a dilute spray to a bubbly flow and in particular the PDF(H) shows a very similar shape in all the studied cases.

Acknowledgements

This work was granted access to the HPC resources of IDRIS, TGCC and CINES under the allocation A0032B06153 made by GENCI (Grand Equipement National de Calcul Intensif).

Babinsky, E., Sojka, P. E., 2002. Modeling drop size distributions. Progress in energy and combustion science 28 (4), 303-329.

  and its subsequent combustion nowadays still represents one of key source for energy production both in aero-engine and automotive context. Considering the difficulties and the high costs of experimental measurements for two phase flow mainly in the dense spray region, in last years numerical simulations have been reliably used to gain a detail insight in liquid evolution. Several numerical methods, associated to different computational costs, have been therefore developed to this end (see for instance Gorokhovski and Herrmann (2008); Demoulin et al. (2013) and references therein) and they are normally aimed at providing predictive tools to be used to understand the various physical mechanisms involved in the atomization context. A standard output normally expected from any atomization model or theory, going from DNS studies till experimental correlations is the drop size distribution (DSD). From a theoretical point of view, several definitions of this function can be introduced. Generally, it is F N (D) in such way that F N (D )dD is the number of droplet per unit of volume with a diameter ∈ [D , D + dD[. In this case, it is the number diameter distribution (NDD) and it may be also normalized to define a probability density function P N (D). The function F N (D)

Figure 1 :

 1 Figure 1: Link between the DSD and the SCD for a known cloud of spherical droplets.

  and the mean surface density Σ = F V (0).

Figure 2 :

 2 Figure 2: Rayleigh Plateau instability at t = 1, 88.10 -4 (Top) and mean curvature in function of Gauss curvature with the first method (left) and with the second method (right).

  Figure 7-left, a step profile at values corresponding to this cylinder is derived. On the graph showing mean curvature in function of Gauss curvature (Figure 5-top-left), all the points are at this time step all gathered around the vertical line at G=0. Then, when such cylinder is deformed (t b ), negative values for Gauss curvature appear (Figure 7-right) related to all the saddle points in the middle of the cylinder. This portion of the cylinder is visible on the mean and Gauss curvatures graph (Figure 5-top-right) with corresponding points that create an arc.

Figure 3 :

 3 Figure 3: Top left : Rayleigh Plateau instability at ta = 3, 92.10 -7 s ; Top right : Rayleigh Plateau instability at t b = 1, 49.10 -4 s ; Bottom left : at tc = 1, 86.10 -4 s ; Bottom right : at t d = 2, 03.10 -4 s

  n d r i c a l S h a p e C y li n d r i c a l S h a p e Spherical Shape : Droplet Spherical Shape : Bubble

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Different spaces for curvature representation.

Figure 7 :Figure 8 :

 78 Figure 7: Cumulative of the Gauss curvature (Left : at ta ; Right : at t b )

Figure 9 :

 9 Figure 9: Left : PDF of the distance function at ta; Right : cumulative of the distance function at ta. The lines with • represent the numerical results, the solid lines represent the theoretical shape for a cylinder and the dashed lines separate gas (on the left) from liquid (on the right)

Figure 10 :

 10 Figure 10: Left : PDF of the distance function at t d ; Right : cumulative of the distance function at t d . The lines with • represent the numerical results, the solid lines represent the theoretical shape for a cylinder and the dashed lines separate gas (on the left) from liquid (on the right)

Figure 11 :Figure 12 :

 1112 Figure 11: Top left : Surface evolution for droplet collision test case ; Top left : at ta = 2.97.10 -5 s ; Top right: at t b = 7.2.10 -5 s; Bottom left : at tc = 9.9.10 -5 s ; Bottom right : at t d = 4.3562.10 -4 s

Figure 13 :Figure 14 :

 1314 Figure 13: Regions considered for the extraction of curvature distribution all along the atomization process.

375Figure 15 :Figure 16 :

 1516 Figure 15: Time averaged PDF of mean curvature.

Figure 17 :

 17 Figure 17: Joint PDF of the surface curvatures for α 1 = 1%. Left : joint PDF(H, G), right: joint PDF(κ 1 -κ 2 , κ 1 + κ 2 ).

Figure 18 :

 18 Figure 18: Left : Visualization of liquid-gas interface for liquid volume fraction of 5%. Right : time averaged PDF(H).

Figure 20 :

 20 Figure 20: Visualization of liquid-gas interface for liquid volume fraction for 10% (left), 50% (center) and 90% (right).

Figure 21 :

 21 Figure 21: Time averaged PDF of mean curvature. Top left : α l = 10%, Top right : α l = 50%

Table 2 :

 2 Physical properties of the Rayleigh Plateau instability

Table 3 :

 3 Properties for DNS simulations of HIT (S.I. units).

acting on the liquid surface, increases the probability of collisions and determines the final spreading of the size distribution.

PDF(H, G) and PDF(κ 1 -κ 2 , κ 1 + κ 2 ) are shown respectively in Figures 21

Oh = 0.015 has been considered.

Cylindrical Shape Spherical Shape

C y li n d r ic a l S h a p e Spherical Shape

C y li n d r ic a l S h a p e Spherical Shape Vaudor, G., Ménard, T., Aniszewski, W., Doring, M., Berlemont, A., Jul. 2017.

A consistent mass and momentum flux computation method for two phase flows.


Pierre Bernhard

Marc Deschamps

Ieqiwhd Niversit¡

The birth of modern system theory

teaching and research institutions in France or abroad, or from public or private research centers.

1 Introduction udolf imil ulmn pssed wy on the Pnd of tuly PHITF ris nme will remin tied to the ulmn (lterF he purpose of the ulmn (lter is to estimte timeEvrying vriles from noisy nd inomplete mesurementsD using @noisyA model of the underlying dynmisF ulmn ojeted tht his (lter elonged to estimtion theoryD rguing tht this term hd ome to nme ody of nonEsieneF re lter lso denied tht the (lter ws sed P on proilitiesY s muh s he refused to redue it to mere lgorithmD euseD for himD the (lter ws to e knowledged s theory on its ownF feuse the (lter llows one to estimte vrile whih is not diretly mesuredD it hs sometimes een puliized y omputer sientists s softE wre sensorF he (lter nd its extensions re now used in mny dominsD e they industril nd tehnologilD or in the helthD soilD iologilD or erth sienesF st is lso implemented in glol positioning systems @qAF sn the sme yer IWTH when ulmn pulished his rtile out the (lter @ulmnD IWTHAD he lso pulished mjor rtile on liner qudrti optiE ml ontrol @ulmnD IWTHAD whih he hd erlier presented t onferene in IWSWD nd n lso mjor rtile on system theory @ulmnD IWTHAD omE pleted y ulmn @IWTPD IWTQAF hese ontriutions were groundEreking for ontrol theoretiinsF hey trnsformed ontrol theory nd liner system theoryD nd their in)uene went eyond the mere domin of (ltering nd preditionF ulmnD in this yer IWTHD lso pulished importnt rtiles on the stility of liner dynmil systems @ulmn nd fertrmD IWTHD AF folstered y the dvent of the digitl omputerD the theories presented in @ulmnD IWTHD A were extensively put into use in the eutomti ynthesis rogrm involved in the epollo lunr lnding progrmD whih strted lso Q in IWTHF e new theory hs seldom een dopted so quikly y prtitioners nd for suh signi(nt projetF his hd een the seD thoughD of the iener (lterD whih ws oneived in IWRH for the ntiEirrft rdrD nd hd een emrgoed until IWRW for its militry sensitivityF e shll present here ulmn9s ontriutions in the erly IWTH9sD with fous on ontrol theory nd with toy pplition in popultion eologyF 2 System theory 2.1 A paradigmatic change PFIFI tte of the theory efore IWTH fefore IWTHD system nd ontrol theories were on(ned to liner timeEinvrint @vsA systemsD mostly soElled monovrile onesD whih re hrterized y slr signlsF e @linerA dynmil system is devie whih is exited y timeEvrying signlD lled the inputD to produe timeEvrying signlD lled the ouputF por mthemtiinsD liner system mounts to liner opertorD trnsE forming time funtion into nother oneF feuse this exittion hppens in rel timeD the trnsformtion must stisfy the priniple of uslity @the R urrent output does not depend on future inputsAF woreoverD desription nd nlysis sed on the theory ville t tht time require linerityD timeEinvrineD nd some other properties out ehvior t in(nityF wthemtiins hd developed n elegnt nd powerful wy for hndling suh trnsformtionsD through the use of the @quite esoteriA vple trnsE formD whih leds to the representtion of system through trnsfer funE tionD whih is rtio of polynomils of omplex vrileF en importnt point is tht the trnsfer funtion is esily otined from liner di'erentil eqution governing the systemF he trnsfer funtion hs the property of trnsforming sde of sysE tems @wherey the output of system is the input of the next oneA into simple produts of their trnsfer funtionsF his llows the nlysis of feedE k systemsD where the output is reEintrodued s omponent of the input of the sme systemF his loop is neessry to servomehnismsF he onept of trnsfer funtion lso led to the iener (lterD whihD t tht epohD ws the stndrd tool in signl proessingF he trnsfer funtion nd the iener (lter were intimtely tied to the frequeny response of the systemD tht isD its ehvior if exited y sinusoidl signls of vrious frequeniesF S e lever trik of iener (ltering ws to onsider noisy signl to e (ltered @net its noiseA s the output of liner system exited y noise with dequte sttistil propertiesF ulmn used this move in the ulmn (lterD dvising one of the uthors @ierre fernhrdAX ke the ulmn (lterD whihD s everyody knowsD ws invented y ienerF PFIFP snnovtions in the yer IWTH ulmn represented the trnsformtion of inputs into outputs y the mediE tion of n internl stte of the systemD onsisting of vetor of rel vrilesD whose temporl vritions re governed y (rstEorder di'erentil eqution in ontinuous time or y (rstEorder di'erene eqution in disrete timeF his is why this desription of the system is lled internlD the lssil one then eing lled externlF sn tht representtionD the input ts on the dynmis of the stteD nd the stte instntneously produes n outE putF feuse the stte is vetor nd ll reltions re linerD liner lger is ppropriteF his mthemtil mhinery llows the onsidertion of vetorEvlued inputs nd outputsF woreoverD if the mtries de(ning the sysE tem vry over timeD the system is no longer time invrintF rt of reliztion theoryD whih lies t the hert of liner system theoryD not only onerns linE T er timeEinvrint systemsD ut lso n e extended to systems tht re not invrint over timeD inluding (ltering nd optiml ontrolF he onept of the stte of system ws impreisely known s set of numers from whih the entire future ehvior of the plnt my e deE termined provided tht the future inputs of the plnt re known @ulmn nd uoepkeD IWSVX IVPIAY for exmple positions nd veloities in mehnE il systemD intensities in indutors nd hrges of pitors in n eletril iruitF ulmn @IWTHA mentioned tht it is the smllest suh setF his property ws present in the @erlierA onept of xerode equivlene of forml lnguges nd used in the theory of utomt @ulmn et lFD IWTWAF he intimte link of the onept of stte with (rstEorder di'erentil equtions ws now knowledgedD ut often on(ned to the derivtion of the trnsfer funtion of the liner timeEinvrint systemF he diret use of di'erentil @or di'ereneA equtions in optimiztion hd eome more ommon in the lte (fties under the in)uene of fellmn9s hynmi rogrmming @fellmnD IWSUAF ulmn9s old move ws to implement this use in the very de(nition of liner system nd to investigte its propertiesF he internl desription voids using the vple trnsformD with the onsequene tht systems no longer need to e time invrintF ignls now U mount to time funtionsD whih justi(es the term time dominD in ontrst to the term frequeny dominF he tools developed y ulmn nd fter him re treted with the lulus of vritions for optiml ontrolD nd with wrkov proesses for (lteringF he dvent of the digitl omputer nd of diret digitl ontrol led ulmn nd others to develop disreteEtime theoryD originlly s theE ory of smpled dt systemsD whih onsists of onsidering ontinuousEtime systems only t disrete instntsF he prllel ws nturl nd elegnt in the new theoryF uolmogorov @IWRIAD independently of ienerD hd erlier pioneered disreteEtime version of the iener (lterF he trnsfer funtion of timeEinvrint system in internl form results from simple lgeri formulF gonverselyD (nding the internl representE tion of system given in externl form is deep questionD whih involves throrough nlysis of liner system in internl formF his is the ojet of reliztion theoryF 2.2 Realization theory uhlmn9s ontriution to system theory egn with his rtile yn the generl theory of ontrol systems @ulmnD IWTHAD ut the founding puliE V tionsD whih we shll follow hereD re ulmn @IWTPD IWTQAF he(nition I e reliztion of liner @or 0neA inputEoutput trnsformE tion is representtion in internl form s expressed in iqF @IA nd @PAD or @QA nd @RA elowF vet x ∈ R n e the stte of the system @n is lled the dimension of the reliztionAD u ∈ R m the inputD or ontrolD nd y ∈ R p the outputF he interesting ses re m ≤ n nd p ≤ nF e use xewton9s nottion for time derivtivesX ẋ = dx/dtF e ontinuousEtime system is of the form

ẋ(t) = F x(t) + Gu(t) , @IA y(t) = Hx(t) , @PA nd in disrete time x(t + 1) = F x(t) + Gu(t) , @QA y(t) = Hx(t) . @RA
ulmn rgued ginst the possiility of ugmenting the output eqution with term +Ju(t)F

por the ske of ompletenessD we mention tht if ll three mtries F D GD nd H re onstntD then the system is liner timeEinvrintD nd its trnsfer W funtion is

H(s) = H(sI -F ) -1 G. @SA
grmer9s rule @firkho' nd w vneD IWTUA implies tht H(s) is mtrix of @stritly properA rtionl frtions of sD with the hrteristi polynomil of F s their ommon denomintorD nd thus the poles of H(s) re the eigenvlues of F F iven without referring to the trnsfer funtionD the trnsformtion from input to output indued y these equtions is unhnged when hnging the sis in the stte speD or equivlently if we use the stte ξ = T xD where T is n invertile mtrixF he ontinuousEtime system eomes ξ = T F T -1 ξ + T Gu, @TA y = HT -1 ξ. @UA rene hnging (H, F, G) into (HT -1 , T F T -1 , T G) leves the system unE hngedD with the sme trnsfer funtion if the system is timeEinvrintF his points out tht the sme inputEoutput system my orrespond to new representtionD whih is then not uniqueF he lk of uniqueness lso ours with higher dimensionl vetorD sy zD mde of x nd vetor ξ of ritrry x(t 0 ) = 0D is onto @tht isD surjetiveA for some t 1 > t 0 F e stte x 0 = 0 is unoservle if the output of the free systemD orreE sponding to u(•) = 0D initilized t x 0 D is y([t 0 , t 1 ]) = 0 for ny t 1 ≥ t 0 F e system is ompletely oservle if no stte is unoservleF gonsequentlyD the system is ompletely oservle if the pplition x(t 0 ) → y([t 0 , t 1 ])D whih is liner if u(•) = 0D is one to one @tht isD injetiveA for some t 1 > t 0 F ulmn @IWTHA lso gve e0ient riteri to exmine these propertiesF sn the se of timeEinvrint systemsD the ulmn riteri re expressed in terms of the rnks of omposite mtriesX

IH dimensionX z =     x ξ     , ż =     F 0 A B     z +     G C     u , y = (H 0)z . @VA AD BD
IP heorem I (F G) ompletely rehle ⇔ rank(G F G F 2 G . . . F n-1 G) = n , @WA (H F ) ompletely oservle ⇔ rank                 H HF HF 2 F F F HF n-1                 = n @IHA
ulmn @IWTHA lso gve simple riteri for systems tht re not timeE invrintF hey re less lgeriD more nlytiD ut they shre the property lled dulity @we use prime for trnsposedAD de(ned sX (F G) ompletely rehle ⇔ (G F ) ompletely oservle @IIA (H F ) ompletely oservle ⇔ (F H ) ompletely rehle. @IPA e shll see tht this dulity onerns optiml ontrol nd (lteringF ulmn @IWTHA nlyzed the dulity etween the iener (lter nd the liner qudrti regultorF re pointed out tht it is lso relted to the known dulity etween the di'erentil equtions ẋ = F x nd ṗ = -F p or the di'erene equtions heorem P e reliztion is miniml @hs stte spe of minimum dimenE sionA if nd only if it is nonilF st is unique up to hnge of the sis in the stte speF ulmn @IWTPA further showed tht the stte spe of ny liner system in inE ternl formD not neessrily timeEinvrint @tht isD system suh s iqF @IA nd @PA with mtries HD F D nd G depending ontinuously on timeA n e deomposed nonilly s the diret sum of four suspes "whih n vry over time if the system is not timeEinvrint" s represented in the pigure IX he nonil deomposition of liner system in internl formF oureX ulmn @IWTPAF tht deomposition yields nonil deomposition of the system mtriesD whih tkes the formX

x(t + 1) = F x(t) nd p(t) = F p(t +
F =             F AA F AB F AC F AD 0 F BB 0 F BD 0 0 F CC F CD 0 0 0 F DD             , G =             G A G B 0 0             , @IQA H = ( 0 H B 0 H D ). @IRA
he suspes BD CD nd D re not uniquely de(nedD ontrry to the deomE position of the mtries of the system (H, F, G)D up to hnges of the sis within eh of the four suspesF IS 2.3 Compensator design by pole placement he ulmn (lter ws the (rst oserverD lso lled oserving system in ulmn @IWTHAD who proposed either n optiml oserver or n optiml oserving system in terms of the totl numer of time steps neessry to retrieve the stte in disreteEtime systemF he term oserver ws oined y vuenerger @IWTRAD who extended the onept in less expliit nd more omplited formD ignoring the stility rgument we use elow in prgrph RFIFPF he equtions of the ulmn (lter reD in disrete time s it ppered originlly in ulmn @IWTHAX

x(t+1) = F x(t) + Gu(t) + K(y(t) -H x(t)), @ISA
or in ontinuous time @ulmn nd fuyD IWTIAX

ẋ(t) = F x(t) + Gu(t) + K(y(t) -H x(t)), @ITA
providing n estimte x(t) of the stteD optiml in the sense tht it minimizes the expeted squred L 2 norm of the error signl x(t) = x(t) -x(t)F he nturl ideD proposed y ulmn @IWTHA for monovrile disreteEtime systemsD is to ssoite suh n oserver with ontrol lw he topi overed here is prtilly ddressed in ulmn @IWTHAD ut the uthorittive rtile is ulmn @IWTHAD whose nottion we dopt for the most prtF sn his introdutionD ulmn @IWTHX IHPA sttedX his prolem dtes kD in its modern formD to iener nd rll t out IWRQF re lso ited xewton tr et lF @IWSUAD s representing the stte of the rtF hereforeD in its in(nite horizon @optiml regultorA formD ulmn9s linE er qudrti optiml ontrol theory ws not newF roweverD the solution provided y the theory then ville @xewton tr et lFD IWSUA ws in terms of spetrl ftoriztionD whih is omprle to the iener pilter nd did IV not esily led to e0ient lgorithmsD prtiulrly so for multivrile proE lemsF st ould not e extended to (niteEhorizon prolems eitherF etullyD one (niteEhorizon nonEhomogeneous slrEontrol linerEqudrti optimizE tion prolem is solved in werrim sss @IWSWAD with the orret iti equE tion nd the liner equtions for the nonEhomogeneous termsD lthough the ltter re di0ult to reognizeF QFIFI piniteEhorizon prolem he theory of linerEqudrti optiml ontrol egn with the investigtion of (niteEhorizon optiml ontrol prolemD whih is then not timeEinvrintF his prolem involves qudrti formsF e express it sX for ny positive de(nite or semiEde(nite

u(t) = -C x(t).
× mtrix M nd z EvetorD z, M z = z M z = z 2 M . @IVA
he prolem is @presenttion in more generl form in ulmn @IWTHAAX viner qudrti optiml ontrol prolem qiven the system of iqF @IA nd @PAD with ll system mtries @possily pieewiseA ontinuous funtions of

timeD nd x(t 0 ) = x 0 D nd given the symmetri n × n mtrix A ≥ 0D nd the symmetri @pieewiseA ontinuous n × n mtrix funtion Q(t) ≥ 0 nd m × m mtrix funtion R(t) > 0D (nd the ontrol lwD if it existsD tht minimizes IW the performne index V (x 0 , t 0 , t 1 ; u(•)) = x(t 1 ) 2 A + t 1 t 0 y(t) 2 Q(t) + u(t) 2 R(t) dt . @IWA
hile ulmn nd uoepke @IWSVA nd ulmn @IWTHA expliitly use fellE mn @IWSUA nd disrete dynmi progrmmingD whih ws strong inentive to use the stteEspe representtionD ulmn @IWTHA for the ontinuousE time prolem used the theory of grth¡ eodory @IWQSAF gonverselyD fellmn @IWSUA did not refer to grth¡ eodoryF etD ontinuousEtime dynmi proE grmming is mostly reEdisovery of grth¡ eodory9s theoryD with fellmn9s return funtion plying the sme role s grth¡ eodory9s prinipl funtionF he(ne symmetri mtrix funtion P (t) s the solutionD if it existsD of the mtrix iti eqution @where ll mtries re timeEdependentA

-Ṗ = P F + F P -P GR -1 G P + H
u(t) = -C(t)x(t) , C(t) = R(t) -1 G (t)P (t) , @PIA
QF nd the optiml vlue of the performne index is

V 0 (x 0 , t 0 , t 1 ) = x 0 2 P (t 0 ) .
@PPA e emphsize tht the iti iqF @PHA nd the optiml feedk gin in iqF @PIA re the duls of the iti eqution nd gin ssoited with the ulmn (lter @see setion QFPA QFIFP yptiml regultor @in(niteEhorizonA prolem por the ske of simpliityD we give here only the liner timeEinvrint versionD whih is the only one tht the previous theory ould hndleF ulmn @IWTHA however gve lso the solution for non timeEinvrint prolemF yptiml regultor prolem qiven the timeEinvrint liner system iqF @IA nd @PA with initil stte x(0) = x 0 D nd nonEnegtive de(nite p × p mtrix PI Q nd positive de(nite m × m mtrix RD (nd the ontrolD if it existsD tht minimizes the performne index

V (x 0 ; u(•)) = ∞ 0 y(t) 2 Q + u(t) 2 R dt . @PQA
his study requires the introdution of oth ontrollility nd oservilityF sndeedD in his introdutionD ulmn sttes tht the prinipl ontriution of the pper lies in the introdution nd exploittion of the onepts of ontrollility nd oservilityF etrospetivelyD he ould lso hve quoted the iti iqF @PHAF e use here its lgeri versionD where ll mtries re now onstntX P F + F P -P GR -1 G P + H QH = 0 . @PRA he full theorem isX heorem T

IF sf the pir (F, G) is ompletely ontrollleD then @A the solution P (t) of the iti iqF @PHA hs limit P s t → -∞D whih solves the lgeri iti iqF @PRAD @A the solution of the optiml regultor prolem in stte feedk form is

u(t) = -Cx(t) , C = R -1 G P , @PSA
PP @A nd the optiml vlue of the performne index is

V 0 (x 0 ) = x 0 2 P @PTA @it is tully su0ient tht the pir (F, G) is stilizleD tht isD ∃D : F -GD stleAY
PF sf furthermore the pir (H, F ) is ompletely oservle nd Q > 0D P is positive de(nite nd the system governed y iqF @PSA is symptotilly stle @for the stility resultD it is su0ient tht

Q ≥ 0 nd (Q 1/2 H, F )
is detetleD tht isD ∃L : F -LQ 1/2 H stleAF he dulity we pointed out in the (niteEhorizon prolem holds hereD mking the optiml regultor dul to the sttionry ulmn (lterD tht isD to reliztion of the iener (lterF QFIFQ hisreteEtime se ulmn @IWTHD AD uilding on ulmn nd uoepke @IWSVX eppendixAD whih delt with smpled dt ontrol @whihD s we mentioned eforeD onE sists of onsidering ontinuousEtime systems only t disrete instntsAD of ontinuousEtime system provides the equivlent disreteEtime resultsF he PQ system is mde of iqF @QA nd @RAD nd the performne index isX

V (x 0 , t 0 , t 1 ; {u(•)}) = x(t 1 ) 2 A + t 1 -1 t=t 0 y(t) 2 Q(t) + u(t) 2 R(t)
. @PUA he iti di'erentil eqution is repled y the soElled disrete iti eqution @the system mtries my ll depend on timeA or its lgeri version @with ll system mtries onstntAD where P (t) = P (t + 1) = P X

P (t) = F P (t+1)F -F P (t+1)G(G P (t+1)G + R) -1 G P (t+1)F + H QH, P (t 1 ) = A .
@PVA he optiml feedk ontrol is e strt with the disreteEtime prolemD fter ulmn @IWTHA @where there is neither dded noise in the mesurement eqution nor ontrolAF e onsider disreteEtime liner system exited y white noise nd known ontrol u(•)X

u(t) = -C(t)x(t)
x(t + 1) = F (t)x(t) + G(t)u(t) + D(t)v(t) , x(t 0 ) = x 0 , @QHA y(t) = H(t)x(t) + w(t), @QIA
where (v(t), w(t)) is qussin rndom vrile with zero men nd known ovrineD independent of ll the (v(τ ), w(τ )) for τ = tF sn the simplest seD v(t) is lso independent of w(t)D ut this is not neessry to the theoryF e possile nonEzero rossEorreltion etween v(t) nd w(t) is dul to the presene of ross term x Su in the qudrti performne index of linerE qudrti ontrolF he noise is hrterized y its ovrine mtrix @with δ t,τ the uroneker symolAX

E     v(t) w(t)     (v (τ ) w (τ )) =     V (t) 0 0 W (t)     δ t,τ .
@QPA he initil stte is lso qussin rndom vrile of known distriutionX

E(x(t 0 )) = x0 , E(x(t 0 ) -x0 )(x(t 0 ) -x0 ) = Σ 0 . @QQA

PS

he prolem is to ompute the onditionl mthemtil expettionX

x(t) = E(x(t)|y(τ ), τ < t) . @QRA
he solution is of the form of iqF @ISA initilized t x(t 0 ) = x0 D where the gin K is given through the error ovrine mtrix

Σ(t) = E(x(t) -x(t))(x(t) -
x(t)) D whih is solution of the disrete iti iqF @QSAD whih is dul of iqF @PVAD nd y the formul in iqF @QTAD whih is dul of iqF @PWAX

Σ(t + 1) = F Σ(t)F -F Σ(t)H (HΣ(t)H + W ) -1 HΣ(t)F + DV D , @QSA K(t) = F Σ(t)H (HΣ(t)H + W ) -1 . @QTA
he timeEinvrint in(niteEhorizon se is the internl form of the uolmogorov (lterF gontinuousEtime qiven ontinuousEtime system in internl form exited y white noises @quotes euse suh thing does not exist rigorouslyY it is n engineering presenttionA involved oth in the dynmis nd in meE surementD s in the previous susetion on disrete timeD ut in ontinuous timeD with

E     v(t) w(t)     (v (τ ) w (τ )) =     V (t) 0 0 W (t)     δ(t -τ )
@QUA PT @we use the hir δ s in ulmn nd fuy @IWTIAD ut proilists nowE dys write the system di'erentlyD with di'erent ontologyD in terms of difE fusionsD st¢ o lulusD nd without white noiseAD the onditionl expettion in iqF @QRA we re looking for is the solution of the ontinuousEtime oserver of iqF @ITA with the dul formuls from linerEqudrti ontrolX

Σ = F Σ + ΣF -ΣH W -1 HΣ + DV D , Σ(t 0 ) = Σ 0 . @QVA nd K(t) = Σ(t)H W -1 . @QWA
he timeEinvrint in(niteEhorizon se oinides with the internl represenE ttion of the iener (lterD whihD like the uolmogorov (lterD ws given in externl form @ienerD IWRWAF 3.3 The separation theorem he ontrol lws in iqF @PIAD @PSAD or @PWA re sed on the ext mesureE ment of the stte x(t)F roweverD the underlying ssumption of this theory is tht only the output y(t) is mesuredF he nturl ideD thenD is to sE soite ulmn (lter with the optiml linerEqudrti ontrol lw in n optimum oserverEontroller u(t) = -C x(t)F his ide ws proposed y PU ulmn @IWTHD AD who suggested tht the dulity priniple legitimizes the ssoition of the ulmn (lter estimte with the optiml linerEqudrti ginF re ws then ntiipting the seprtion theoremF e now relte the results pulished in the yer IWTH with the results tht were to ome lterF yptimlity is onsequene of the seprtion theoremF his theorem ws (rst proved y toseph nd ou @IWTIA in disrete timeD with no oservtion noiseD s in ulmn @IWTHAD ndD dullyD no ontrol ostD for nonEqudrti performne indexF st ws proved in ontinuous time y onhm @IWTVAF en erly proofD more spei( to the liner qudrti qussin seD ws lso due to purre @IWTVAF gertinty equivlene resultsD in some prtiulr ses with perfet stte informtion nd without formultion in terms of system theoryD hd ppered erlier in the eonomi literture @imonD IWSTY heilD IWSUAF he ontinuousEtime prolem is muh more di0ultD euse prt of the prolem is to stte the prolem preiselyF his involves ontinuousEtime frownin motionsD st¢ o lulusD (ltrtionsD nd mesurilityF es ulmn in onferene on pplitions of the ulmn (lter in hydrogeologyD hydrulisD nd wter resoures @IWUVA put itX here re three types of (ltersX @iA those whih keep te leves from flling into the te upD @iiA those we re tlking PV out todyD @iiiA those whih re so fny tht only topologists use themF e will not ttempt to present this sttement in modern termsD ut rther stik to the form devised for engineers in the erly sixtiesF yur purpose is to set (rm grounds to devise feedEk dynmi ompenstor @prgrph RFIFQ elowFA viner qudrti qussin @vqA stohsti optiml ontrol proE • the e'ort U i produed y individul i to improve the @ommonA enviE QS ronmentD nd u i = √ U i the deresing mrginl e'et of e'ortD

• the forging e'ort v i of individul iD

• the e0ieny ftor K of ny individul9s forging e'ortD in terms of energy intkeD nd the @positiveA prmeters

• the popultion size nD

• the nturl rte f of dey of the environment qulityD

• the oe0ient g of totl e'ort k u k nd h of totl e'ort k v k in the dynmis governing environmentl qulityD

• the ost q of one unit of forging e'ort nd the ost r of one unit of e'ort devoted to improve the environment per time unitD in terms of energy onsumptionD

• the oe0ients a nd b in the 0ne funtion K @see elowAD nd α = a -q ssumed positiveD

• the oe0ient ρ of the exponentil mortlity rteF QT he model isX

ẋ = -f x + g k u k -h k v k , @RHA K = a + x -b k v k , @RIA
ih individul @eh everA seeks to mximize its performne indexD whih hs the the sme expression for the sesonl or for the lifetime optiE miztionF sn the (rst seD ρ my @optionllyA e tken equl to zero @if the possiility of deth during the seson my e ignoredAY in the seond seD T = ∞ nd the sttionry theory ppliesX

J i = T 0 e -ρt (Kv i -qv i -ru 2 i ) dt .
@RPA e nlyze this prolem in the ppendixF e del with the non homoE geneous hrter y introduing nonEhomogeneous lue funtionD with exponentil disounting y introduing the sme exponentil term in the lue funtionD nd with the mnyEplyer prolem y looking for xsh equiliriumD whih we otin y solving n optimiztion prolemD where eh individul ssumes tht ll other plyers use the sme optiml strtegyF he result is tht oopertion is sustinedF e introdue two oe0ients P nd p @vrile in the (niteEhorizon seD (xed in the in(niteEhorizon seAD QU nd (nd tht the equilirium strtegies reX

u = g r (P x + p) , @RQA v = 1 (n+1)b (1 -2hP )x -2hp + α @RRA
por suitle vlues of the prmetersD nd in prtiulr if n is lrge enoughD iqF @RRA yields positive ontrolsD with positive oe0ients P nd

(t, x) = P (t)x 2 + 2p(t)x + π(t) , @RSA
husD the rtgsf eqution is

-Ẇ + ρW = max u i ,v i 2(P x + p) (-f x(t) + g n k=1 u k (t) -h n k=1 v k ) +v i (α + x(t) -b n k=1 v k (t)) -ru 2 i .
@RTA RT e derive the equtions for the mximizing strtegies u nd v s ru = g(P x + p) , @RUA (n + 1)bv = (1 -2hP )x + α -2hp . @RVA utting these expressions k into the rtgsf iqF @RTA nd equting the oe0ients of the sme powers of xD we otin

-Ṗ = (2n -1) g 2 r + 4n 2 h 2 (n+1) 2 b P 2 -2 f + ρ 2 + (n 2 +1)h (n+1) 2 b P + 1 (n+1) 2 b , -ṗ = -f -ρ + (2n-1) g 2 r P + 4n 2 h 2 P-(n 2 +1)h (n+1) 2 b p + 1-(n 2 +1)hP (n+1)2b α , -π = -ρπ + (2n -1) g 2 r -2hα b + (α+2nhp) 2 (n+1) 2 b .
@RWA hese di'erentil equtions re to e integrted kwrd @numerillyA from the terminl onditions @t T AD whih re equl to zeroF he eqution for P is deoupled from those for p nd for πF yne P is omputedD the eqution for p is linerD nd exp(-ρt)π(t) is otined s n integrlF eFIFP enlysis in (nite horizon @SPA hereforeD x remins positive provided tht α is smller thn limit vlueD whih is esy to (ndF @etullyD if α < 0 ut x + α > 0D lthough p < 0D we still hve u > 0FA eFIFQ enlysis in in(nite horizon sn in(nite horizonD if the lgeri iti iqF @PRA hs rel rootsD then ll oe0ients P D pD nd π re onstnt nd equl to their limit vlues given y the lgeri equtions otined y setting the derivtives to zeroF glosed form solutions re esy to write nd useful for n e0ient numeril impleE menttionD ut otherwise of little interestF ixistene of the positive limit vlue p requires tht 1 -(n 2 + 1)h P > 0D

whih is ondition we hve lredy exminedD nd the twoEinput multivrile systemF he simplest wy to hieve this tsk is with the iti iqF @PRAD revived y ulmn @illemsD IWUIAF he (niteE horizon sesonl prolems were out of reh for this theoryF SI

f + ρ + (n 2 + 1)h (n + 1) 2 b > (2n -

  leftEhnd side digrm of pigure ID or more lssillyD in lok digrmD s in ulmn @IWTQAD nd reprodued in ll textooks sine thenF A is the suE IR spe of rehle ut unoservle sttesY B the suspe of rehle nd oservle sttesY C the suspe of unrehle nd unoservle sttesY D the suspe of unrehle ut oservle sttesF sing sis dpted to g

  the hoie of the gins K nd CF ith hnge of vriles from (x, x) to (x, x)D one n prove the prinE iple of seprtion of the dynmisX heorem Q @eprtion of dynmisA he set of the eigenvlues of the dynmi mtrix of the losed loop oserverEontroller is the union of the eigenvlues of the ontroller F -GC nd those of the oserver F -KHF onhm @IWTUA extended the pole shifting theorem from singleEinput sysE tems @vuenergerD IWTRA to multiEinput ones @erly lterntive proofs re in reymnn @IWTVA nd hvison @IWTVAAX heorem R sf the pir (F, G) is ompletely rehleD thenD given ny moni polynomil p(z) of degree nD there exists mtrix C suh tht the hrterE isti polynomil of F -GC is pF hullyD if the pir (H, F ) is ompletely oservleD the hrteristi polynomil of F -KH n e ssigned to ny desired moni polynomil of degree n y the hoie of KF e dedue n rgument of pure system theory in fvor of the proposed ontrol strutureD nd mens of hoosing C nd K @see prgrph RFIFP elowAF woreoverD in the disreteEtime seD onsequently lso in the smpledEdt IU prolem of ny digitl ontrolD the oserver ehves s oneEstep preditorD euse the system estimtes x(t+1) for x(t+1) with the dt y(τ )D τ ≤ tF he ontrol lgorithm then hs one time step t its disposl to ompute the ontrol u(t + 1) = -C x(t+1)D with the knowledge of the gin CD whih hs een omputed o' lineF 3 Linear quadratic Gaussian (LQG) theory 3.1 Linear quadratic (LQ) optimal control

  and C(t) = (G P (t + 1)G + R) -1 G P (t + 1)F. @PWA foth the (niteE nd the in(niteEhorizon results follow extly s for the ontinuousEtime seD with the sme ontrollility nd oservility onE ditionsF 3.2 The Kalman Filter por the ske of ompletenessD nd to stress dulityD we rie)y review the ulmn (lter @ulmnD IWTHY ulmn nd fuyD IWTIAF ixistene nd stility properties for oth the (nite nd the in(nite horizon ses derive diretly from those of the dul linerEqudrti ontrol prolemF PR QFPFI hisrete time

  lem qiven liner system in internl form with dditive qussin white rnE dom disturnes in the dynmis nd output equtionsD (nd ontrol lw u(.)D if it existsD where u(t) depends only on pst outputs y(τ )D for τ < tD tht minimizes the mthemtil expettion of qudrti performne index mong ll suh ontrol lwsF he nswer is the seprtion nd ertinty equivlene theoremD holding true for disrete time nd ontinuous timeD (niteEhorizon prolemsD nd stE tionry in(niteEhorizon prolemsX heorem U he solution of the liner qudrti qussin stohsti optiE ml ontrol prolem exists nd is otined y repling the stte x(t) y the ulmn (lter estimte x(t) in the liner qudrti deterministi optiml onE trol stte feedk lwF @ositive terms however must e dded to iqF @design in engineering RFIFI viner nd linerized ontrol systems elredy in the preEulmnin erD the liner timeEinvrint theory ws used in vrious physil systemsF ome hd resonle liner physil modelsD whih were timeEinvrint when onsidered t the stedy stteF roweverD most industril systems suh s trnsporttion or energy systemsD demnd nonliner modelF he engineering prtieD thenD is to de(ne desired or nominl output trjetoryD nd tke the error signl s the output of the ontrol systemD tht isD the di'erene etween the tul nd the desired outputsF he ojetive of the ontrol system is then to keep this error signl lose to zero thnks to dynmi ompenstorD whih is dynmi system whose input is the mesured error signlD nd the output the ontrol input of the toEeEontrolled system @tht isD feedk systemD s understood y ienerFA sn order to hieve this golD one uilds liner model s the lineriztion QH of the nonliner model for smll devitions round the nominl trjetoryF his n e done either on the sis of n nlyti nonliner model linerized y (rstEorder expnsionD or of experiments using further prts of the theory @suh s the onsidertion of rossEorreltions etween input nd output pseudoErndom smll devitionsAF feuse the model is n pproximtion of the rel physil systemD the vriles in the model re pproximtions of the physil vrilesF his is the prdox of liner ontrol theoryY the error stys lose to zero euse the ontrol is e0ientY the ontrol is e0ient euse the system is well pproximted y the modelY the system is well pproximted y the model euse the liner pproximtion is goodY the liner pproximtion is good euse the error stys lose to zeroF he errorsD whih re pproximted y the vriles ppering in the ontrol modelD re then to e kept lose to zeroD in spite of disturnes in the dynmisD mesurement errorsD lk of diret mesurement of some key vrilesD not to mention modeling errors nd isesF hese errors nd ises di'er from noises "here not wellE de(ned onept"D nd were t the ineption of roust ontrol theoriesD nd for our purposeD H ∞ Eoptiml ontrol y mes @IWVIAF QI RFIFP yserverEontrollerX lgeri proedure ueeping stedyEstte vrile lose to zero is hieved y foring the system to e su0iently stleF sn tht proessD the poles of the trnsfer funtionD tht isD the eigenvlues of the internl desription of the overll systemD re key informtionF he reson is thtD for ontinuousEtime systemsD the rel prts of these eigenvluesD whih must e negtive to insure stilityD re)et the degree of stilityF he imginry prts mesure the osilltory ehvior of the response of the systemF @sn disrete timeD their modulus must e less thn one to insure stilityFA es we mentioned in susetion PFIFID engineers hd developed dvned tools eyond the mere inspetion of the poles of the trnsfer funtionF hese poles however remin of prmount importneD whih justi(es the use of the theorems on the seprtion of dynmis nd on pole shiftingD s well s the ft tht the prtitioner hooses the poles of the oserver nd the ontrollerD suh tht these poles re independent of one nother @this is the priniple of seprtion of the dynmisAF e rule of thum is tht the oserver must e n order of mgnitude fster thn the ontrollerF rils nd errors with the loliztion of the poles using simultion models @either linerized or nonliner if villeA llow one to onstrut n e0ient ontrol devieF QP RFIFQ udrti synthesis edvnement of modeling sieneD whih ws lrgely due to digitl omputE ersD llowed the tretment of n inresingly high dimension of the modelsD to point tht simple methods sed on the lotion of poles were no longer prtilF por exmpleD the rigid ody dynmis of lnding irplne re deE sried y IPthEorder systemF edding engines nd ontrol surfe dynmis nd )exile modes in modern irliners hs pushed the order wy eyond the IPth oneF woreoverD some prolemsD suh s utomti lnding of n irplne or ontrol of n industril king yle y heting nd ooling n ovenD intrinE silly hve their time horizon (nite nd re not liner time invrintD with emphsis sometimes on terminl error ontrolF hese prolems re eyond the reh of lgeri methodsF ingineers my hve firly preise ide of the origin nd mgnitude of the disturnes ourring in the dynmis nd in the soures of error during mesurementsD let lone in ises inherent in the modelingF his provides sound sis for using the ulmn (lter to ompute n oserver9s ginF es for the ontrollerD the ontrol gin is omputed through qudrti performne index nd the ulmn optiml ginF st is tuned y serhing QQ for stisftory weighting mtries of the performne index y trils nd errorsF his opertion is mde esy nd e0ient if we understnd the gin s the vlue minimizing this performne indexF henD one knows how the di'erent stte nd ontrol vriles will respond to hnge in the weighting mtriesF his proessD known s qudrti synthesisD is wht mde the new theory so populr mong engineersD to the point of eing the min tool used in designing the ontrol systems for the epollo missionF 4.2 Population ecology he use of ulmn (ltering in estimting popultionD suh s in ensusD is nothing newF sn the line of our presenttionD we provide now n exmple of the use of liner qudrti optiml ontrol in rudimentry model of oopertion within n niml popultionF e show how the individul quest of self interest my sustin ollortive ehviorF e use this model to present severl extensions of the liner qudrti theoryD minlyX • the prolem omprises nonEhomogeneous termsD • there is timeEdisounted riterionD • severl plyers re involvedD rther thn single optimizing oneD nd QR there my exist xsh equiliriumF e lrge (nite popultion of nimls live in n re where individuls hve two tions t their disposlX improve the forging environment or forgeF e hve olony of evers onstruting dm in mindD for exmpleF he food intke individuls otin for given forging e'ort is n inresing funtion of environmentl qulityD nd deresing one of totl forging e'ortD whih depletes the resoureF he qulity of the environment deys exponentilly when left unttendedD ut improves when nimls invest in itD with deresing mrginl e'etF st my lso e slightly degrded y the forging e'ortF e investigte the ehvior t the xsh equilirium mong sel(sh inE dividuls either over seson [0, T ] or over lifetimeD with n exponentilly distriuted rndom time of dethF wny models of evolutionry dynmis onverge towrd xsh equiliriumD minly so whenD s hereD this equiliE rium is unique @inent nd frownD PHHSY ndholmD PHIHX hpF TAF e therefore introdue the slr vriles • the qulity x of the environmentD

  (1 -2hP )D s well s stle dynmis nd positive vlues for the qulity xF 5 Carrying on ulmn9s ontriutions in the yers IWTH nd IWTI were powerful stimulus for system theory nd ontrol reserhF wny reserhers followed suitD oth to get further theoretil dvnes @suh s those mentioned in setions QFQ nd RFPFPAD nd to develop lgorithms implementing those theoretil resultsF elgeri system theory ttrted mny reserhers suh s onhm @IWTUA nd remined ulmn9s min reserh reF st ene(ted from dE vned lgeri tools nd ides oming from the theory of utomt @ulmnD IWTSY ulmn et lFD IWTWY ulmnD IWUPAAF he liner qudrti theory of ontrol nd oservtion ws deeply renewed y the theory of H ∞ Eoptiml onE trolD pioneered y mes @IWVIA in the lssil externl frequenyEdomin QV desriptionD nd lter formulted in terms of time domin l ulmn @f sr nd fernhrdD IWWSAD with minimx proilityEfree tretment of unerE tintiesD where the sme dulity shows upD in more omplex setupD nd it mysteriouslyF he si theory quikly ppered in ll textooks of engineering ontrolD nd lter in textooks of mny other domins @eerD PHIIAF he lgorithms were oded into pulily ville softwre pkges suh s wtl or ilD nd they hve een used in wide rnge of pplition dominsD well eyond the industril nd trnsporttion systems of the erly timesD enompssing ll rnhes of engineering s well s nturl nd ioEmedil sienesF viner plus qudrti performne indies were used in the study of niml oopertion @frown nd inentD PHHVAF etD the dynmi linerE qudrti optimiztion theory is mostly sent from the literture on popuE ltion dynmisF wost models used re nonliner @glrk nd wngelD PHHHY storD PHHVAD suh s the logisti growth modelD utD s mentioned in prE grph RFIFID this is lso true in other domins where the linerEqudrti optiE miztion theory hs hd mny pplitionsF e hve skethed here possile use of the linerEqudrti optiml ontrol theory in popultion eologyF QW References f srD F nd fernhrdD F @IWWSAF H ∞ Eyptiml gontrol nd elted winiE mx hesign rolemsX hi'erentil qmes pprohF fostonX firkh¤ userD seond editionF fellmnD F iF @IWSUAF hynmi rogrmmingF rinetonX rineton niE versity ressF firkho' nd w vne @IWTUAF viner elgerF vondonX ew ghelse uE lishingF frownD tF F nd inentD F vF @PHHVAF ivolution of oopertion with shred osts nd ene(tsF roeedings ot the oyl oiety fD PUSX IWVS!IWWRF grth¡ eodoryD gF @IWQSAF ritionsEehnung und rtielle hi'erentilE qleihungen erster yrdnungF veipzigX unerF rnslted into grth¡ eodoryD gF @IWWWAF gluls of ritions nd rtil hi'erenE til iqutions of the pirst yrder @Qrd editionAF xew orkX ghelse ulishing gompnyF glrkD gF F nd wngelD wF @PHHHAF hynmi tte rile wodels in iolE RH ogyD wethods nd epplitionsF yxford eries in iology nd ivolutionF xew orkD eX glrendon ressF hvisonD iF tF @IWTVAF gomments on yn pole ssignment in multiEinput ontrollle liner systemsF siii rnstions on eutomti gontrolD egEIQX URU!URVF purreD F @IWTVAF gommnde optimle stohstique et prinipe de s¡ eprtion de l9estimtion et de l ommndeF evue du girihigD ITX IPW!IQSF reymnnD wF @IWTVAF gomments on yn pole ssignment in multiEinput onE trollle liner systemsF siii rnstions on eutomti gontrolD egEIQX URV!URWF tosephD F hF nd ouD tF F @IWTIAF yn liner ontrol theoryF rnstions of the eFsFiFiF ssD epplitions nd sndustryD VHX IWQ!IWTF ulmnD F iF @IWTHAF gontriutions to the theory of optiml ontrolF foE letin de l oiedd wtemti wexinF roeedings of the symposium on ordinry di'erentil equtionsD wexio gity IWSWD SX IHP!IIWF ulmnD F iF @IWTHAF e new pproh to liner (ltering nd predition RI prolemsF rnstions of the ewi E tournl of fsi ingineeringD VPX QS!RSF ulmnD F iF @IWTHAF yn the generl theory of ontrol systemsF roeedings of the pirst snterntionl gongress of the spegD ID RVI!RWIF ulmnD F iF @IWTPAF gnonil struture of liner dynmil systemsF roE eedings of the xtionl edemy of ienes of the nited ttes of emerE iD RVX SWT!THHF gommunited y F vefshetzD tnury PQD IWTIF ulmnD F iF @IWTQAF wthemtil desription of liner dynmil systemsF sew tournl on gontrolD IX ISP!IWPF ulmnD F iF @IWTSAF elgeri struture of liner dynmil systemsD s" he module of ΣF roeedings of the xtionl edemy of ienes of the eD SRX ISHQ!ISHVF ulmnD F iF @IWUPAF uroneker invrints nd feedkF sn veonrd eiss @edFAD roeedings of the xvEwg gonferene on yrdinry hi'erentil iqutionsD shington hFgFD tune IR!PQ IWUID RSW!RUIF xew orkD vonE donX edemi ressF ulmnD F iF nd fertrmD tF iF @IWTHAF gontrol system nlysis nd RP design y the seond method of vypunovF s"gontinuousEtime systemsF tournl of fsi ingineeringD VP hX QUI!QWQF ulmnD F iF nd fertrmD tF iF @IWTHAF gontrol system nlysis nd deE sign y the seond method of vypunovF ss"hisreteEtime systemsF tournl of fsi ingineeringD VP hX QWR!RHHF ulmnD F iF nd fuyD F F @IWTIAF xew results in liner (ltering nd preE dition theoryF rnstions of the ewi E tournl of fsi ingineeringD IWVQX WS!IHUF ulmnD F iFD plD F vFD nd eriD wF eF @IWTWAF opis in wthemtil ystem heoryF xew orkX wqrwErillF ulmnD F iF nd uoepkeD F F @IWSVAF yptiml synthesis of liner smE pling ontrol systems using generlized performne indexesF rnstions of the ewiD VHX IVPH!IVPTF uolmogorovD endrei @IWRIAF snterpoltion nd ixtrpoltionF fulletin of the edemy of ienesD @hokldy ekdemii xuk A wthemtis SX QEIRF @sn ussinFA RQ vuenergerD hF qF @IWTRAF yserving the stte of liner systemF siii rnstions on wilitry iletronisD VX UR!VHF wgeirouD iF pF @IWUTAF lues nd trtegies for sn(nite ime viner udrti qmesD siii rnstions on eutomti gontrolD egEPID SRU! SSHF werrim sssD gF F @IWSWAF e lss of optimum ontrol systemsF tournl of the prnklin snstituteD PTUX PTU!PVIF xewton trD qF gFD qouldD vF eFD nd uiserD tF pF @IWSUAF enlytil hesign of viner peedk gontrolsF xX tohn iley 8 onsF storD tF @PHHVAF wthemtil iologyF xX ileyEflkwellF oosD gF pF @IWPSAF wthemtil theory of ompetitionF emerin tournl of wthemtisD RTX ITQ!IUSF ndholmD F rF @PHIHAF opultion qmes nd ivolutionry hynmisF gmridge @weAX he ws ressF imonD rF eF @IWSTAF rogrmming under unertinty with qudrti riE terion funtionF ionometriD PRX UR!VIF RR heilD rF @IWSUAF e note on ertinty equivlene in dynmi plnningF ionometriD PSX QRT!QRWF inentD F vF nd frownD tF F @PHHSAF ivolutionry qme heoryD xturl eletionD nd hrwinin hynmisF gmridgeX gmridge niversity ressF eerD F eF @PHIIAF yptiml gontrol with epplitions in ionomisF gmE ridge @weAX he ws ressF ienerD xF @IWRWAF ixtrpoltionD snterpoltion nd moothing of ttionry ime eriesF xew orkX tohn iley nd onsF illemsD tF gF @IWUIAF vest squres sttionry optiml ontrol nd the lgeE ri iti equtionF siii rnstions on eutomti gontrolD eg ITX TPI!TQRF onhmD F wF @IWTUAF yn pole ssignment in multiEinput ontrollle linE er systemsF siii rnstions on eutomti gontrolD egEIPX TTH!TTSF onhmD F wF @IWTVAF yn the seprtion theorem of stohsti ontrolF sew tournl on gontrolD TX QIP!QPTF RS mesD qF @IWVIAF peedk nd optiml sensitivityX model referene trnsE formtionsD multiplitive seminormsD nd pproximte inversesF siii rnstions on eutomti gontrolD eg!PTX QHI!QPHF A A model in population ecology A.1 The Nash equilibrium eFIFI he rmiltonEtoiEgrth¡ eodoryEsssEfellmn @rtgsfA eqution e ssumed tht ll individuls re identilF herefore they shre unique fellmn @or rther sss lueA funtionD whih we seek of the form V (t, x) = exp(-ρt)W (t, x)D with W nonEhomogeneous qudrti funtion of xX W

  ith this elementry exmpleD our min ojetive ws to show IF tht ulmn9s stndrd linerEqudrti theory n e extended in vrious diretionsD nd PF tht some dynmi popultion eology models my led to suh proE lemsF he in(niteEhorizon prolem in oneEplyer version with h = 0 ould e treted with the preEulmnin theory of xewton tr et lF @IWSUAF sn tht se indeedD the v i hve no dynmi e'et nd my e hosen so s to mxE imize the integrndF hus the prolem eomes monovrileF he gme SH version ws more prolemtiD lthough oos @IWPSA otined the gournotE xsh equilirium of liner qudrti di'erentil gme y the stndrd vritionl methods of the lulus of vritionsF he se with h = 0 ws solvle in prinipleD provided tht one performs spetrl ftoriztion for

  nd C re ritrry mtriesF hey ply no roleD euse ξ does

	re equivlent for ontinuousEtime systemsD ut not for disreteEtime onesD
	unless the mtries F (t) re invertile for ll tF gonsequentlyD system is
	ompletely rehle if the pplition u([t 0 , t 1 ]) → x(t 1 )D whih is liner if
	not in)uene yD neither diretly nor through xF hereforeD reliztions of
	di'erent dimensions my represent the sme inputEoutput systemF sn the
	se of iqF @VAD it is trivil ut if hnge of the sis suh s the previous
	one mixes x nd ξD then the dimension in exess my e more di0ult to
	detetF woreoverD other ses re possileF
	he solution of this prolem involves ontrollility nd oservilityF e
	stte is ontrollle if there exists ontrol funtion u(•) tht drives the
	system from this stte to the originF he system is sid to e ompletely
	ontrollle if every stte is ontrollleF
	por the theory to enompss oth ontinuousE nd disreteEtime systemsD
	we now use the lter onept of rehilityF e stte is rehle if there
	exists ontrol funtion tht drives the system from the origin to tht stteF
	e system is ompletely rehle if every stte is rehleF he two onepts
	II

  he eqution for P is iti eqution Ṗ = -Π(P )D where Π is seondE degree polynomil with two positive roots if these roots re relF he disE riminnt of Π is lwys positive if side goes to zero s n inresesA nd positive otherwise if h is lrger thn limit vlue tht dereses s n inresesF et terminl time T D P rehes H with negtive slopeY s t → -∞D P (t) onverges towrd the smllest @relA root P of Π if it existsF st diverges to +∞ if notF @he possiility of P divergingD ontrry to wht is stted in heorem TD is due to ross term xv in the riterionFA herefore P (t) remins positiveD nd smller thn P when it existsF P goes to zero s 1/n 2 when n inreses to in(nityF por h lrger thn limit h @whih goes to zero s n goes to ∞AD the Π hs rel rootsD nd A 1 -(n 2 + 1)hP nd fortiori 1 -2hP re positive for ll tF he eqution for -ṗ in iqF @RWA ensures tht p is lso positive for ll t nd smller thn limit vlue pD whih goes to zero s 1/n when n → ∞F reneD for n lrge enoughD if x is positiveD so re u nd v F pinllyD the dynmis of environment qulity under the equilirium strteE

	gies reX @the rightEhnd quntity ẋ = -f -nh (n + 1)b	+	ng 2 r	+	2nh 2 (n + 1)b	P x+		ng 2 r	+	2nh 2 (n + 1)b	p-	nhα (n + 1)b	.
	Π	1 (n 2 + 1)h	=	(2n -1)g 2 (n 2 + 1) 2 r	×	1 h 2 -2	f + ρ 2 n 2 + 1	×	1 h	-	1 (n + 1) 2 b	@SIA
	f + is negtiveD ensuring tht A RV ρ 2 ≥ (2n -1)g 2 (n + 1) 2 br RU	1 2	,	@SHA

  is stis(ed provided tht n is lrge enoughF st lso ensures tht e -ρt x 2 is symptotilly stleD nd tht this stility is preE served if some of the plyers ply u i = 0F foth re neessry onditions for the xsh equilirium in in(nite horizon to exist @wgeirouD IWUTAF A.2 Time domain versus frequency domain

	he dynmis of environment qulity re stle provided tht
	f +	nh (n + 1)b	>	ng 2 r	+	2nh 2 (n + 1)b	P ,	@SRA
	gin ondition whih
					1)	g 2 r	+	4n 2 h 2 (n + 1) 2 b	P ,	@SQA
	gin ondition whih is stis(edD provided n is lrgeF
				RW