Okan Yılmaz
email: okan.yilmaz@epfl.ch

Jérémy Bleyer

Jean-François Molinari

Influence of heterogeneities on crack propagation

Keywords: pre-strained plate, dynamic fracture, cohesive elements, crack propagation, branching

The influence of material heterogeneities is studied in the context of dynamic failure. We consider a pre-strained plate problem, the homogeneous case of which has been widely studied both experimentally and numerically. This setup is used to isolate the effects of the elastic field resulting from pre-straining and stress wave interactions throughout the crack propagation by adding stiffer and denser regions in the plate. While the crack tip is pushed away by stiffer inclusions, it is attracted to the denser ones. With the presence of denser media, only a portion of the total elastic energy in the system is effectively used to drive crack propagation, leading to a drop in the velocity of its tip in comparison to the homogeneous case. The crack velocity is shown to be an invalid criterion for crack branching. Instead, we introduce an effective stored energy to analyze the crack velocity and the emergence of crack branching instabilities.

Introduction

Fracture prevention of materials and structures plays a key role in modern engineering designs. Understanding the physical mechanisms governing the dynamics of fracture in brittle and quasi-brittle media is thus crucial to improve existing procedures and design better products. However, this is still a challenge due to the presence of material heterogeneities at different length scales, complexity of fracture process zone, and dynamic stress redistribution caused by the emitted stress waves from the moving crack tip and boundary reflections [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. Although elastodynamic linear fracture mechanics theory provides the time variation of the dynamic stress intensity factor as well as other field quantities, it is not adequate to capture the dynamic response of brittle materials [START_REF] Ravi-Chandar | Dynamic fracture[END_REF]. Experiments indicate several phenomena including limiting crack velocities below the Rayleigh wave speed of the material and crack branching at micro and macroscale that the linear theory fails to describe.

Especially for the evolution of crack dynamics, much of the theoretical and experimental work is done in homogeneous settings. Heterogeneity is present in materials at different scales, which are likely to introduce modifications to the homogeneous medium solutions [START_REF] Sagy | Dynamic fracture of granular material under quasi-static loading[END_REF][START_REF] Albertini | Supershear transition of dynamic mode II fracture in heterogeneous elastic media[END_REF]. Material inhomogeneities have been studied at all scales from kilometers [START_REF] Fisher | Statistics of earthquakes in simple models of heterogeneous faults[END_REF][START_REF] Ben-Zion | Collective behavior of earthquakes and faults: Continuumdiscrete transitions, progressive evolutionary changes, and different dynamic regimes[END_REF] to nanometers [START_REF] Zhang | Wavy cleavage fracture of bulk metallic glass[END_REF][START_REF] Murali | Atomic scale fluctuations govern brittle fracture and cavitation behavior in metallic glasses[END_REF], the applications of which range from earthquake faults to the nanocomposition of materials.

These various alterations to the dynamic process of fracturing are governed by the changes in the elastic field and stress wave emissions and reflections in the medium. Here, we focus on the evolution and propagation of fractures and the resulting crack patterns.

The objective of this paper is to examine the influence of distant heterogeneities on the dynamics of fracture. For this purpose, we numerically analyze the well-studied pre-strained plate configuration [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF][START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF][START_REF] Livne | Universality and hysteretic dynamics in rapid fracture[END_REF][START_REF] Goldman | Acquisition of inertia by a moving crack[END_REF][START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF][START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. Stiffer and denser regions are introduced within this rapid fracture scheme to study their effect on the crack propagation velocity and crack branching instabilities.

In Section 2, the cohesive-element methodology, which is used to model fracture is explained. Section 3 presents a pre-strained rectangular plate problem using a homogeneous medium and the influence of the heterogeneities are studied in Section 4. Stiffer inclusions are added to alter the elastic field in Section 4.1. The strip is narrowed with denser regions in Section 4.2 to investigate the influence of wave interactions while preserving the homogeneous elastic properties. Wave interactions are further studied with circular dense inclusions in Section 4.3, allowing comparisons with the cases having stiffer inclusions.

Conclusions are given in Section 5.

Fracture model

Dynamic fracture is modeled using the cohesive-element methodology. The method is based on the cohesive crack model of Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF].

To represent the decohesion, we bury all complex debonding processes in a phenomenological, simple, cohesive law relating the traction and opening displacement. This constitutive response is called traction-separation law:

T = T(∆ ∆ ∆), ( 1
)
where T is the traction acting on the separating surfaces and ∆ ∆ ∆ is the relative opening displacement vector. We use a simple linear irreversible softening law [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF][START_REF] Ortiz | Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis[END_REF]. The existing free potential energy φ, is assumed to depend only on one effective scalar displacement δ of the following form:

δ = ∆ 2 n + β 2 ∆ 2 s , (2)
where ∆ n and ∆ s are the normal and tangential displacement components while β is the parameter that couples these two displacements. This parameter can be estimated by imposing lateral confinement on specimens subjected to high strain rate axial compression [START_REF] Chen | Static and dynamic compressive behavior of aluminum nitride under moderate confinement[END_REF] and roughly defines the ratio of the critical stress intensity factors of Mode II to Mode I of the material [START_REF] Ruiz | Three-dimensional finite-element simulation of the dynamic brazilian tests on concrete cylinders[END_REF]. The derivation of free potential energy with respect to the opening displacement gives the cohesive tractions:

T = ∂φ ∂∆ ∆ ∆ = T δ ∆ n n + β 2 ∆ s s , (3)
where n and s are the unit vectors in normal and tangential directions. Crack opening condition is denoted with δ = δ max and δ > 0, and crack closure and reopening case is stated as δ < δ max . Thus, the traction is

T =      f ct 1 -δ δc for δ = δ max and δ > 0 δ δmax T max for δ < δ max , (4)
where f ct is the critical stress associated to the cohesive element, δ c is the effective relative displacement value beyond which complete decohesion occurs, δ max is the maximum value of the achieved effective displacement up to that instant and T max is the associated traction value for maximum relative displacement.

The cohesive law and the areas that are associated to the energy dissipation are illustrated in Figure 1. Cohesive elements are dynamically inserted at the inter-element boundaries when and where the local stress exceeds a critical stress threshold, which is the tensile strength of the corresponding material.

The irreversible part of dissipated energy is shown with the dark grey area and denoted by E irrev while the reversible part is the light grey area named E rev . Damaging of the cohesive element is an irreversible process while tractionseperation law allows recovering energy upon unloading. Damage in a cohesive element is denoted as:

D = min δ max δ c , 1 . (5)
A simple penalty model is incorporated in the cohesive law for the case of compression:

T n = α p ∆ n if ∆ n < 0, (6)
where α p is the stiffness parameter and T n is the normal component of the traction vector.

The total dissipated fracture energy during the separation process is given by the entire grey region. When the dissipated energy is equal to the total area under the curve shown in the figure, which is 1 2 f ct δ c , full decohesion occurs and two free surfaces are formed at the interelement boundaries.

Pre-strained rectangular plate

The pre-strained rectangular plate problem was studied experimentally by [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] using PMMA and by [START_REF] Livne | Universality and hysteretic dynamics in rapid fracture[END_REF][START_REF] Livne | Breakdown of linear elastic fracture mechanics near the tip of a rapid crack[END_REF][START_REF] Livne | The near-tip fields of fast cracks[END_REF][START_REF] Goldman | Acquisition of inertia by a moving crack[END_REF] using polyacrylamide gels. Several numerical investigations can be found using different fracture modeling approaches such as cohesive elements [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF], non-local damage modeling [START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF], and phase-field modeling [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. The model experiment is conducted as follows. First, the plate is pre-strained by moving the upper and lower boundaries to the desired displacement. In this strained state, a sharp crack is created in the middle of one side of the sample using a razor edge. Pre-notched version of the identical setup is also popular and widely studied both experimentally [START_REF] Ramulu | Mechanics of crack curving and branching -a dynamic fracture analysis[END_REF][START_REF] Fineberg | Instability in dynamic fracture[END_REF][START_REF] Sharon | Local crack branching as a mechanism for instability in dynamic fracture[END_REF][START_REF] Sharon | Confirming the continuum theory of dynamic brittle fracture for fast cracks[END_REF] and numerically [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Song | A comparative study on finite element methods for dynamic fracture[END_REF][START_REF] Ha | Studies of dynamic crack propagation and crack branching with peridynamics[END_REF].

δ max δ c T max f ct δ T E rev E irrev α p 1

60

The only difference with the pre-strained condition is that the crack is seeded in the specimen before applying the load. For the plate geometry, pre-straining conditions, and crack size, see Figure 2.

PMMA is used as the plate material, see Table 1 for the material parameters.

The coupling parameter β for the cohesive-element method is chosen to be 1 for this uniaxial tension case. Using the finite-element method, we solve the equation of motion

Mü + f int = f ext , ( 7
)
where M is the mass matrix, u is the displacement vector, and f int and f ext are the internal and external force vectors respectively. This equation is solved at every time step using an explicit time integration scheme and the displacement field is updated. For further details regarding the implementation see [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. The finite-element code used for the study is based on a Newmark explicit time integration scheme, which is conditionally stable when the time step is less than a critical value set by the Courant-Friedrichs-Lewy condition [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]:

∆t crit = α min e=1...N l e c , (8)
where c is the pressure wave speed of the associated material, l e is the size of element e, and N is the total number of elements. The smallest value of the calculated time step over all elements is the chosen time step for the simulation.

A security coefficient, α, is also used (0.2 for this study) to ensure stability.

A two-dimensional model with plane-stress assumption is used. Rayleigh wave speed of the material is calculated using the following formula:

c R = 0.862 + 1.14ν 1 + ν c s , (9)
where c s is the shear wave speed, c s = µ/ρ. For PMMA, c s = 984.8 m/s and c R = 919.9 m/s.

65

Cohesive-zone modeling introduces another length scale set by the elastic and fracture properties. This length scale can be expressed for simple cases (mode I loading, homogeneous material) as [START_REF] Palmer | The growth of slip surfaces in the progressive failure of over-consolidated clay[END_REF][START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]:

l z = EG c f 2 ct . ( 10
)
The condition l m < l z between the length scales has to be satisfied in our simulations where l m is the biggest element size. The cohesive zone must contain several elements (typically around four [START_REF] Snozzi | Influence of the meso-structure in dynamic failure simulation of concrete under tensile loading[END_REF]) and should be small compared to the sample size for mesh independency. The cohesive zone length for PMMA material is calculated as 0.23 mm. A mesh of average mesh size l av = 0.05 mm is used to correctly resolve the cohesive zone in front of the crack tip. Fracture energy, G c (J/m 2 ) 300

Tensile strength, f ct (MPa) 75

Coupling parameter, β 1

Figure 3 illustrates the crack patterns at displacements from ∆u = 0.035 mm to ∆u = 0.10 mm. At lower displacements, there exists only a single crack propagating from the initial cut. However, increasing the level of pre-straining, we observe first microbranching and an increase in the damage band. Microbranches are defined as the local branching events where main crack sprouts side branches that propagate for a short distance and die [START_REF] Sharon | Local crack branching as a mechanism for instability in dynamic fracture[END_REF]. This is a numerical instability observed beyond a critical velocity, which [START_REF] Fineberg | Instability in the propagation of fast cracks[END_REF] observed around 0.4c R . The limiting velocity is in fact a result of this instability as it prevents the crack to further accelerate while energy dissipation is not limited to a single propagating crack tip [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. In addition, local branching events turn into mac-robranches, which do not die right away and can propagate in the medium for some time when ∆u ≥ 0.06 mm.

The cohesive-element approach is known to exhibit a mesh dependency on the dissipated energy, thus the amount of microcracking, especially for coarse meshes [START_REF] Molinari | The cohesive element approach to dynamic fragmentation: the question of energy convergence[END_REF]. In this study, we use meshes that are fine enough to fully resolve the cohesive zone in front of the crack tip with the help of parallel computing.

However, a convergence analysis of the dissipated energy with the mesh size is missing. Nonetheless, salient experimental features of the dynamic crack propagation scheme such as the limiting crack velocity and branching at higher pre-straining levels [9, 10] may be captured using cohesive elements.

The snapshots in Figure 3 are taken when the crack fully propagates from one end to the other. The crack needs less time to reach the right end of the plate at higher pre-straining levels. Crack propagation velocities as a function of x-position of the crack tip are plotted in Figure 4. The instantaneous crack velocity is calculated by fitting a linear function for several consecutive time steps versus crack tip position around that particular location.

The velocities are higher than the experimentally observed limiting velocities. It was shown by [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] that the rate-independent cohesive law wrongly predicts Rayleigh wave speed of the material, c R as the limiting velocity, which is reproduced by our model. However, experiments show that the limiting velocities are in the range of 0.5 -0.65c R for glass, 0.58 -0.62c R for PMMA, and 0.3 -0.45v/c R for Homalite-100, see [START_REF] Ravi-Chandar | Dynamic fracture[END_REF] and references therein. This velocity also depends on the experimental setup for a given material. Zhou et al. [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] also showed that the experimentally observed limiting velocities can be obtained with a phenomenological rate-dependent cohesive law. Since this study focuses on the effects of heterogeneities, the rate-independent scheme is preserved.

Crack propagation in this setting can be altered by two mechanisms. First is the changes in the elastic field, which can be obtained by introducing inclusions of different elastic properties. Wave dynamics can also influence crack propagation. To have different stress wave interaction in an identical elastic field, regions with same elastic modulus but different mass density will be introduced.

Influence of heterogeneities

Stiff inclusions

The pre-strained elastic field is changed by adding a series of circular inclusions, see Figure 5. Elastic modulus of inclusions is ten times the plate material, E inc = 10E 0 = 30.9 GPa. Boundary conditions are identical to the homogeneous case. l x is the x-distance from the initial crack tip to the centerline of the first inclusion, l y is the y-distance from the plate centerline to the bottom tangent of the inclusions and l c is the distance between centers of the inclusions.

The first inclusion is l x = 3 mm away from the initial crack tip and the distance between the inclusions is l c = 3 mm. These two dimensions are kept constant while the inclusion size and distance to the plate centerline are varied.

Figure 6 shows crack propagation for different inclusion size and centerline proximity. The displacement is ∆u = 0.035 mm is applied for all cases to have a single propagating crack.

We observe that the crack is pushed away by the stiff inclusions. This deflection increases when the diameters of the inclusions are larger. Total deflections Pushing the inclusions closer to the centerline (l y from 0.5 to 0.1) does not have a prominent influence to the total deflection, they are measured in this case 0.5 mm and 0.7 mm respectively. No significant change in crack velocity is observed with the stiff inclusions.

Dense strip

Crack propagation may also be influenced by the elastic wave interactions during its movement. The pre-strained plate is a suitable setup to control stress wave dynamics as the boundaries are fixed at a given displacement and only the propagating crack tip emits stress waves. Here, we introduce regions ten times denser than PMMA material, ρ d = 10ρ 0 in the upper and lower part of plate keeping the elastic modulus constant. Thus, the initial elastic field is identical to the homogeneous case and only wave interactions can modify the crack. The

calculated for the given material values as

R = c d ρ d c 0 ρ 0 = ρ d ρ 0 E d ρ 0 ρ d E 0 = ρ d ρ 0 = √ 10. ( 11
)
Assuming one-dimensional stress-wave propagation of magnitude f I , reflected and transmitted waves are found using the impedance ratio:

f R = f I R -1 R + 1 = 0.52f I , f T = f I 2R R + 1 = 1.52f I . (12)
At the boundary of dense region, a compressive wave is reflected back as a compressive wave of half the magnitude and transmitted to the dense region with 145 a higher magnitude. Transmitted compressive wave is then reflected back from the plate boundary as a tensile wave, however travels at a lower speed since the medium is denser. Entering back to the PMMA region, the impedance ratio is reversed. With R = 1/ √ 10, the reflection coefficient is -0.52 and transmission coefficient is 0.48. Hence, approximately 73% of the compressive wave emitted 150 returns back as a tensile wave when it travels to the simulation-box boundary and back. Note that this is valid for only one wave in a one-dimensional setting, thus with continuous wave emission from a fast moving crack tip, the wave interaction picture is much more complex.

Figure 8 shows crack velocities for various PMMA widths (b parameter).

155

The plots are grouped under the levels of pre-straining. In all cases we observe that crack velocity decreases as b gets narrower. Apart from this global decrease, local drops are observed. This setup is closer to the infinite strip problem [START_REF] Marder | New dynamical equation for cracks[END_REF][START_REF] Marder | Adiabatic equation for cracks[END_REF] than the infinite plate problem. As we decrease the width of the PMMA material with denser media, we can draw analogies with infinite strip case. The system is bounded vertically in this case hence the interactions with external boundaries are not neglected. Once time t significantly exceeds 2b/c s , which is the time that an emitted wave travels from the crack tip and back, crack velocity approaches a steady state value, which can be significantly lower than c R [START_REF] Goldman | Acquisition of inertia by a moving crack[END_REF]. The crack tip can emit both longitudinal and shear waves thus the time scale of interest is defined by the slower wave speed to remain in the conservative side. In our configuration, the proposed time scale 2b/c s is much lower than the crack propagation time from one end to the other, which explains the decrease in the steady-state velocities as we decrease b.

Figure 9 shows the crack configurations at different dense region configurations at ∆u = 0.06 mm. Despite the slowing of the crack tip, branching occurs at b = 1 mm at a later stage of propagation. We observe that local branching events can be suppressed by decreasing b (which decreases crack tip velocity), see the homogeneous case in Figure 3d for comparison. However, at extremely low b values, macrobranching occurs at a velocity much lower than c R . Hence, the crack velocity alone does not seem to be a valid criterion for branching phenomena. The efforts to relate the crack velocity to its branching dates back to the analysis of [START_REF] Yoffe | The moving griffith crack[END_REF] with the argument that a critical velocity of 0.6c s causes the crack tip to branch at an angle 60 • to the initial direction of crack propagation. Instead, crack branching can be viewed as a consequence of an excess of available energy flowing to the crack tip, which cannot be dissipated by a single crack propagation, suggesting the existence of a critical value of the stress intensity factor or energy release rate [START_REF] Seelig | On the interaction and branching of fast running cracks -a numerical investigation[END_REF][START_REF] Ravi-Chandar | Dynamic fracture[END_REF][START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. Furthermore, the experiments of [START_REF] Ravi-Chandar | An experimental investigation into dynamic fracture: III. on steady-state crack propagation and crack branching[END_REF] show that the cracks propagate at an identical speed after and prior to branching event, justifying that the crack velocity by itself cannot be a meaningful criterion.

The initially stored energy is a well-defined property of the system and can be related to dynamic parameters such as the crack tip velocity or fracture energy dissipated by crack advance [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF]. For the pre-strained plate configuration, total stored energy per specimen width is calculated as

W = 2E(∆u) 2 h . (13)
Figure 10a shows the crack velocities versus the total stored energy. Filled points indicate the homogeneous plate cases while hollow points are the cases with b < 16 mm. Although the set of points corresponding to a particular 190 pre-straining level indicate different crack velocities, total stored energy is identical. Introducing dense regions do not change the total energy available in the system as it is only a function of stiffness (constant in all cases) and applied displacement.

One can attain much lower velocities when the propagating crack is surrounded with denser regions. and this is proportional to their amount. In other words, These regions limit the elastic energy, which can be used to propagate the crack tip. To explain the difference in the failure dynamics, we introduce an effective stored energy, which is the multiplication of total stored energy by b/h, the ratio of actual PMMA volume in the plateto the denser material:

W eff = 2E(∆u) 2 h b h . (14)
Using the effective stored energy instead of the total one, the points in Fig-

ure 10a fall remarkably on the same curve, see Figure 10b. Thus, the given density ratio is able to trap the corresponding portion of the mechanical potential energy of the system to the denser regions where mechanical waves travel slower, therefore limiting the amount to be used for crack propagation. Dividing the total dissipated energy by the distance from the initial crack 200 tip to the end of the plate (28 mm), we get the specific fracture energy, which is plotted in Figure 11 for different b and ∆u with respect to the crack velocity and effective stored energy. The outlier point is the case ∆u = 0.06 mm, b = 1 mm, at which crack branching is observed at a low velocity, see Figure 9d for the pattern. The material property G c constitutes a lower bound for the normalized dissipated energy in the system. Above this limit, local branching events start, which is found to be around 650 m/s. With the increasing energy dissipation capability of the system, we observe macrobranches, which can also occur at low velocities. (a) A lower bound of 300 J/m 2 , which is the material property of the traction-separation law is shown with a horizontal dashed line. The vertical dashed line is the Rayleigh wave speed of the material. (b) Both the horizontal and vertical dashed lines denote the limit Gc = 300 J/m 2 . It is possible to have crack propagation in the grey area (W eff < Gc) as the total stored energy is higher than the effective one.

Figure 11b shows the relation between the dissipated and effective stored energies. In this case, both the horizontal and vertical limits denote the material fracture energy, G c = 300 J/m 2 . An explanation for the initiation of the local branching events is provided by the effective stored energy. We observe that when W eff > G c , microbranching starts. It is possible to have crack propagation when W eff < G c (the grey-colored area), since the total available elastic energy in the plate is higher as shown before in Figure 10a.

Dense inclusions

Lastly, we consider circular shaped dense inclusions for direct comparison with stiff inclusions. The scheme in Figure 5 is used but with identical elastic modulus and density of inclusions ρ inc = 10ρ 0 . Only l y = 0.1 mm case is compared. The applied displacement is ∆u = 0.035 mm to prevent any branching. The attraction of the crack tip can be explained by the stress wave interactions creating energy favored directions. Figure 13 shows the evolution of strain energy density.

We observe that the strain energy density is higher at the upper portion of the crack tip vicinity as it propagates. The crack tip can be located in each snapshot by tracking the localized region of strain energy density. The deflection in the energy-favorable direction is expected as the Young's modulus is constant throughout the plate. The contribution of the kinetic energy is found to be insignificant compared to the strain energy.

The velocity of the crack tip is plotted in Figure 14. When the crack is

235

passing by an inclusion it slows down and this effect is also proportional to the inclusion size. The crack tip is attracted to the dense inclusions and it passes through them during the analysis. The drop in the velocity is considerably higher while the crack propagates in the denser inclusion as the material wave speed is ρ inc /ρ 0 times lower.

Conclusions

We have numerically analyzed the pre-strained plate configuration to study the influence of distant material heterogeneities on dynamic failure. PMMA is used as the plate material. With this setup, the effects of elastic field resulting from pre-straining and stress wave interactions throughout the crack propagation are isolated by defining stiffer and denser media in the plate. The crack tip is repelled by the stiffer inclusions, while it is attracted to the identical-shaped denser ones.

Lower steady-state crack velocities are obtained when the PMMA strip is narrowed down by defining denser media at the top and bottom of the plate.

250

We have shown that a portion of the total stored elastic energy can be trapped in the dense regions, thereby limiting the energy flow to the crack tip. The portion that is still effective in the process of crack propagation is calculated by multiplying the total stored energy by the volume fraction of actual PMMA material. Using the effective stored energy instead of the total elastic one, the crack velocities at different applied displacements collapse on the same curve.

We observed crack branching at much lower velocities than the limiting velocity of our system, providing further evidence that the crack velocity is not a valid criterion for the branching phenomena.

This study limits itself to the theoretical considerations of distant heterogeneities. There are engineering materials such as metals, ceramics, polymers and elastomers, the combinations of which can give the discussed stiffness and density ratios by various designs, see an Ashby plot of Young's modulus and density, for example [START_REF] Fleck | Micro-architectured materials: past, present and future[END_REF]. However, the recommendations for possible material combinations are outside the scope of this work as this requires a thorough analysis of considered materials, design and production phases.

Figure 1 :

 1 Figure 1: Graphical representation of linear irreversible cohesive law [16, 17]. Arrows denote loading and unloading directions.

Figure 2 :

 2 Figure 2: Pre-strained PMMA plate problem with an applied displacement ∆u to the top and bottom boundaries.

Figure 3 :

 3 Figure 3: Crack propagation with different pre-strain loadings. The snapshots are taken at times (a) t = 4.53 × 10 -5 s (b) t = 4.30 × 10 -5 s (c) t = 4.07 × 10 -5 s (d) t = 3.85 × 10 -5 s (e) t = 3.62 × 10 -5 s (f) t = 3.40 × 10 -5 s. The applied displacements are (a) ∆u = 0.035 mm (b) ∆u = 0.04 mm (c) ∆u = 0.05 mm (d) ∆u = 0.06 mm (e) ∆u = 0.08 mm (f) ∆u = 0.10 mm.

Figure 4 :

 4 Figure 4: Crack velocity normalized with the Rayleigh wave speed of the material as a function of crack tip position for different levels of pre-straining the homogeneous PMMA plate given in Figure 2.

Figure 5 :

 5 Figure 5: Pre-strained PMMA plate problem with the addition of stiff inclusions of Young's modulus E inc = 10E 0 and diamater d inc . The density of the inclusions is identical to PMMA. lx is the x-distance from the initial crack tip to the centerline of the first inclusion, ly is the y-distance from the plate centerline to the bottom tangent of the inclusions and lc is the distance between centers of the inclusions. The applied displacement is ∆u and boundary conditions are identical to Figure 2.

Figure 6 :Figure 7 :

 67 Figure 6: Crack propagation in different inclusion pattern configurations with the inclusions having the same density but different Young's modulus E inc = 10E 0 . The applied displacement is ∆u = 0.035 mm for all cases. The snapshots are taken at times (a) t = 4.67 × 10 -5 s (b) t = 4.67 × 10 -5 s (c) t = 4.55 × 10 -5 s (d) t = 4.55 × 10 -5 s. The diameters and distance from the centerline are (a) d inc = 1 mm, ly = 0.5 mm (b) d inc = 1 mm, ly = 0.1 mm (c) d inc = 2 mm, ly = 0.5 mm (d) d inc = 2 mm, ly = 0.1 mm. The first inclusion is lx = 3 mm away from the initial crack tip and the distance between the inclusions is lc = 3 mm.

Figure 8 :

 8 Figure 8: Crack velocities as a function of crack tip position for different dense region configurations, b-widths (mm) in the legends. The applied displacements are (a) ∆u = 0.035 mm (b) ∆u = 0.04 mm (c) ∆u = 0.05 mm (d) ∆u = 0.06 mm. The horizontal dashed line is the Rayleigh wave speed of the material. In all cases it is observed that decreasing b, the crack velocity decreases. Apart from this global decrease, local drops are also observed, see cases b = 8 and 10 mm in (a).

Figure 9 :

 9 Figure 9: Crack propagation in different dense plate configurations. Grey areas denote the PMMA material and black areas are for the denser region with ρ d = 10ρ 0 . The applied displacement is ∆u = 0.06 mm for all cases. The snapshots are taken at times (a) t = 4.01 × 10 -5 s (b) t = 4.32 × 10 -5 s (c) t = 6.91 × 10 -5 s (d) t = 1.03 × 10 -4 s. The width of the confined PMMA regions are (a) b = 10 mm (b) b = 5 mm (c) b = 3 mm (d) b = 1 mm.

Figure 10 :

 10 Figure 10: Stored energy per unit length W (a) as a function of crack velocity. Legend shows the applied displacement. Filled circles represent the homogeneous cases while hollow circles are the cases where b < 16 mm. Using the effective stored energy per unit length W eff (b) we observe that the points fall on the same curve. The vertical dashed line is the Rayleigh wave speed of the material.

Figure 11 :

 11 Figure 11: Total dissipated energy normalized with the distance from the initial crack tip to the plate end (0.028 m) versus the crack velocity (a) and effective stored energy (b). Filled circles represent the homogeneous cases while hollow circles are the cases where b < 16 mm.

Figure 12

 12 Figure 12 illustrates the crack propagation for two different inclusion sizes. Unlike the original stiffer inclusions case (see Section 4.1), the crack tip is attracted to the denser inclusions and this attraction is proportional to the inclusion size.

Figure 12 :

 12 Figure 12: Crack propagation in different inclusion pattern configurations with the inclusions having the same Young's modulus but a density of ρ inc = 10ρ 0 . The snapshots are taken at times (a) t = 6.74 × 10 -5 s (b) t = 9.93 × 10 -5 s. The diameters are (a) d inc = 1 mm (b) d inc = 2 mm and the distance to the centerline is ly = 0.1 mm in both cases. The first inclusion is lx = 3 mm away from the initial crack tip and the distance between the inclusions is lc = 3 mm.

Figure 13 :

 13 Figure 13: Strain energy density N/mm 2 at times (a) t = 0 s (b) t = 1.75 × 10 -5 s (c) t = 2.34 × 10 -5 s (d) t = 4.68 × 10 -5 s with the inclusions having the same Young's modulus but a density of ρ inc = 10ρ 0 . The diameters are d inc = 2 mm and the distance to the centerline is ly = 0.1 mm. The first inclusion is lx = 3 mm away from the initial crack tip and the distance between the inclusions is lc = 3 mm (the configuration shown in Figure 12b).

Figure 14 :

 14 Figure 14: Crack velocity as a function of crack tip position for different inclusion pattern configurations with the inclusions having the same Young's modulus but a density of ρ inc = 10ρ 0 . Inclusion positions are colored with grey. The diameters are (a) d inc = 1 mm (b) d inc = 2 mm and the distance to the centerline is ly = 0.1 mm in both cases. The first inclusion is lx = 3 mm away from the initial crack tip and the distance between the inclusions is lc = 3 mm. The horizontal dashed line is the Rayleigh wave speed of the material. The grey areas are the x-locations of the inclusions.

Table 1 :

 1

Material properties of PMMA [10, 1]. PMMA Young's modulus, E (GPa) 3.09 Poisson's ratio, ν 0.35 Density, ρ (kg/m 3 ) 1180


Haris Haralambous
email: h.haralambous@frederick.ac.cy

Antonis Ioannou

Harris Papadopoulos
email: h.papadopoulos@frederick.ac.cy

A Neural Network tool for the interpolation of foF2 data in the presence of sporadic E layer

Keywords: Ionosphere, HF communications, F2 layer critical frequency

This paper presents the application of Neural Networks for the interpolation of (critical frequency) foF2 data over Cyprus in the presence of sporadic E layer which is a frequent phenomenon during summer months causing inevitable gaps in the foF2 data series. This ionospheric characteristic (foF2) constitutes the most important parameter in HF (High Frequency) communications since it is used to derive the optimum operating frequency in HF links and therefore interpolating missing data is very important in preserving the data series which is used in long-term prediction procedures and models.

Introduction

The ionosphere is defined as a region of the earth's upper atmosphere where sufficient ionisation can exist to reflect radio waves in the frequency range 1 to 30 MHz therefore facilitating long distance HF (High Frequency) communication by skywave propagation. It ranges in height above the surface of the earth from approximately 50 km to 1000 km. The influence of this region on radio waves is accredited to the presence of free electrons.

The uppermost layer of the ionosphere is the F2 region which is the principal reflecting region for long distance HF communications [START_REF] Goodman | HF Communications, Science and Technology[END_REF][START_REF] Maslin | The HF Communications, a Systems Approach[END_REF][START_REF] Mcnamara | Grid The Ionosphere: Communications, Surveillance, and Direction Finding[END_REF]. The maximum frequency that can be reflected at vertical incidence by this layer is termed the F2 layer critical frequency (foF2) and is directly related to the maximum electron density of the layer (see figure 1). The F2 layer critical frequency is the most important parameter in HF communication links since when multiplied by a factor which is a function of the link distance, it defines the optimum usable frequency of operation. The maximum electron density of free electrons within the F2 layer and therefore foF2 depend upon the strength of the solar ionising radiation which is a function of time of day, season, geographical location and solar activity [START_REF] Goodman | HF Communications, Science and Technology[END_REF][START_REF] Maslin | The HF Communications, a Systems Approach[END_REF][START_REF] Mcnamara | Grid The Ionosphere: Communications, Surveillance, and Direction Finding[END_REF]. The E-region lies below the F-region extending approximately from 90 km to 150 km (figure 1). Typically the E-region is characterised by lower electron densities than the F-region (continuous profile in figure 1). However sometimes an irregular layer of extremely high density layer called sporadic E (Es) layer appears (dashed part in the profile in figure 1) in the E-region which prevents any radio signals from reaching the F-region. Es can occur during daytime or night-time, and its characteristics vary with latitude. These regions can apparently be caused by several mechanisms, and have a wide variety of characteristics. They have been correlated in literature with solar activity, geomagnetic activity and associated with thunderstorms and meteor showers [START_REF] Mathews | Sporadic E: current views and recent progress[END_REF][START_REF] Whitehead | Recent work on mid-latitude and equatorial sporadic-E[END_REF].

This paper describes the development of a neural network tool to interpolate missing foF2 data over Cyprus. The tool development is based on 10342 hourly foF2 measurements recorded above Cyprus from September 2008 to August 2010. The practical application of this model lies in the fact that in the presence of Es above Cyprus, due to the inability of the ionosonde to perform foF2 measurements, this model can provide an efficient interpolation method to fill the gaps in the foF2 data series that is used to drive long-term prediction models of foF2. There have not been any other efforts to apply interpolation techniques on foF2 datasets apart from a specific dataset over UK [START_REF] Francis | Nonlinear Prediction of the Hourly foF2 Times Series in Coonjunction with the Interpolation of Missing Data Points[END_REF]. This can be partly explained by the fact that most ionospheric stations are situated at middle latitudes where Es is not a frequent phenomenon and so data gaps may only appear because of system failure. In the case of Cyprus which is the lowest European ionospheric station in operation, the Es phenomenon is very frequent causing data gaps under normal system operation. A possible candidate for the interpolation procedure would have been the recently developed neural network approach from modelling foF2 [START_REF] Haralambous | A Neural Network model for the critical frequency of the F2 ionospheric layer over Cyprus[END_REF]. However the key novelty of the proposed interpolation tool is the use of TEC measurements as an additional input resulting in improved performance compared to the existing longterm (without TEC as input) model in terms of interpolation. Measurements of foF2 are conducted by ionosondes. which sweep through the HF frequency range transmitting short pulses. These pulses are reflected at various layers of the ionosphere, and their echoes are received by the receiver giving rise to a corresponding plot of reflection altitude against frequency (called an ionogram-figure 2) which is further analysed to infer the ionospheric plasma height-electron density profile (figure 1). The maximum frequency at which an echo is received is called the critical frequency of the corresponding layer. Since the F2 layer is the most highly ionised ionosperic layer its critical frequency foF2 is the highest frequency that can be reflected by the ionosphere (figure 2b). The presence of a Es layer (figure 2a) which does not allow for ionosonde signals to reach F2 region altitudes does not allow a useful ionogram to be obtained and therefore gaps in the data series of foF2 are formed (figure 3b). These gaps have to be interpolated in a way to preserve the inherent variability of foF2 data. This study suggests that one way to achieve this is to incorporate Total Electron Content (TEC) measurements which is another type of ionospheric characteristic the measurement of which is not affected by Es layer. The reason is that the signals involved in the measurement of TEC are high enough (> 1 GHz) to pass through the Es layer practically unaffected as opposed to the much lower frequency of ionosonde signals (<20 MHz). Dual-frequency GPS data recorded by GPS (Global Positioning System) receivers enable an estimation of TEC because of the frequency dependent delay imposed on the signal due to the ionosphere [START_REF] Wild | Estimating the Ionosphere Using One or More Dual Frequency GPS Receivers[END_REF]. By processing code and phase measurements on two frequencies in the L-band ( L1=1575.42 MHz, L2=1227.60 MHz) it is possible to extract an estimate of the TEC measured in total electron content units, (1 TECU = 10 16 el m-2). This is the total amount of electrons along a particular line of sight between the receiver and a GPS satellite in a column of 1 m 2 cross-sectional area. TEC is therefore the integral of the electron density profile (figure 1) from the Eregion up to the height of the GPS satellite encompassing the electrons in the F-region as well.

The analysis used in the present work to estimate TEC from GPS transmissions was carried out by means of the procedure developed by Ciraolo [START_REF] Ciraolo | Evaluation of GPS L2-L1 biases and related daily TEC profiles[END_REF][START_REF] Ciraolo | Results and problems in GPS TEC evaluation[END_REF]. It can be shown that the measurable GPS ionospheric delay is expressed in terms of TEC as [START_REF] Davis | Ionospheric Radio[END_REF]:

TEC cf t ion 2 3 . 40  (1) and   2 1
). (

h h dh h N TEC (2)
t ion … GPS ionospheric delay c… velocity of light in free space f… frequency of radio wave propagating through the ionosphere (GPS radio signal) N(h) … ionospheric electron density at the height h above the earth's surface TEC … total electron content along the satellite ray path.

Characteristics of the F2 layer critical frequency and Total Electron Content

Solar activity has an impact on ionospheric dynamics which in turn influence the electron density of the ionosphere. The electron density of the F2 layer and TEC exhibit variability on daily, seasonal and long-term time scales in response to the effect of solar radiation. It is also subject to abrupt variations due to enhancements of geomagnetic activity following extreme manifestations of solar activity disturbing the ionosphere from minutes to days on a local or global scale.

The most profound solar effect on foF2 and TEC is reflected on their daily variation as shown in figure 4. As it is clearly depicted, there is a strong dependency on local time which follows a sharp increase of foF2 and TEC around sunrise and gradual decrease around sunset. This is attributed to the rapid increase in the production of electrons due to the photo-ionization process during the day and a more gradual decrease due to the recombination of ions and electrons during the night. The long-term effect of solar activity on both parameters follows an eleven-year cycle and is clearly shown in figure 5(a),(b) where all values of foF2 and TEC are plotted against time as well as a modeled monthly mean sunspot number R which is a well established index of solar activity. We can observe a marked correlation of the mean level of both parameters and modeled sunspot number. In addition to the effects of solar activity on both parameters mentioned above we can also identify a strong effect on the diurnal variability as solar activity gradually increases through its 11-year cycle. This is demonstrated again in figure 4 where the diurnal variation of foF2 and TEC is plotted for three different days corresponding to low, medium and high sunspot number periods. It is evident from this figure that the night to day variability is increased as sunspot number increases for both parameters. There is also a seasonal component in the variability of foF2 and TEC which can be attributed to the seasonal change in extreme ultraviolet (EUV) radiation from the Sun. This can be clearly identified in figure 6 for noon values of both parameters for two different years. From the evidence drawn from this section the incorporation of TEC as an input to the neural network interpolation tool for foF2 data is justified based on its established correlation with TEC. However as shown in other studies this correlation varies diurnally and seasonally and also it is established that foF2 is not linearly related to TEC [START_REF] Spalla | TEC and foF2 comparison[END_REF].

Model parameters

The diurnal variation of foF2 is clearly evident by observing figure 4. We therefore include hour number as an input to the model. The hour number, hour, is an integer in the range 0 ≤ hour ≤ 23. In order to avoid unrealistic discontinuity at the midnight boundary, hour is converted into its quadrature components according to:

sin 2 24 hour sinhour         (3)
and

cos 2 24
hour coshour      (4)
A seasonal variation is also an underlying characteristic of foF2 as shown in figure 6 and is described by day number daynum in the range 1 ≤ daynum ≤ 365. Again to avoid unrealistic discontinuity between December 31st and January 1st daynum is converted into its quadrature components according to:

sin 2 365

daynum sinday         (5)
and cos 2 365

daynum cosday         (6)
Long-term solar activity has a prominent effect on foF2. To include this effect in the model specification we need to incorporate an index, which represents a good indicator of solar activity. In ionospheric work the 12-month smoothed sunspot number is usually used, yet this has the disadvantage that the most recent value available corresponds to foF2 measurements made six months ago. To enable foF2 data to be modelled as soon as they are measured, and for future predictions of foF2 to be made, the monthly mean sunspot number values were modeled using a smooth curve defined by a summation of sinusoids (figure 5). TEC is the last parameter used as an input to the neural network tool. It corresponds to the value of TEC measured simultaneously as the foF2 value. So the tool makes use of the available TEC values to generate the missing values of foF2 when the Es layer is present.

Experiments and Results

To test the performance of this approach in interpolating missing foF2 values we followed a 10-fold cross-validation process; the 10342 foF2 values (recorded above Cyprus from September 2008 to August 2010) and their corresponding inputs were split into 10 continuous parts of almost equal size and the predictions for each part were obtained by training a neural network on the other 9 parts. We chose to split the data into blocks of continuous values in order to test the performance of the approach in a worst-case scenario where the missing values cover a relatively big period of time.

The neural networks used were fully connected two-layer networks, with 6 input, 5 hidden and 1 output neuron. Both its hidden and output neurons had hyperbolic tangent sigmoid activation functions. The number of hidden neurons was determined by trial and error. The training algorithm used was the Levenberg-Marquardt backpropagation algorithm with early stopping based on a validation set created from 20% of the training examples. In an effort to avoid local minima ten neural networks were trained on each training set with different random initializations and the one that performed best on the corresponding validation set was selected for application to the test examples. The inputs and target outputs of the network were normalized setting their minimum value to -1 and their maximum value to 1. This made the impact of all inputs in the model equal and transformed the target outputs to the output range of the neural network activation functions. The results reported here were obtained by mapping the outputs of the network for the test examples back to their original scale. The whole process was implemented in Matlab with the aid of the Matlab Neural Network toolbox.

The obtained Root Mean Squared Error (RMSE) over all data was 0.508 MHz and the Mean Absolute Error (MAE) was 0.383 MHz. In order to see the improvement resulting from the inclusion of TEC in the inputs of the neural network (i.e. by using the long-term prediction approach of [START_REF] Haralambous | A Neural Network model for the critical frequency of the F2 ionospheric layer over Cyprus[END_REF]) we performed exactly the same experiment without including TEC. In this case the resulting RMSE was 0.712 MHz and the resulting MAE was 0.545 MHz. This shows that by using TEC we achieve a much more accurate interpolation of missing foF2 values. Finally examples of the approach are demonstrated graphically in figure 7.

Conclusions and future work

In this paper we have presented the development of a neural network tool for the interpolation of missing critical frequency data of the F2 ionospheric layer (foF2) above Cyprus due to sporadic E layer conditions. The model has been developed based on a data set obtained during a period of approximately two years. The tool has produced a good interpolation capability of the gaps in the foF2 data series therefore providing a method to preserve the variability of foF2, a fact which is essential in the development of long-term time-series prediction models and procedures. The tool will be further evaluated as we collect more foF2 measurements in the next few years.

Fig. 1 .

 1 Fig. 1. Typical electron density altitude profile of the ionosphere superimposed with enhanced density in the E-region due to Es layer.

Fig. 2 .

 2 Fig. 2. Ionogram (a) with Es blocking F2 layer reflections (b) with normal F2 reflections in the absence of Es layer.

Fig. 3 .

 3 Fig. 3. Data series of foF2 (a) without gaps (b) with gaps due to ionograms with Es layer.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Diurnal variability of (a) foF2 and (b)TEC for low, medium and high solar activity.

Fig. 6 .

 6 Fig. 6. Seasonal variability of (a) foF2 and (b)TEC.

Fig. 7 .

 7 Fig. 7. Examples of measured and interpolated foF2 values.


Leon Bobrowski
email: leon@ibib.waw.pl

Induction of Linear Separability Through the Ranked Layers of Binary Classifiers

Keywords: Learning sets, linear separability, formal neurons, binary classifiers, ranked 1

The concept of linear separability is used in the theory of neural networks and pattern recognition methods. This term can be related to examination of learning sets (classes) separation by hyperplanes in a given feature space. The family of K disjoined learning sets can be transformed into K linearly separable sets by the ranked layer of binary classifiers. Problems of the ranked layers deigning are analyzed in the paper.

Introduction

The Perceptron model and the error-correction learning algorithm of formal neurons played a fundamental role in the early neural networks [START_REF] Rosenblatt | Principles of neurodynamics[END_REF], [START_REF] Minsky | Perceptrons[END_REF]. The key concept here was the linear separability of learning sets. The convergence in a finite number of steps of the classical error correction algorithm depends on the linear separability of learning sets. Introducing a positive margin into error correction procedure allowed to ensure the convergence in a finite number of steps of the modified learning algorithms, also when the learning sets are not linearly separable [START_REF] Duda | Pattern classification[END_REF], [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF]. The perceptron criterion function linked to the error error-correction learning algorithm belongs to the family of convex and piecewise linear (CPL) criterion functions. Minimizing the CPL functions allows for, inter alia, effective designing of linear classifiers, linear forecasting models designing or carrying out selection of features subsets by using the relaxed linear separability (RLS) method [START_REF] Bobrowski | Feature selection based on relaxed linear separability[END_REF].

The most popular algorithms currently used in data mining are the support vector machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF]. The SVM algorithms are also linked to linear separability of learning sets. An essential part of these algorithms is the linear separability induction through the application of kernel functions. The selection of the appropriate kernel functions is still an open and difficult problem in many practical problems.

A family of K disjoined learning sets can be transformed into K linearly separable sets as a result of the transformation by ranked layer of formal neurons [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF], [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF]. This result can be extended on the ranked layer of arbitrary binary classifiers. The proof of the theorem concerning induction of linear separability through ranked layer of binary classifiers is given in this paper.

Separable Learning Sets

Let us assume that m objects O j (j = 1,….,m) are represented as the so called feature vectors x j = [x j1 ,...,x jn ] T , or as points in the n-dimensional feature space

F[n] (x j  F[n]
). Components (features) x i of the feature vector x represent numerical results of different measurements on a given object O (x i {0,1} or x i  R).

We assume that the feature vector x j (k) (j = 1,......, m) has been labelled in accordance with the object

O j (k) category (class)  k (k = 1,….,K). The learning set C k contains m k feature vectors x j (k) assigned to the k-th category  k C k = {x j (k)} (j  I k ) (1)
where I k is the set of indices j of the feature vectors x j (k) assigned to the class  k .

Definition 1. The learning sets C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF] are separable in the feature space

F[n], if they are disjoined in this space (C k  C k' = , if k  k').
This means that the feature vectors x j (k) and x j ' (k) belonging to different learning sets C k and C k' cannot be equal:

(k  k') (j  I k ) and (j I k ) x j (k)  x j' (k) (2)
We are also considering separation of the sets C k (1) by the hyperplanes H(w k , k ) in the feature space F[n]:

H(w k , k ) = {x: w k T x =  k }. ( 3
)
where w k = [w k1 ,....,w kn ] T  R n is the weight vector,  k  R 1 is the threshold, and (w k ) T x is the inner product.

Definition 2. The feature vector x j is situated on the positive side of the hyperplane

H(w k , k ) (3) if and only if (w k ) T x j   k .
(k  {1,...,K}) ( w k , k ) (x j (k)  C k ) w k T x j (k)   k . and (x j (i) C i , i  k) w k T x j (i) <  k (4)
In accordance with the inequalities (4), all the vectors x j (k) from the set C k are situated on the positive side of the hyperplane H(w k , k ) (3) and all the vectors x j (i) from the remaining sets C i are situated on the negative side of this hyperplane.

Layers of Binary Classifiers

The binary classifiers BC i (v) operate on feature vectors x (x F[n]) and are characterized by such a decision rule r i (v; x) which depends on the vector of parameters v (vR N ) and has the binary output of

r i = r i (v; x) (r i  {0, 1}, i = 1,…,m).
Definition 4. The activation field A i (v) of the i-th binary classifier BC i (v) with the decision rule r i (v;x) is the set of such a feature vectors x which activate this classifier

A i (v) = {x: r i (v; x) = 1} (5)
A layer composed of L binary classifiers BC i (v i ) with the decision rules r i (v i ;x) produces output vectors r( T . The layer of L binary classifiers BC i (v i ) transforms feature vectors x j (k) from the learning sets C k (1) into the sets R k of the binary output vectors r j (k) = r(V;x j (k)).

V;x) = [r 1 (v 1 ;x),…,r L (v L ;x)] T , where V = [v 1 T ,…,v L T ]
R k = {r j (k)} (j  I k ) (6)
Under certain conditions, it is possible to achieved the linear separability (4) of the sets R k [START_REF] Vapnik | Statistical Learning Theory[END_REF]. These conditions are analyzed in the paper.

Consider a few examples of binary classifiers BC i (v i ).

i. Formal neurons FN(w i , )

The formal neuron FN(w, ) with the weight vector w = [w 1 ,…,w n ] T  R n and the threshold  R) can be characterized by the below decision rule r FN (w, ; x):

if w T x  then r FN (w, ; x) = 1, else r FN (w, ; x) = 0 (7)
The formal neuron FN(w,) is activated (r FN (w, ; x) = 1) if and only if the weighed sum w 1 x 1 +...+ w n x n of n inputs x i x i  R) is greater than the threshold .

This decision rule r FN (w, ; x) depends on n + 1 parameters w i (i = 1,…,n) and . The activation field

A i (v) (5) is in this case the half space A FN (w,) = {x: w T x   ii. Logical elements LE i (w, ) The i-th logical element LE i (w,) (i = 1,…, L) can be defined as the formal neuron (7) with only one entry x k(i) x k(i)  R),
where k(i) is the number of the distinguished entry k(i){1,...,n}). The decision rule r LEi (w k(i) , i ; x) of the i-th logical element LE i (w,) can be given similarly to [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF] as:

if w k(i) x k(i)  then r LEi (w k(i) , i ; x) = 1, else r LEi (w k(i) , i ; x) = 0 (8)
The above decision rule r LEi (w k(i) , i ; x) depends on two parameters: w k(i) and

 i . The hyperplane H(w k(i) , i ) = {x: w k(i) x k(i) =  i } (3) in the feature space F[n]
is parallel to all axes x l except the axis x k(i) (l = 1,...,n and l ≠ k(i)) and perpendicular to the axis x k(i) .

iii. Radial binary classifiers RC(w 0i , i )

The decision rule r RC (w 0 ,; x) of the radial classifier RC(w 0 ,) can be defined in the below manner by using the distance (w 0 ,x) between the point x and the center of the ball

w 0 = [w 01 ,…,w 0n ]. if (w 0 ,x) ≤  then r RC (w 0 ,; x) = 1, else r RC (w 0 ,; x) = 0 (9)
The radial classifier RC(w 0 ,) is excited (r RC (w 0 ,; x) = 1) if and only if the distance (w 0 ,x) is no greater than the radius . The decision rule r RC (w 0 ,; x) (9) of the radial classifier RC(w 0 ,) depends on the n + 1 parameters w 0 = [w 01 ,…,w 0n ] T and . The activation field A RC (w 0 ,) (5) of the radial classifier RC(w 0 ,) is the ball with the center w 0 and the radius . The shape of this activation field A RC (w 0 ,) depends of the choice of the distance function (w 0 ,x). A few examples of distance functions (w 0 ,x) are given below [START_REF] Duda | Pattern classification[END_REF]:

 E (w 0 ,x) = ((x -w 0 ) T (x -w 0 )) 1/2 -the Euclidean distance  M (w 0 ,x) = ((x -w 0 ) T  -1 (x -w 0 )) 1/2 -the Mahalanobis distance  L1 (w 0 ,x) =  |w 0i -x i | -the L 1 distance i = 1,…,N (10)

Properties of the Ranked Layers of Binary Classifiers

The strategy of the ranked layers designing allows to find a layer with such a property that the sets R k (6) of transformed vectors r j (k) = r(V;x j (k)) become linearly separable (4). The proposed multistage designing procedure involves finding a sequence of admissible cuts by binary classifiers BC i (v i ) with the decision rules r i (v i ; x) [START_REF] Bobrowski | Feature selection based on relaxed linear separability[END_REF].

Definition 5. The binary classifier BC i (v i ) with the decision rule r i (v i ; x) (5) is admissible in respect to a given learning set C k from the family (1) if and only if n k (n k  0) elements x j (k) of this set activate the classifier and none of elements

x j (k) from other sets C k (k ≠ k) activate the classifier BC i (v i ) (k {1,…,K}) (x j (k) C k ) r i (v i ; x j (k)) = 1, and
(x j (k)  C k , where k ≠ k) r i (v i ; x j (k)) = 0 (11)
The activation field A i (v i ) (4.1) of the admissible binary classifier BC i (v i ) contains some elements x j (k) of only one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF].

Let us describe the multistage procedure of the ranked layer designing from binary classifiers BC i (v) on the basis of the learning sets C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF].

Ranked layer designing procedure.

v i(n) of classifier BC i(n) (v i(n) ) admissible (11) to some set C k(n) [n]. Stage iii. (Admissible reduction (cut) of the set C k(n) [n]
) -delay such feature vectors x j which activate the classifier It can be seen that each binary classifiers BC i (v i ) added to the layer reduces (13) the set C k(n) [n] by at least one vector x j (k(n)). Based on this property, it can be proved that if the learning sets C k (1) are separable (2), then after finite number L of steps the procedure will be stopped. The following Lemma results.

BC i(n) (v i(n) ) C k(n) [n+1] = C k(n) [n] -{x j (k(n)) C k(n) [n]: r i(n) (v i(n) ; x j (k(n)) = 1} (
Lemma 2. The number L of binary classifiers BC i (v i ) in the ranked layer is no less than the number K of the learning sets C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF] and no greater than the number m of the feature vectors x j (k) in these sets.

K  L  m (14)
The minimal number L = K of binary classifiers BC i (v i ) appears in the ranked layer when whole learning sets C k (1) are reduced (13) during successive steps n. The maximal number L = m of elements appears in the ranked layer when only single elements x j (k) are reduced during successive steps n.

In accordance with the postulate of a large admissible reduction, the number n k(n) [n] of the reduced elements x j (k(n)) (x j (k(n)C k(n) [n]) during each step n (13) should be as large as possible [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF]. In order to fulfill this postulate the below rule can be introduced:

The set C k(n) [n] reduced during the n-th designing step should be characterized by the largest number n k(n) [n] of reduced elements x j (k(n)):

(k {1,…,K}) n k(n) [n]  n k [n] (15)
The sequence of the parameters v i(1) appears as the result of the procedure (12):

v i(1) , v i(2) ,..., v i(L) (16)
Additionally, the class indices k(n) of the sets C k(n) [n] reduced (13) during successive steps n can be obtained

k(1), k(2),..., k(L) (17)
The parameters v i(n) (16) define the decision rule

r i(n) = r i(n) (v i(n) ;x) (r i(n)  {0, 1}) of the binary classifiers BC i(n) (v i(n) ) in the ranked layer BC i(1) (v k(1) ), BC i(2) (v i(2) ),..., BC i(L) (v i(L) ) (18)
Feature vectors x j (k) from the learning sets C k (1) are transformed by the ranked layer of L classifiers (18) into vectors r j (k) = [r j1 ,…,r jL ] T with L binary components r ji (r ji  {0, 1}) which are related to the category  k .

(k {1,…,K}) (x j (k)  C k ) r j (k) = [r i(1) (v i(1) ;x j (k)),…, r i(L) (v i(L) ;x j (k))] T (19)
The transformed vectors r j (k) (19) can be represented in the below manner by Ldimensional vectors q n = [q n1 ,…, q nL ] T (l = 1,…, L) with the binary components q lj

q 1 = [1, q 12 ,…,r lL ] T q 2 = [0,1, q 23 ,…,r 2L ] T ……. ……. q L = [0,0,…,0,1] T (20)
It is assumed in the above representation that each component q nj with j  n can be set to any binary value q nj = 1 or q nj = 0. In result, each vector q n (20) can represent more than one transformed vector r j (k) (19).

The vectors q n = [q n,1 ,…,q n,L ] T (20) can be used in the below decision rule, which results from the procedure (12) of the ranked layers designing if (i: 0  i  n) q n,i = 0 and q n,n = 1, then q n   k(n) (21) where k(n) is the index (17) of the set C k(n) reduced during the n-th stage.

The rule (21) links the vectors q n (20) to particular classes  k . This rule is based only on the content of the n-th admissible cut (13) of the set C k(n) . In this manner, each class  k can be represented by the set Q k of the related vectors q n (k) (20).

Q k = {q n (k): k(n) = k}, ( 22
)
where k  {1,…, K} and the distinguished index k(n) is determined by ( 17).

Theorem 1. The sets Q k (22) are linearly separable (4).

Proof: The proof is based on the example of such parameters w k ,  k , which fulfil the inequalities (4) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF]. Let us introduce the L-dimensional vector a = [a 1 ,…,a L ] T with the components a i specified below

(i  {1,…, L}) a i = 1/2 i (23)
The weight vector w k = [w k1 ,…,w kL ] T is defined on the basis of the sequence (17) in accordance with the below rule

(i {1,…, L}) if k(i) = k, then w ki = a i else w ki = -a i ( 24
)
It can be directly verified that all the inequalities (4) are fulfilled by the vectors w k with the components w ki (24) and the threshold  k = 0. This means that the sets Q k (22) are linearly separable (4).

Corollary. The sets R k = {r j (k)} ( 6) of the transformed vectors r j (k) ) related to particular categories  k are linearly separable [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF].

(k  {1,…, K}) (w k ) (r j (k) R k ) (w k ) T r j (k)  0 and (r j (k) R k , k ≠ k) (w k ) T r j (k) 0 (25)
The transformation (25) of the feature vectors x j (k) by the ranked layer allows to replace the separable (2) learning sets C k (1) by the sets R k = {r j (k)} ( 6) which are linearly separable (25). Let us remark that the Theorem 1 is also valid if the ranked layer is composed of different types of binary classifiers BC i (v i ).

Examples of Implementations of the Ranked Layer Designing Procedure

Let us consider at the beginning, the procedure ranked layer designing from logical elements LE i (w,) (8). The decision rule r LEi (w k(i) , i ;x) (8) of the i-th logical element LE i (w,) depends on only one feature x k(i) . The decision rule r LEi (w k(i) , i ;x) (8) can be decomposed into two rules r i + (x) and r i -(x) and represented in the below manner for each feature x i (i = 1,…,n):

(i  {1,…,n}) if x i   i + , then r i + (x) = 1 else r i + (x) = 0 (26) and (i  {1,…,n}) if x i   i -, then r i -(x) = 1 else r i -(x) = 0 (27)
where  i + and  i -are two numerical parameters ( i

+  R 1 and  i - R 1 ).
In accordance with the admissible reduction principle (11), the parameters  i + and  i -should ensure the homogenous reduction of the maximal numbers n + (i) and n -(i) of elements x j (k) belonging to only one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF]. The term homogenous reduction (11) means that the condition r i + (x) = 1 (26) is fulfilled by n + (i) elements x j (k) of only one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF]. Similarly, the condition r i -(x) = 1 (27) should be fulfilled by n -(i) elements x j (k) of only one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF].

In order to satisfy the postulate of a large admissible reduction the below rules can be used

(i  {1,…,n}) n(i) = max{ n + (i), n -(i)}, and n(i) = max{ n(1),…, n(L)}, ( 28
)
where n(i) is the maximal number of elements x j (k) from one learning set C k (1), which can be reduced by using the i-th feature x i (coordinates x i )and the rule (26) or (27), i is the index of such a coordinate x i which allows to reduce the maximal number n(i) of elements x j (k) belonging to only one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF].

The rule (28) allows to determine the class indices k(n) (17) of the sets C k(n) [n] reduced (13) during successive steps n of the Ranked layer designing procedure (12). The optimal threshold value (n) can be determined on the basis of the rules (26) and ( 27). The optimal threshold value (n) should result in the maximal number n(i) (28) of elements x j (k) of one learning set C k (1) reduced during each step n. The rules (28) and include the selection of the optimal threshold values  i + and  i -for each coordinate x i . The optimal threshold values  i + (26) and  i -(27) should result in the maximal numbers n + (i) and n -(i) of the reduced elements x j (k) of one learning set C k [START_REF] Rosenblatt | Principles of neurodynamics[END_REF]. The implementation of the ranked layers design rules is relatively simple in the case of logical elements LE i (w,) (8). For small data sets [START_REF] Rosenblatt | Principles of neurodynamics[END_REF], such rules can even be implemented without a computer, based on handwritten calculations.

Let us consider now the layer of formal neurons FN(w i ,  i ) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF] which can be linked to the hyperplanes H(w i , i ) (3). Definition 6. The hyperplane H(w k , k ) (3) is admissible (Definition 5) to a given set C k [n] from the family {C k [n]} (13) if and only if n k (n k  0) elements x j (k) of this set is situated on the positive side of the hyperplane H(w k , k ) and all elements x j (k) from other sets C k (k ≠ k) are situated on the negative side of this hyperplane:

(k {1,…,K}) (x j (k) C k [n]) w k T x j (k)   k , and
(x j (k)  C k [n], where k ≠ k) w k T x j (k) <  k (29)
The sequence of admissible hyperplanes H(w i , i ) (3) should be consisted with the Ranked layer designing procedure (12). The parameters w i and  i of admissible hyperplanes H(w i , i ) (3) define the formal neurons FN(w i , i ) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF] of the ranked layer.

The hyperplane H(w k , k ) (3) admissible (Definition 5) to a given set C k [n] from the family {C k [n]} (13) can be obtained through minimization of the convex and piecewise linear (CPL) criterion functions  k" (w,) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF]. The preceptron criterion function belongs to the CPL family [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF]. The criterion function  k (w,) can be defined on the basis of the positive G k + and the negative G k -sets of feature vectors x j from the sets {C k [n]} (13) G k + = {x j :

x j  C k [n]} and G k -= {x j : x j   C k [n]} k ≠ k (30)
Each element x j of the set G k + defines the positive penalty function  j + (w,) The basis exchange algorithms which are similar to the linear programming allow to find efficiently minimum  k' (w k  , k  ) of the criterion function  k (w,) (33) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF]. It can be shown that the optimal parameters (w k * , k * ) (34) obtained from the function  k (w,) (39) with a sufficiently large parameter  define the admissible hyperplane H(w k  , k  ) (29) [START_REF] Bobrowski | Eksploracja danych oparta na wypukłych i odcinkowo-liniowych funkcjach kryterialnych (Data mining based on convex and piecewise linear (CPL) criterion functions)[END_REF]

(x j  G k + ) if w T x j -  1, then  j + (w,) = 1 -w T x j + , else  j + (w,) = 0 (

Concluding Remarks

Separable learning sets C k (2) can always be transformed into such sets R k (6) which are linearly separable (4). The linear separability is induced through the transformation of elements x j (k) of the learning sets C k (1) by the ranked layer of binary classifiers BC i (v i ). The paper contains a proof of this property.

The paper also contains examples of designing ranked layers of formal neurons FN(w i ,  i ) [START_REF] Bobrowski | Design of piecewise linear classifiers from formal neurons by some basis exchange technique[END_REF] and logical elements LE i (w i , i ) (8). Efficient designing ranked layers of radial classifiers RC(w 0 ,) (9) is still an open problem from computational point of view.

  Stage i. (Start) -put n =1 and define the sets C k [n]: (k  {1,...,K}) C k [n] = C k Stage ii. (Admissible classifier) -find parameters

13 )

 13 Stage iv. (Stop criterion) if all the sets C k [n+1] are empty then stop else increase the index n by one (n  n+1) and go to the Stage ii. 

  31)Similarly, each element x j of the set G k -defines the negative penalty function  j -(w,)(x j  G k -) if w T x j -  -1, then  j -(w,) = w T x j -  -1, else  j -(w,) = 0 (32)The criterion functions  k  (w,) is the sum of the functions  j + (w,) and  j -(w,) k (w,) =   j + (w,) +    j -( w,)(33)where  ( > 0) is a positive parameter (price).Minimization of the function  k (w,) allows to find optimal parameters ( k (w,) =  k (w k  , k  )  0 (34)

  . The optimal parameters (w k * , k * ) allows to define the formal neurons FN(w k * , k * ) (7) of the ranked layer in accordance with the Ranked layer designing procedure (12).

  Similarly, vector x j is situated on the negative side of H(w k , k ) if and only if (w k ) T x j <  k .

	Definition 3. The learning sets (1) are linearly separable in the n-dimensional feature
	space F[n] if each of the sets C k can be fully separated from the sum of the remaining
	sets C i by some hyperplane H(w k , k ) (3):

Acknowledgment. This work was supported by the by the NCBiR project N R13 0014 04, and partially financed by the project S/WI/2/2011 from the Białystok University of Technology, and by the project 16/St/2011 from the Institute of Biocybernetics and Biomedical Engineering PAS.