H Haddad

W Leclerc

G Alhajj

A Ammar

C Pélegris

M Guessasma

E Bellenger

Numerical investigation of heat conduction in heterogeneous media with a discrete element method approach

Keywords: Discrete element method, heat conduction, effective thermal conductivity, thermal heat flux

Composite materials have been widely used across industry sectors. However, they are characterized by a variability of thermal conductivity with the architecture and manufacturing processes. Hence, thermal transfer in composite materials requires an improved fundamental understanding. From numerical purposes, the Finite Element Method (FEM) seems a robust method to study heat transfer in composite material. However, it does not establish a high-fidelity with the real life of material since it is difficult to take into account potential damage matrix or interface imperfection by this method. As such, this paper discusses the development of numerical approach based on the Discrete Element Method (DEM) to study heat transfer in composite materials. For that purpose, we consider a hybrid lattice model based on the equivalence between a particulate domain and a continuous medium. Several works have used the DEM to study heat transfer in a homogeneous and continuous medium. Through this contribution, we aim to extend this approach to investigate composite materials. The model is then validated in terms of temperature by comparison with numerical and experimental results through several applications. Furthermore, a special care taken in the evaluation of the heat flux density fields. To the knowledge of the authors, previous works did not interest to the examination of heat flux density when using the DEM. Indeed, this sensitive to the packing configuration and consequently always heterogeneous even if there typically homogeneous. To overcome this problem, an original smoothing technique called Halo is proposed and discussed in this work. Results exhibit the relevance of the proposed approach to evaluate both temperature and heat flux density fields with a good degree of precision compared with th FEM, the Fast Fourier Transformer (FFT) based method and experimental data.

Introduction

The number of applications of composite materials in the industry has increased significantly in recent years, due to low manufacturing costs as well as the customization of their mechanical and thermal properties. For example, when silica particles are added into a polymer matrix, they play an important role in improving electrical, mechanical and thermal properties of resulting composites [START_REF] Buchdahl | Mechanical properties of polymers and composites -vols. i and ii, lawrence e. nielsen, marcel dekker, inc., new york[END_REF]. Due to their complex microstructures, the detailed study of heat transfer turns out to be an arduous task. Therefore, technical analyses potentially focus on the macroscopic behaviour of such materials, dictated by effective properties such as ETC. The determination of the latter according to microstructural characteristics, individual phase properties and other relevant physical parameters is consequently of scientific and practical importance. For example, Chauhan et al. [START_REF] Chauhan | Effect of geometry of filler particles on the effective thermal conductivity of two-phase systems[END_REF] studied the effect of geometry of filler particles on the effective thermal conductivity in the context of polymer composites. Mishra et al. [START_REF] Mishra | An experimental investigation on the effect of particle size on the thermal properties and void content of solid glass microsphere filled epoxy composites[END_REF] investigated on the effect of particle size on the thermal properties and void content of solid glass microsphere filled epoxy composites. Thus, considerable experimental works have been reported on the subject of thermal conductivity in composites. Among this vast literature, numerous studies have focused on the influence of filler volume fraction. For example, Tavman [START_REF] Tavman | Effective thermal conductivity of isotropic polymer composites[END_REF] experimentally evaluated ETC of reinforced polyethylene composites as function of tin powder filler. Wong et al. [START_REF] Wong | Thermal conductivity, elastic modulus, and coefficient of thermal expansion of polymer composites filled with ceramic particles for electronic packaging[END_REF] led experimental tests to determine ETC and effective elastic properties of epoxy resin with different kinds of inclusions as silica, alumina and aluminum fillers. Liang et al. [START_REF] Liang | A new heat transfer model of inorganic particulate-filled polymer composites[END_REF] measured experimentally ETC of aluminum powder-filled phenol-aldehyde composites and graphite powder-filled phenol-aldehyde composites with filler volume fraction from 0 to 50 %. Boudenne et al. [START_REF] Boudenne | Anomalous behavior of thermal conductivity and diffusivity in polymeric materials filled with metallic particles[END_REF], investigated combined effects of inclusion concentration and size in the context of a polypropylene matrix filled with aluminum and copper fillers.

Many theoretical models [START_REF] Pietrak | A review of models for effective thermal conductivity of composite materials[END_REF] have also been developed to predict the ETC of composites. Among them, the Maxwell scheme [START_REF] Maxwell | A Treatise on Electricity and Magnetism[END_REF] is considered to be one of the best homogenization technique due to its applicability to cases of multiphase composites. It predicts effective conductivities for dilute random sphere systems and a low volume fraction typically less than 25%. The method developed by Gandarilla-Pérez et al. [START_REF] Gandarilla-Pérez | Extension of maxwell homogenization scheme for piezoelectric composites containing spheroidal inhomogeneities[END_REF] enabled to extent the classical Maxwell's homogenization to different geometrical shapes of inclusions. Authors also refer today to the Rayleigh sphere model [START_REF] Sec | On the influence of obstacles arranged in rectangular order upon the properties of a medium[END_REF] which includes thermal interaction between the particles to enhance Maxwell's model. Besides, Bruggeman [START_REF] Bruggeman | Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. i. dielektrizitätskonstanten und leitfähigkeiten der mischkörper aus isotropen substanzen[END_REF] developed the Effective Medium Theory (EMT) which covers a wide range of materials, and Agari et al. [START_REF] Agari | Estimation on thermal conductivities of filled polymers[END_REF] developed a model based on a generalization of parallel and series conduction models. Wong et al. [START_REF] Wong | Thermal conductivity, elastic modulus, and coefficient of thermal expansion of polymer composites filled with ceramic particles for electronic packaging[END_REF] verified that both schemes are more accurate than Maxwell model for a high ratio of filler in the context of SiO 2 /epoxy composite. Nielsen [START_REF] Nielsen | The thermal and electrical conductivity of two-phase systems[END_REF] proposed a modified Halpin-Tsai equation based on the generalized Einstein coefficient and a constant related to the relative conductivity of the phases. EMT , Agari model and Nielsen model were compared in the case of a Si 3 N 4 /epoxy composite [START_REF] He | Thermal conductivity of ceramic particle filled polymer composites and theoretical predictions[END_REF]. Results exhibit that EMT and Nielsen model better predict ETC for a filler concentration lower than the range of 15-20%, 25-30% respectively but Agari model is more suited for a higher volume fraction. Khan et al. [START_REF] Khan | Effective thermal conductivity of two-phase composites containing highly conductive inclusions[END_REF] developed a theoretical model to predict ETC of two-phase particulate composites containing highly conductive inclusions. More recently, in [START_REF] Seppälä | Efficient method for predicting the effective thermal conductivity of various types of two-component heterogeneous materials[END_REF], is presented an efficient method for predicting the effective thermal conductivity of various types of two-component heterogeneous materials. Finally, it is acknowledged that the third order model discussed by Beasley and Torquato [START_REF] Beasley | Bounds on the conductivity of a suspension of random impenetrable spheres[END_REF] enables to compute ETC of 2-phase particulate composite with a higher degree of precision, whatever the volume fraction of filler and the relative conductivity. Nevertheless, such a scheme takes benefit of a threepoint parameter which requires very costly numerical computations to be accurately determined [START_REF] Gillman | Third order model of thermal conductivity for random polydisperse particulate materials using well-resolved statistical descriptions from tomography[END_REF].

Moreover, several numerical approaches treated heat transfer by conduction and most of them focused on the ETC of composites [START_REF] Thiele | Effective thermal conductivity of three-component composites containing spherical capsules[END_REF][START_REF] Zhang | Effectiveness of the heat conduction reinforcement of particle filled composites[END_REF]. Typically, researchers take benefit of the Finite Difference Method (FDM) [START_REF] Yinping | Numerical analysis of effective thermal conductivity of mixed solid materials[END_REF][START_REF] Spittle | Numerical prediction of the effective thermal conductivity of dendritic mushy zones[END_REF] and the FEM [START_REF] Jiang | Apparent thermal conductivity of periodic two-dimensional composites[END_REF][START_REF] Nayak | A computational and experimental investigation on thermal conductivity of particle reinforced epoxy composites[END_REF][START_REF] Agrawal | Computational, analytical and experimental investigation of heat conduction of particulate filled polymer composite[END_REF] to attend their objectives. Thus, the FDM presented in [START_REF] Yinping | Numerical analysis of effective thermal conductivity of mixed solid materials[END_REF] is suitable for calculating the ETC of mixed solid materials for different shapes of phases. It was also developed in [START_REF] Spittle | Numerical prediction of the effective thermal conductivity of dendritic mushy zones[END_REF] to predict the ETC of dendritic mushy zones as a function of fraction solid. The FEM was investigated in [START_REF] Saini | Numerical study using finite element method for heat conduction on heterogeneous materials with varying volume fraction, shape and size of fillers[END_REF] to study heat conduction in heterogeneous materials with varying volume fraction, shape and size of fillers. The finite approach was also used in [START_REF] Agrawal | Computational, analytical and experimental investigation of heat conduction of particulate filled polymer composite[END_REF] to calculate the ETC of polymer composites with filler content ranging from 2.35 to 26.8 %. In , a numerical study using the FEM is proposed to However, other approaches also exist to determine ETC. For example, a numerical method based on FFT was developed to evaluate the effective properties of composite with periodic microstructure composed of constituents with linear or nonlinear mechanical behaviour [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Leclerc | A numerical investigation of effective thermoelastic properties of interconnected alumina/al composites using fft and fe approaches[END_REF]. Besides, the DEM has also been applied is such a context [START_REF] Vargas | Heat conduction in granular materials[END_REF][START_REF] Haddad | Application of the discrete element method to study heat transfer by conduction in particulate composite materials[END_REF][START_REF] Terreros | Simulation of continuum heat conduction using DEM domains[END_REF][START_REF] Hahn | Discrete element representation for the thermal field: Proof of concept and determination of the material parameters[END_REF]. Indeed, this latter naturally accounts for discontinuities and is all well-suited to deal with a large spectrum of multi-scale and multiphysical phenomena. This approach originally developed by Cundall [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] to treat issues in rock mechanics field has been later applied to numerous problems, from the silo discharge [START_REF] Nicot | On the definition of the stress tensor in granular media[END_REF] to the tribology of the wheel-rail contact [START_REF] Descartes | A new mechanical-electrical approach to the wheel-rail contact[END_REF] through the defect diagnosis in ball bearings [START_REF] Machado | Electromechanical modeling by DEM for assessing internal ball bearing loading[END_REF][START_REF] Machado | Electromechanical prediction of the regime of lubrication in ball bearings using discrete element method[END_REF] and the simulation of crack propagation in continuous media [START_REF] Chen | A new thermo-mechanical coupled DEM model with non-spherical grains for thermally induced damage of rocks[END_REF][START_REF] Hentz | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF][START_REF] Radi | Elasticity and fracture of brick and mortar materials using discrete element simulations[END_REF]. Besides, a DEM based model was also investigated to simulate heat transfer by conduction in particulate systems using a contact conductance model introduced at the scale of the elementary inter-particular bond [START_REF] Vargas | Heat conduction in granular materials[END_REF][START_REF] Feng | Discrete thermal element modelling of heat conduction in particle systems: Basic formulations[END_REF][START_REF] Zhang | A DEM study on the effective thermal conductivity of granular assemblies[END_REF]. A similar approach based on the Finite Discrete Element Method (FDEM) is presented in [START_REF] Joulin | Capturing heat transfer for complex-shaped multibody contact problems, a new fdem approach[END_REF]. Other works used a similar paradigm to consider heat transfer by conduction in continuous media using the concept of representative element [START_REF] Haddad | Application of the discrete element method to study heat transfer by conduction in particulate composite materials[END_REF][START_REF] Terreros | Simulation of continuum heat conduction using DEM domains[END_REF][START_REF] Hahn | Discrete element representation for the thermal field: Proof of concept and determination of the material parameters[END_REF]. Thus, some authors [START_REF] Haddad | Application of the discrete element method to study heat transfer by conduction in particulate composite materials[END_REF][START_REF] Hahn | Discrete element representation for the thermal field: Proof of concept and determination of the material parameters[END_REF] explored the benefit of a Voronoï tessellation to associate a well-defined polyhedron to each particle while other researchers considered Platonic solids [START_REF] Terreros | Simulation of continuum heat conduction using DEM domains[END_REF] for the same purpose. Furthermore, DEM has also proved its ability to simulate the thermo-elastic behavior of continuous materials as well as the thermal induced damage due to thermal expansion mismatch in continuous heterogeneous media [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Leclerc | On the suitability of a discrete element method to simulate cracks initiation and propagation in heterogeneous media[END_REF][START_REF] Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF].

The objective of this contribution is to study the ability of the DEM to simulate heat transfer by conduction in a continuous and heterogeneous 3D media. For this purpose, we consider a model based on the equivalence between a particulate system and a continuous media. In such an approach, the duality between particulate and lattice descriptions enables to introduce heat transfer by conductionbetween each pair of particles in contact using the natural inter-particular link. Thus, section and thermal properties associated to each contact manage the heat transfer by conduction through the entire 3D domain. Besides, we use the concept of representative element in which each discrete element is represented by a solid polyhedron with a number of faces equal to its number of neighbors. Thus, it is possible to provide proper transmission surface to each contact and a suitable mass to each particle. Due to limitations exhibited by Platonic solids which do not ensure conservation of mass and the high computational cost of the generation of Voronoï tessellations, we propose a third approach based on the calibration of transmission surface under the hypothesis of the equality between local and global conductivities. Furthermore, we are also interested in determining the heat flux field in 3D heterogeneous media. In DEM, this latter is typically heterogeneous even if this is theoretically homogeneous. In order to tackle such a difficulty, we propose to handle and adapt the Halo approach [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D discrete element method simulation : Validation and application to flax/bio based epoxy composite[END_REF] recently developed to better evaluate the distribution of local stress in the context of linear elasticity. Please notice that, in this work, granular packings are generated using Lubachevsky-Stillinger Algorithm (LSA) [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] which enables an accurate control of intrinsic parameters such as the compacity, the size of particles and the coordination number which is defined as the average number of contacts per DE. All simulations are carried out using the house-built parallel MULTICOR3D++ code [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF].

The present paper is outlined as follows. In Section 2, the numerical modeling based on the DEM is described and validated for a homogeneous media by comparison with analytical solutions. The estimation of the heat transmission surface between two particles based on a numerical calibration is also presented. In Section 3, we introduce the Halo approach to better evaluate the heat flux field distribution in a homogeneous media. In Section 4, we perform validation tests in the context of a single-inclusion composite. Comparisons are performed with analytical and other numerical approaches in terms of ETC, temperature and heat flux fields using the Halo approach. In Section 5, we discuss the case of a heterogeneous media, namely epoxy resin filled with ceramic fillers. ETC of the heterogeneous media is evaluated using the discrete approach and compared with analytical, numerical and experimental results. Afterwards, the heat flux is investigated, and comparisons with FEM are set up.

Heat transfer in continuous domain

In this section, we present a numerical approach based on the DEM to model thermal conduction in a continuous and homogeneous 3D media. The latter consists of a polydisperse particulate packing which can be considered as an Equivalent Continuous Medium (ECM) under several assumptions related to the arrangement of particles. For this purpose, intrinsic parameters such as the coordination number, the compacity and the polydispersity are controlled by the LSA. Besides, the isotropy of the particulate system is preliminary verified according to the process described by Leclerc et al. [START_REF] Leclerc | A numerical investigation of effective thermoelastic properties of interconnected alumina/al composites using fft and fe approaches[END_REF][START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF] using statistical tools as the polar plots or the 2-point probability function to ensure that the effective properties are independent of direction. For information purposes, a volume fraction close to 0.64 and an average coordination number Z of 6.2 are typical of a 3D system of monodisperse particles. However, in the present work, we have deliberately chosen an average coordination number equal to 7.5 to ensure that each particle is at least in contact with 2 particles locally. This corresponds to the generation of new contacts and allows for a denser network of contacts. Moreover, the particle's radius follows a Gaussian distribution law and the dispersion is described by a ratio between the standard deviation and the average radius set to 0.3 in order to avoid possible directional effects.

DEM-based approach

In the proposed approach, we associate a representative element to each particle of the particulate system. Please notice that we are not interested in the shape of this element but rather in its characteristics such as volume and heat transmission surfaces. Thus, the volume of each representative element is chosen in order to satisfy the mass conservation hypothesis. In this purpose, the sum of the volumes of the representative elements must be equal to the volume of the studied domain. To do this, the mass of each particle is adjusted to the mass of the corresponding representative element by φ the volume fraction of the particulate system. If ρ c [kg/m 3 ] is the density of the continuous material, the density of the particulate packing ρ d [kg/m 3 ] can then be related to ρ c via the following relation:

ρ c = φρ d (1)
The heat flux W p,q [W] transmitted by the transmission surface S t p,q [m 2 ] between two representative elements (Figure 1) associated to particles p and q in contact, is defined as follows:

W p,q = H p,q c (T q -T p ) (2)
where T p , T q are the temperatures of particles p, q and H p,q c [W/K] is the thermal conductivity coefficient given by:

H p,q c = λS t p,q d p,q (3)
with λ the thermal conductivity of material and d p,q the distance between the centers of p and q. Thus, the heat transfer by conduction equation associated to each particle is:

C d p dT p dt + np q=1 W p,q = Q tot p ( 4
)
with

C d p = Cρ d V p /φ (5)
where Q tot p [W] represents the external heat flux associated to particle p, n p is the number of particles in contact with particle p, V p is the volume of particle p and C represents the specific heat of the constitutive material. The temporal resolution of Equation ( 4) can be achieved by two types of numerical schemes (implicit or explicit). Implicit numerical patterns are stable for large steps of time. However, because of the large number of particles in contact, implicit numerical schemes are too expensive. As a result, our choice was to use an explicit numerical scheme that does not require an iterative procedure. The time discretization of Equation ( 4) leads to: Calibration process Unlike the FEM, for which the local properties at the finite element scale are identical to the macroscopic properties for a homogeneous media, local properties and parameters in discrete approach must be correlated to macroscopic properties. In the present work, the estimation of the heat transmission surface between two particles in contact is obtained by calibration. This choice is less expensive compared to approaches that are based on the generation of representative elements from a Voronoï decomposition of the domain [START_REF] Haddad | Application of the discrete element method to study heat transfer by conduction in particulate composite materials[END_REF]. It also ensures the conservation of the mass contrary to the approach based on Platonic solids [START_REF] Terreros | Simulation of continuum heat conduction using DEM domains[END_REF]. In our case, the thermal conductivity of the equivalent medium depends on S t p,q the transmission surface between two representative elements. We chose to correlate S t p,q to a coefficient C t and a surface S p,q = πR 2 m deducted from R m the mean radius of the contact particles so that the expression of S t p,q is:

T t+∆t p = T t p + ∆t Cρ d V p [Q p + np q=1 S t p,q λ d p,q (T t q -T t p )] Q tot p (6)
S t p,q = C t S p,q (7)
Thus, the thermal conductivity of the material is directly related to C t coefficient without any dimensional effect. The latter is calibrated in order to obtain an equality between the thermal conductivity of the material and that considered at the contact scale. To do this, we proceed as follows. A heat flux ϕ is applied to the upper surface of the sample of length L, in the Y direction while a constant temperature T 0 is applied to the lower surface. The conditions under which the test is performed are described in Figure 2. For simplicity purposes, calculations are first performed using an arbitrary C t (C arb t ) set to 1. At the stationary state, a temperature difference ∆T is determined at the top surface and the ETC of the material λ e is deducted from the following equation:

λ e = ϕL ∆T (8)
Finally, C t coefficient is deducted from the following equality:

C t C arb t = λ λ e (9)
We have evaluated C t for a large range of density number which reaches 8000000 DE for cubic and cylindric configurations. Each data is obtained considering a number of 4 realizations. The variation of C t is shown in Figure 3. C t is found between 0.75 and 0.84. It depends on the number of particles and there is little effects of the geometry. The dispersion of results is characterized by the Coefficient of Variation (CoV) which is defined as the ratio of the standard deviation to the mean value. Since the CoV is less than 0.2% whatever the number of particles, we can conclude that we don't need calibration for a granular packing with a volume fraction close to 0.64 and an average coordination number Z of 7.5.

Time step calculation

Based on the Equation ( 6), the explicit time integration for an interior particle p can be expressed as:

T t+∆t p = T t p + α∆t V p np q=1 S t p,q d p,q (T t q -T t p ) (10)
where α = λ Cρ d designates the thermal diffusivity of the constitutive material. The explicit scheme is not unconditionally stable and the largest permissible value of the time step ∆t is limited by a stability criterion. The time integration step is then limited by a critical time step ∆t cr . We assume that S t p,q ≈ C t πR 2 m , V p ≈ 4 3 πR 3 m and d p,q ≈ 2R m with R m the mean radius of particles. Furthermore, for a local coordination number n p , Equation [START_REF] Maxwell | A Treatise on Electricity and Magnetism[END_REF] simplifies to :

T t+∆t p ≈ T t p + α∆t 4 3 πR 3 m np q=1 C t πR 2 m 2R m (T t q -T t p ) (11)
Moreover, Equation ( 11) can be rearranged as:

T t+∆t p ≈ T t p + 3C t α∆t 8R 2 m np q=1 (T t q -T t p ) (12)
We now define a dimensionless Fourier number τ as:

τ = 3C t α∆t 8R 2 m ( 13
)
It characterizes the part of the heat flux transmitted to a body at a given time t compared to the heat stored by the body. Then, we can reduce Equation ( 4) to:

T t+∆t p ≈ T t p + τ np q=1 (T t q -T t p ) = (1 -τ n p )T t p + τ np q=1 T t q ( 14
)
The stability criterion requires that the coefficient multiplied by T t p in the T t+∆t p expression must be greater or equal to zero [START_REF] Cengel | Heat and mass transfer: a practical approach, 3rd edn[END_REF]. Thus, (1 -τ n p ) ≥ 0 or τ ≤ 1 np . The integration time step is deduced :

∆t ≤ 8R 2 m 3C t αn p ( 15
)
The critical time step ∆t cr for the solution of a thermal problem can be estimated by : ∆t cr = 2(d min p,q ) 2 3αn p C t [START_REF] He | Thermal conductivity of ceramic particle filled polymer composites and theoretical predictions[END_REF] with d min p,q the minimal interparticle distance.

Validation of the DEM-based approach

In this section, we propose to validate the model in the context of a homogeneous media by comparison with an analytical solution. To do this, we consider a cubic pattern of length L = 10cm composed of 100, 000 spherical and polydisperse particles. This density of DE is a good compromise to ensure the isotropy of the particulate system and the accuracy of results while keeping a reasonable computational cost. C t coefficient for this configuration is set to 0.806 according to the previous study. The pattern is subject to the following thermal conditions:

   T 1 : T (y = 0) = 25 • C T 2 : T (y = L) = 35 • C t = 0 : T (y) = T 0 = 25 • C 0 < y < L
Besides, the lateral boundaries are under adiabatic conditions (Figure 4) and the material properties are given in Tab. 1:

Density ρ c 7800 kg/m 3
Thermal conductivity λ 33 W/(mK) Specific heat C 0.9 J/(kgK)

Tab. 1: Material properties The full development of the analytical resolution is detailed in [START_REF] Weigand | Analytical Methods for Heat Transfer and Fluid Flow Problems[END_REF]. The temperature field at any point is given by the following equation:

T (y) = T 1 + y∆T L + ∆T ∞ n=1 C n sin nπy L e -( nπ L )
2 λt ρC [START_REF] Khan | Effective thermal conductivity of two-phase composites containing highly conductive inclusions[END_REF] where

C n = 2 × (-1) n ×(1-θ0)+θ0
nπ and θ 0 is defined as a dimensionless quantity: θ 0 = T0-T1 T2-T1 .

The temperature variation obtained by the analytical solution and the DEM after 0.05s, 0.1s and 0.2s is represented graphically in Figure 5-a. The results provided by the DEM are in good agreement with the analytical solution with a maximum relative error equal to 0.7 % (Figure 5-b). This validates the proposed formulation in the context of a homogeneous media. A second validation test is performed in order to prove the efficiency of the proposed approach. We have considered a spherical pattern with diameter D = 10cm composed of aluminum. The initial temperature is T 0 = 20 • C. The sphere is subjected to a constant temperature of T i = 200 • C at the external surface. The specific heat is assumed to be 900 J/(kgK) . The thermal conductivity of the material is equal to 185 W/(mK) while the considered density is 2700 Kg/m 3 . The conditions under which the test is performed are described in the figure below. C t = 0.806 when considering granular packing with 100000 DE and C t = 0.834 for granular packing composed of 1 million DE. The analytical solution is given by the following equation:

T (r, t) = T 0 + 2(T 0 -T i ) πr ∞ n=1 (-1) n n sin( nπr R )e -a(πn) 2 t R 2 ( 18
)
where r is the distance of the studied point from the center of the sphere and R is the radius of the sphere. The temperature T c at the center of the sphere at time t is obtained when r approaches zero. It is expressed as:

T c (t) = T 0 + 2(T 0 -T i ) ∞ n=1 (-1) n e -a(πn) 2 t R 2 (19)
Comparisons are made between the DEM with 100000 and 1 million particles, and analytical solution and FEM calculation with 100000 tetrahedral finite elements. In Figure 7. We have plotted the evaluation of the temperature at the center of the sphere. Results indicate that DEM matches closely to both analytical solution and the FEM calculations.

Halo approach to evaluate the heat flux

In a stationary state, the expression of the heat flux ϕ i p applied to a particle p is defined as below:

ϕ i p = φ 2V p q∈Zp λS t p,q (T q -T p )
φp,q e i p,q [START_REF] Gillman | Third order model of thermal conductivity for random polydisperse particulate materials using well-resolved statistical descriptions from tomography[END_REF] where V p is the volume of the particle, Z p is the set of particles linked to the particle p and e i p,q is the component of the inter-particle normal vector corresponding to i direction. In discrete simulations, due to the limited size and the randomness of the set of inter-particular contacts associated to each DE, the heat flux field is always heterogeneous even if it is theoretically homogeneous. In a first study, the case of a homogeneous media is addressed in order to determine the level of dispersion due to fluctuations at the particle scale. For this purpose, the same thermal conditions as described in Subsection 2.2 are imposed with a thermal conductivity set to 33 W/(mK). Thus, for a temperature difference ∆T equal to 10 • C and a length L = 10cm, we should obtain an average heat flux equal to 3,300 W/m 2 according to the equation below:

ϕ = λ∆T L ( 21
)
For a set of particulate packings ranging from 100,000 to 500,000 DE, the process is carried out and the heat flux density field is determined in order to evaluate the CoV corresponding to each test.Whatever the density of particles, a CoV about 27 % is found which describes a high level of dispersion. As a solution, we propose to adapt the Halo approach introduced by Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D discrete element method simulation : Validation and application to flax/bio based epoxy composite[END_REF] in the context of the stress field determination to evaluate the heat flux (Figure 9). More precisely, at the mesoscopic scale, a Halo of spherical shape is introduced for every DE. The center of DE is supposed to also be that of its Halo. This latter has a radius R Halo and contains a pre-defined number of DE. The heat flux of a particle p is evaluated at the scale of Ω p the Halo volume in order to take into account the contributions of the particles in the neighborhood of p. The expression of the heat flux applied to p is then:

ϕ i p = 1 2Ω p r∈Ωp q∈Zr ϕ r,q e i r,q (22)
with

Ω p = 1 φ r∈Ωp V r ( 23
)
where Z r is the set of particles linked to the particle r and V r is the volume of particle r. One of the important questions deals with the choice of the suitable size of Halo. In this context, for a large spectrum of Halo-DE radius ratios, CoV is evaluated for two particulate packings of 100,000 and 500,000 DE. Results illustrated in Figure 10 From this scope of results, we can see that whatever the Halo size, mean values are very close to the analytical one. One can notice a convergence of the maximum and minimum values to the analytical one. The lower fluctuations correspond to the larger Halo size. In addition, the normalized heat flux is illustrated in Figure 12 for different Halo radii. We can notice that the level of dispersion reduces for a larger Halo radius.

Single inclusion problem

We now aim to study the validation of the proposed approach in the context of a continuous and heterogeneous material. To attend this objective, we consider the example of a biphasic material modeled by a cubic pattern composed of a single inclusion of spherical shape embedded in a matrix. We consider a cubic pattern of length L = 10cm and a spherical inclusion of radius a = L 3 (Figure 13). Thus, the volume fraction of inclusion is equal to 15.51%.

ETC

The numerical test is performed using a particulate packing composed of 150,000 particles which is a good compromise between accuracy and computational cost. λ i and λ m are respectively the thermal conductivities of the inclusion and matrix phases. The case where two particles in contact are not located in the same phase, which typically corresponds to contacts located at the matrix-inclusion interface Γ, is treated specifically. From a numerical standpoint, this interface is supposed to be perfect and without thermal barrier. In this context, the thermal conductivity associated to the contact zone between both particles is averaged. In this work, we consider the following inverse average:

λ Γ = 2λ i λ m λ i + λ m ( 24
)
This choice is justified by preliminary studies which showed that the arithmetic mean overestimates the ETC [START_REF] Haddad | Application of the discrete element method to study heat transfer by conduction in particulate composite materials[END_REF]. The specific heat is assumed to be 0.9 J/(K.kg) for both phases but this has little importance as long as this section focuses on the stationary state results. c λ = λ i λ m refers to the contrast of thermal properties between the matrix and the inclusion. Several numerical tests are performed for property contrasts varying from 0.01 to 100. Comparisons, in terms of ETC, with the FEM and the analytical Maxwell model are performed and presented in Figure 14. Regardless of the contrast, less or greater than 1, the results obtained by the DEM are close to those obtained by the FEM and the Maxwell model. The maximum relative error is in the order of 0.65%. This highlights the ability of the DEM to estimate the ETC in this context.

Temperature and heat flux fields

The thermal conductivities of the matrix and the inclusion are now fixed respectively to 33 and 165 W/(mK), so that the contrast of properties c λ is equal to 5. For the same configuration and thermal conditions, comparisons with FEM in terms of temperature and heat flux are set up. For information purposes, FEM calculations are performed using a structured mesh composed of 100,000 4-node tetrahedral elements. Figure 15 illustrates the temperature field obtained by (a) the DEM and (b) the FEM. One can assume that both methods represent a quasi-identical temperature field. Figure 15-c illustrates the relative differences on temperature, based on linear interpolation functions, with respect to the FEM calculations. Numerical comparisons exhibit a quite good agreement between the DEM and the FEM calculations with relative differences less than 1%. Heat flux is also studied using Equation ( 22) with a Halo-DE radius ratio of 11 which corresponds to a 2% indicator level according to Figure 10. Figure 16 illustrates the heat flux fields obtained by the DEM (a) and the FEM (b). From a qualitative standpoint, they show that the results obtained by the DEM are in good agreement with the FEM. Please notice that the results presented in this section are taken in the XY cutting plane with Z = 0.5L. From a quantitative standpoint, the temperature and heat flux values extracted at positions A(0.5L,0.5L,0.5L), B(0.08L,0.5L,0.5L) and C(0.5L,0.92L,0.5L) are very close (Table 2). The relative differences on heat flux with respect to FEM calculation are presented in Figure 16-c with a maximum relative difference of 19.6%. This one is due to the approximation of heat flux at the inclusion/matrix level and adiabatic boundary conditions. The results obtained show that DEM offers the possibility of estimating the heat flux field in a heterogeneous continuous media. According to these results, we can estimate that a particulate packing composed of 150,000 DE is suitable to represent the studied single-inclusion composite. Therefore, we estimate that a number of particles close to 23,000 DE is sufficient to represent a spherical inclusion.

Application to the case of a multi-inclusion composite material

This part focuses on a multi-inclusion composite material composed of an epoxy resin reinforced with ceramic fillers which was previously studied in [START_REF] Wong | Thermal conductivity, elastic modulus, and coefficient of thermal expansion of polymer composites filled with ceramic particles for electronic packaging[END_REF]. Thermal conductivity values were determined using a utility software program at a stabilized mean sample temperature of 70 • C. Thus, ETC of composite was evaluated for a range of volume fraction of inclusions f v from 0 to 50 % using different kinds of fillers. In this application, we are interested in fillers of silica with spherical shape. We should precise that the fillers of silica studied by Wong et al. [START_REF] Wong | Thermal conductivity, elastic modulus, and coefficient of thermal expansion of polymer composites filled with ceramic particles for electronic packaging[END_REF] were spherical polydisperse inclusions. For information purposes, thermal conductivity of epoxy is 0.195 W/(mK) and that of silica is 1.5 W/(mK) for a contrast of conductivities of 7.7. In a first step, we aim to compare the data obtained experimentally in term of ETC for different f v to results provided by analytical and numerical methods. In a second step, numerical simulations are carried out to determine the heat flux distribution for a defined volume fraction of silica.

ETC

We first aim at evaluating ETC using the DEM and comparing results with other analytical, numerical and experimental values. We consider the following procedure. In a first step, a random set of dilute spherical inclusions is generated using LSA [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF]. We limit our investigations to monodisperse inclusions and set the scale ratio between the length of the cubic pattern and the diameter of inclusions to 6. Thus, the number of fillers is directly governed by the targeted volume fraction. For information purpose, the length of the pattern L is set to 10cm but this has no influence on our results. The second steps consists in generating a dense enough particulate system to model the continuous medium. Based on the findings of Section 4, 23,000 DE are required to represent a spherical inclusion. Therefore, we can suppose that 8,000,000 DE can be considered enough to represent the medium where each spherical inclusion is modeled by about 20,000 DE. In a third step, the same thermal conditions are those cited in Subsection 2.2 are imposed. Finally, at the stationary state, the ETC is determined using the heat flux extracted from the numerical simulations and averaged. This process is carried out for a range of f v from 10 to 50 % considering the thermal properties of the above-mentioned epoxy/silica composite. Figure 17 illustrates the results of the discrete approach with those obtained experimentally by Wong et al. [START_REF] Wong | Thermal conductivity, elastic modulus, and coefficient of thermal expansion of polymer composites filled with ceramic particles for electronic packaging[END_REF] and compared with analytical results given by Maxwell and third-order models, and numerical ones obtained by FFT calculations using the Eyre-Milton scheme [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] and a regular grid composed of 134 millions of voxels (512×512×512 resolution). Relative errors between experimental and DE results are, at worse, close to 7 %. However, we hypothesize two explanations. First, the interface is assumed perfect without thermal barrier in our numerical simulations. Second, the polydispersity is not respected in numerical tests which could affect our results. Besides, relative differences with respect to the values given by the FFT-based approach are less than 0.5 % which highlights quite good adequation between DE and FFT approaches. From this scope of results, one can conclude that the DEM is able to predict the ETC of heterogeneous continuous media even with a high volume fraction of inclusions.

Temperature and heat flux fields

In this study, we consider the same configuration as previously but we limit our investigations to a volume fraction of silica of 20 % which corresponds to 83 spherical inclusions. For comparison purposes, numerical simulations are performed with DEM and FEM. In discrete approach, we handle the same particulate packing of 8,000,000 DE. In the Finite Element (FE) simulations, we consider a structured mesh composed of about 2,300,000 4-node tetrahedral elements to represent the multi-inclusion composite. At the stationary state, the heat flux is determined by both approaches and results in the YZ plane with X = 0.5L are extracted and illustrated in Figure 18. For information purposes, the Halo-DE radius ratio is again set to 11 in discrete simulation. The results at positions A(0.5L,0.6L,0.6L), B(0.5L,0.02L,0.5L), C(0.5L,0.3L,0.55L) and D(0.5L,0.82L,0.87L) in terms of temperature and heat flux are reported in Table 3. This scope of results exhibits a quite good adequation between FE and DE approaches. Thus, relative errors with respect to the values given by the FEM are less than 3 % in terms of heat flux and 0.6 % for temperature values. This conclusion exhibits the ability of the DEM to simulate the heat conduction throughout such material at the steady state.

Conclusion

The numerical approach presented in this work showed the potential of the DEM to model heat transfer by conduction in a heterogeneous continuous media. In a first step, an approach based on the DEM was developed and validated in the case of a homogeneous continuum media. In a second step, the discrete approach was applied to the context of a single inclusion composite. Some comparisons were done with the FEM in terms of ETC and heat flux field. In a final step, a multi-inclusion composite, namely an epoxy resin filled with silica inclusions was modeled. DEM results were validated with experimental data found in the literature, the FFT-based approach and the FEM as function of volume fraction of silica. In all studied configurations, DEM showed its ability to model the heat transfer by conduction and to predict correctly even in heterogeneous media with a high fraction of inclusions its ETC. Furthermore, two original contributions can be highlighted. First, a calibration process was introduced to determine the transmission surfaces. Such an approach ensures the equality between local and global conductivities as well as the conservation of mass while avoiding the costly step of Voronoï tessellation. Secondly, due to local fluctuations inherent to the DEM, we investigated a method to better evaluate the heat flux. Thus, the Halo approach was introduced to evaluate the heat flux at the scale of DE taking into consideration the contribution of DE located in its neighborhood. Results exhibited the capability of this concept. These fundamental conclusions pave the way of our future works. We look to combine Fourier and Fick law in order to model heat and mass transfer in heterogeneous media. Our ultimate goal is to provide an efficient and reliable numerical tool to simulate swelling and shrinkage mechanisms by adsorption and desorption of water. Furthermore, this numerical tool will be able to model thermal and hydric stresses and local damages in heterogeneous media.

Fig. 1 :

 1 Fig. 1: Heat transfer by conduction at the contact scale

Fig. 2 :

 2 Fig. 2: Numerical test for the calibration process

Fig. 3 :

 3 Fig. 3: Variation of Ct coefficient as function of DE number and geometry shape

Fig. 4 :

 4 Fig. 4: Thermal conditions applied to the cubic pattern

Fig. 5 :

 5 Fig. 5: Comparison between analytic and discrete model solutions at several times for temperature as function of position y/L (a) and relative errors (b)

Fig. 6 :

 6 Fig. 6: Thermal conditions applied to a spherical sample

Fig. 7 :

 7 Fig. 7: Evolution of the temperature at the center of the sphere.

Figures 8 -

 8 Figures 8-a and 8-b illustrate respectively the temperature fields at time t = 10s of DEM model with 100000

Fig. 8 :

 8 Fig. 8: Temperature profile of DEM model at time t = 10s with 100000 DE (a) and 1 million DE (b) and absolute difference between the DEM with 100000 DE and analytical solution (c)

Fig. 9 :

 9 Fig. 9: Illustration of the Halo approach at the mesoscopic scale

  exhibit that the level of dispersion is reduced for a higher Halo radius but is not influenced by the global density of particles. The user can easily determine the Halo radius corresponding to an expected heat flux distribution. A first indicator level at 5% is reached for a Halo-DE radius ratio close to 4. A second indicator at 2 % is reached for a Halo radius eleven times the radius of the DE. Besides, maximum, minimum and mean values in the discrete simulations are showed in Figure11compared with the expected value of heat flux.

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Evolution of CoV as function of Halo size

Fig. 12 :

 12 Fig. 12: Density functions of heat flux normalized by the average heat flux using different R Halo

Fig. 13 :

 13 Fig. 13: Single-inclusion composite

(a) c λ < 1 (b) c λ > 1 Fig. 14 :

 1114 Fig. 14: Single inclusion problem: non-dimensional thermal conductivity as a function of the contrast of properties

Fig. 15 :

 15 Fig. 15: Single inclusion problem: temperature profile of DEM (a) and FEM (b) models and relative differences (c)

Fig. 16 :

 16 Fig. 16: Single inclusion problem: Heat flux profile [W/m 2 ] of DEM (a) and FEM (b) models and relative differences (c)

Fig. 17 :

 17 Fig. 17: ETC (W/(mK)) as a function of the volume fraction of silica

Fig. 18 :Tab. 3 :

 183 Fig. 18: Heat flux (W/m 2 ) in the case of multi-inclusion composite using DE and FE approaches

Funding

The authors would like to gratefully acknowledge the European Union for the financial support under the INTERREG V France-Wallonie-Vlaanderen Program FWV N • 1.1.22.

Conflicts of interest

None declared.


Aiwen Lai
email: aiwenlai@xmu.edu.cn

Sébastien Lahaye
email: sebastien.lahaye@univ-angers.fr

Zhiwu Li

Initial-State Detectability and Initial-State Opacity of Unambiguous Weighted Automata

Keywords: Weighted automaton, initial state estimation, initial-state detectability, initial-state opacity

In this paper, we investigate the verification problem of initial-state detectability (I-detectability) and initial-state opacity (Iopacity) in discrete event systems modeled by unambiguous weighted automata. An I-observer is constructed so as to derive necessary and sufficient conditions for checking strong I-detectability, weak I-detectability, and I-opacity, with exponential complexity. In addition, an approach based on diagnosability analysis is proposed for verifying strong I-detectability. Compared with an I-observer-based approach, the diagnosability-based approach has a lower complexity, and in the case where all the unobservable events in an unambiguous weighted automaton are represented by a unique symbol, the diagnosability-based approach has polynomial complexity.

1 Introduction I-detectability and I-opacity problems are tightly related to state estimation that has been extensively studied in discrete event systems (DESs) framework. In general, state estimation can be divided into two categories: current-state estimation, e.g., [START_REF] Shu | Generalized detectability for discrete event systems[END_REF], [START_REF] Keroglou | Verification of detectability in probabilistic finite automata[END_REF], [START_REF] Lai | State estimation of max-plus automata with unobservable events[END_REF] and initial-state estimation, e.g., [START_REF] Shu | I-detectability of discrete-event systems[END_REF], [START_REF] Saboori | Verification of initial-state opacity in security applications of discrete event systems[END_REF], [START_REF] Yin | Initial-state detectability of stochastic discrete-event systems with probabilistic sensor failures[END_REF]. In this paper, we focus on the problem of initial-state estimation and its applications to I-detectability and I-opacity. [START_REF] Shu | I-detectability of discrete-event systems[END_REF] formulate the initial-state estimation problem of non-deterministic finite state automata (NFAs) as I-detectability. It characterizes whether the initial state of an NFA can be uniquely determined after a finite number of labels have been observed. An exponential (resp. a polynomial-time) complexity algorithm is introduced for checking weak (resp. strong) I-detectability based on the construction of an I-observer (resp. I-detector). [START_REF] Yin | Initial-state detectability of stochastic discrete-event systems with probabilistic sensor failures[END_REF] introduces the notion of stochastic initial-state detectability (SI-detectability) by considering the probabilistic sensor failures. SIdetectability characterizes the convergence of the probability for determining the initial-state of a probabilistic finite state automaton. An algorithm with PSPACEcomplete complexity is proposed for verifying the SIdetectability.

A system is said to be initial-state opaque if the set of initial states estimated by an intruder for any observation contains at least one element that does not belong to the secret. In [START_REF] Saboori | Verification of initial-state opacity in security applications of discrete event systems[END_REF], an initial-state estimator is constructed, which can be used to verify I-opacity of an NFA for both invariant-secret and varying-secret. The initial-state estimator has up to 2 n 2 states, where n is the number of states in the NFA. In addition, when the secret is fixed, the complexity of verifying I-opacity is reduced to O(4 n ) by introducing the notion of verifier. Recently, [START_REF] Wu | Comparative analysis of related notions of opacity in centralized and coordinated architectures[END_REF] prove that the initial state of an NFA can be estimated by the observer of its reverse automaton, and as a result, the verification complexity of I-opacity is reduced to O(2 n ). More work on (initial-state) detectability or opacity can be found in [START_REF] Yin | A new approach for the verification of infinite-step and k-step opacity using two-way observers[END_REF], [START_REF] Shu | Delayed detectability of discrete event systems[END_REF], [START_REF] Zhang | The problem of determining the weak (periodic) detectability of discrete event systems is pspace-complete[END_REF], Masopust & Yin (2019a), [START_REF] Zhang | Opacity of nondeterministic transition systems: A (bi) simulation relation approach[END_REF], [START_REF] Yin | Infinite-step opacity and k-step opacity of stochastic discrete-event systems[END_REF], [START_REF] Masopust | Deciding detectability for labeled petri nets[END_REF], [START_REF] Bryans | Modelling opacity using Petri nets[END_REF], [START_REF] Sasi | Detectability of networked discrete event systems[END_REF], [START_REF] Ji | Enforcing opacity by insertion functions under multiple energy constraints[END_REF].

Weighted automata (WAs) represent a well studied class of DES models [START_REF] Gaubert | Performance evaluation of (max,+) automata[END_REF]. Unlike classical automata, transitions in WAs carry weights belonging to a semiring. The weight associated with a transition can model, e.g., the cost, the energy, the time needed for executing the transition. WAs have spurred much interest in Computer Science due to their elegant and sound algebraic framework as well as their relevance in practical applications, such as natural language processing, speech recognition and image compression [START_REF] Droste | Handbook of weighted automata[END_REF]. From a different perspective, in this paper, we investigate I-detectability and I-opacity verification problems for WAs, which are core problems with the control community. So far, we have found a formal approach to verify these properties for an important class of WAs, namely, unambiguous weighted automata (UWAs), where no two or more paths are labeled by a given string leading to the same state.

To explore the modeling power of UWAs, it is worth mentioning that the relation between max-plus automata (MPAs) and timed Petri nets (TPNs) has been investigated. From a safe TPN under the preselection policy, it is possible to derive an MPA with the same timed behaviour [START_REF] Gaubert | Modeling and analysis of timed petri nets using heaps of pieces[END_REF][START_REF] Lahaye | Compositions of (max,+) automata[END_REF]. If the race policy is considered, it has been shown that bounded TPNs can be represented by deterministic max-plus automata (DMPAs) [START_REF] Komenda | Determinization of timed petri nets behaviors[END_REF][START_REF] Triska | Behaviour equivalent max-plus automata for a class of timed petri nets[END_REF]. Note that DMPAs constitute a subclass of unambiguous MPAs and the proposed I-detectability and I-opacity approach can be applied. It should also be noted that determinization procedures can be used to transform MPAs into DMPAs with the same timed behavior [START_REF] Gaubert | Performance evaluation of (max,+) automata[END_REF][START_REF] Mohri | Finite-state transducers in language and speech processing[END_REF][START_REF] Kirsten | A burnside approach to the termination of mohri's algorithm for polynomially ambiguous min-plus-automata[END_REF], Lahaye et al. 2020). Using such a transformation, a fairly large class of MPAs can be considered, that is, at least the polynomially ambiguous MPAs having the clones property.

Because of the influence of the transition weights, the I-detectability and I-opacity of a UWA are different from that of its support, i.e., the corresponding logical automaton obtained by removing all the weights in the UWA. For instance, considering the UWA G in Fig. 1, its support is not strongly I-detectable. On the other hand, taking into account quantitative information, it is shown later that G is strongly I-detectable. Besides, if we assume that the set of secret states is Q s = {1, 3}, it can be shown that G is not initial-state opaque while its support is initial-state opaque with respect to Q s .

The above aspects motivate the work in this paper, i.e., the verification of I-detectability and I-opacity for UWAs. The main contributions are as follows. 1) We extend the notion of I-detectability and I-opacity from logical DESs to the framework of WAs. Two types of I-detectability, namely strong I-detectability and weak I-detectability are defined. 2) Given a UWA, a formal procedure is first proposed to construct its initialstate estimator (called I-observer). Then, necessary and sufficient conditions with exponential complexity are developed for verifying I-detectability and I-opacity of the studied UWA based on the constructed I-observer.

3) An approach based on diagnosablility analysis is presented for checking strong I-detectability of a UWA. In the case that the unobservable events are not distinguished by an external agent, i.e., all unobservable events in a UWA are represented by a single symbol, the approach is of polynomial complexity. This paper is organized as follows. Section 2 recalls some basics of WAs. Section 3 formulates the notion of I-detectability and I-opacity. In Section 4.1, the Iobserver based approach is proposed to verify strong and weak I-detectability of a UWA. Section 4.2 presents a diagnosability-based approach for verifying strong Idetectability of a UWA. In Section 5, a necessary and sufficient condition is presented for checking I-opacity of a UWA. Finally, conclusions are drawn in Section 6.

Preliminaries

In this section we recall some basics of weighted automata [START_REF] Droste | Handbook of weighted automata[END_REF], where transitions carry weights belonging to a semiring S = (D, ⊕, ⊗, ε, e), where ε (resp. e) denotes the neutral element for addition ⊕ (resp. multiplication ⊗). Typical examples are the tropical semiring (R ∪ {+∞} , min, +, +∞, 0) and maxplus semiring (R ∪ {-∞} , max, +, -∞, 0).

Let E be an alphabet, i.e., a non-empty set of labels, and E * be the set of all the finite strings over E including the empty string λ. The set of matrices with m rows and n columns over semiring S is denoted by S m×n .

Definition 1 A weighted automaton over a semiring

S = (D, ⊕, ⊗, ε, e) is equivalently defined by G = (Q, E, α, µ) or G = (Q, E, t, Q i , )
, where

• Q and E are respectively a non-empty finite set of states and an alphabet; • α ∈ S 1×|Q| is a row vector specifying the initial weights. A state q ∈ Q is said to be an initial state iff α q = ε, where α q is the initial weight of q. We denote by Q i the set of initial states, and : Q i → D is the function of initial weights (q) α q for q ∈ Q i ; • µ: E → S |Q|×|Q| is a morphism representing the state transitions given by the family of matrices µ (a) ∈ S |Q|×|Q| , a ∈ E. More precisely, we have µ(a) qq = ε 1 if, and only if, there is a transition labeled by a with a weight of µ(a) qq from state q to q . For any string

ω = e 1 • • • e k ∈ E * , we have µ(ω) = µ(e 1 ) ⊗ µ(e 2 ) ⊗ • • • ⊗ µ(e k ). t : Q × E × Q → D
is the transition function with t(q, a, q ) µ(a) qq , and for string ω = e 1 e 2 • • • e n ∈ E * , t(q, ω, q ) = µ(ω) qq .

Definition 2 Given a WA G, a path of length k is defined as a sequence of transitions π = (q 0 , e 1 , q 1 )(q 1 , e 2 , q 2 ) • • • (q k-1 , e k , q k ), where q j ∈ Q, j = 0, . . . , k, e j ∈ E and µ(e j ) qj-1qj = ε for j = 1, . . . , k.

Such a path π leads from state q 0 to q k , and is labeled by e 1 e 2 • • • e k ∈ E * . Besides, π is said to be a circuit if q 0 coincides with q k . We use p ω q to represent the set of paths labeled by string ω ∈ E * from state p to q. For P, R ⊆ Q, we denote by P ω R the union of p ω q for all p ∈ P and q ∈ R.

Definition 3 A WA G is said to be unambiguous if for all q ∈ Q, for all ω ∈ E * , |Q i ω {q} | ≤ 1.
In simple words, the unambiguity implies that for any state q of G and any string ω in E * , there is at most one path labeled by ω leading from an initial state to q.

Remark 1 If Q i = Q, a WA is unambiguous iff any state has no two or more input transitions labeled by the same symbol. In the case of Q i = Q, the above characterization provides only a sufficient condition for unambiguity. Besides, a WA is said to be deterministic if it has a unique initial state and from any state, no two or more output transitions are labeled by the same symbol. Determinism implies unambiguity, but the reverse is not true.

2

Example 1 Fig. 1 depicts a UWA G = (Q, E, α, µ), where Q = {1, 2, 3, 4, 5, 6, 7}, E = {u, b, c, d}, µ(u) 1,2 = 1, µ(u) 4,5 = 2, µ(u) 5,6 = 3, µ(b) 2,3 = 6, µ(b) 6,3 = 4, µ(c) 3,3 = 2, µ(c) 7,7 = 2, µ (
= (q 0 , e 1 , q 1 )(q 1 , e 2 , q 2 ) • • • (q k-1 , e k , q k ) of UWA G with q 0 ∈ Q i , the weighted sequence σ(π) ∈ (E × D) * generated by π is defined as: σ(π) = (e 1 , τ 1 )(e 2 , τ 2 ) • • • (e k , τ k ) where τ 1 = α q0 ⊗ µ(e 1 ) q0q1 , τ j = τ j-1 ⊗ µ(e j ) qj-1qj for j = 2, . . . , k.
1 We assume that the states of G are ordered by positive integers, and by a slight abuse of notation, αq (resp. µ(a) qq ) is used to denote the q th element of α (resp. the element in the q th row and q th column of matrix µ(a)). Here, σ(π) specifies a sequence of labels and their occurrence weights. We use q 0 σ(π) q k to represent that weighted sequence σ(π) leads from q 0 to q k . In the rest of this paper, we focus on semirings where ⊗ represents the usual addition, which fits with WAs where weights represent the durations or consumed energy of state transitions. In other words, the weights along an evolution of a WA are interpreted as time or energy accumulated by addition, and an external agent observes this value in addition to events occurrences.

Definition 5 Given a UWA G = (Q, E, α, µ), the generated weighted language L(G) of G is defined as:

L(G) = {y ∈ (E × D) * | ∃q ∈ Q, ∃ω ∈ E * , ∃π ∈ Q i ω {q} : σ(π) = y}. (1)
For any sequence σ ∈ L(G), we denote by

lab f (σ) ∈ E the last label in σ. Example 2 Consider path π = (1, u, 2)(2, b, 3)(3, c, 3) in UWA G in Fig. 1. By Def. 4, we have σ(π) = (u, e ⊗ µ(u) 1,2 )(b, e ⊗ µ(u) 1,2 ⊗ µ(b) 2,3 )(c, e ⊗ µ(u) 1,2 ⊗ µ(b) 2,3 ⊗ µ(c) 3,3 ) = (u, e ⊗ 1)(b, e ⊗ 1 ⊗ 6)(d, e ⊗ 1 ⊗ 6 ⊗ 2) = (u, 1)(b, 7)(c, 9), and lab f (σ(π)) = c. 2
In this paper, alphabet E is partitioned into two disjoint parts: the unobservable part E uo and the observable part E o . The projection operator P : E * → E * o is defined as: P (λ) = λ, where λ represents the empty string; for each a ∈ E, ω ∈ E * , P (ωa) = P (ω)a if a ∈ E o , otherwise, P (ωa) = P (ω). We extend P :

E * → E * o to weighted sequences P : (E × D) * → (E o × D) * .
We denote by P (L(G)) the set of all observable weighted sequences for G, and by σ 1 σ 2 the concatenation of weighted sequences σ 1 and σ 2 . In addition, σ 2 is a prefix of σ 1 if there exists another sequence σ 3 such that σ 1 = σ 2 σ 3 .

Notion of I-Detectability and I-Opacity

As in [START_REF] Lai | Verification of detectability for unambiguous weighted automata[END_REF], without loss of generality, we restrict our attention to a WA G = (Q, E, α, µ), where the weights of all initial states are equal to e. We make the following assumptions on the studied WA G: (1) G is unambiguous; (2) G is deadlock free, that is, for any state of the system, there exists at least one output transition, i.e., (∀q ∈ Q)(∃a ∈ E, q ∈ Q)µ(a) qq = ε; (3) There is no circuit labeled only by unobservable labels in G.

(2) implies that the length of a generated weighted sequence becomes infinite as the system evolves indefinitely. (3) implies that the generated sequences of unobservable labels have finite length. Besides, as mentioned in the introductory section, (1) defines the broadest class of WAs for which we have been so far able to handle I-detectability and I-opacity verification.

Initial-State Detectability

Definition 6 Given a UWA G, the set of possible initial states after observing σ o ∈ P (L(G)) is defined as

I(σ o ) = {q ∈ Q i | ∃q ∈ Q, ∃σ ∈ L(G) : P (σ) = σ o , q σ q }. (2)
Lemma 1 Consider a non-empty observation σ o ∈ P (L(G)). Let

I (σ o ) = {q ∈ Q i | ∃q ∈ Q, ∃σ ∈ L(G) : P (σ) = σ o , q σ q , lab f (σ) = lab f (σ o )}. ( 3
) Then I(σ o ) = I (σ o ).
Proof: Let q i be an initial state with q i ∈ I (σ o ). According to Eqs. ( 2) and ( 3), q i ∈ I(σ o ). On the other hand, consider an initial state q 0 ∈ I(σ o ). By Eq. ( 2), there must exist a path π from q 0 ending with lab f (σ o ) or an unobservable label such that P (σ(π)) = σ o . If π ends with lab f (σ o ), then it is trivial that q 0 ∈ I (σ o ). When π ends with an unobservable label, we assume that it can be represented as π = (q 0 , v 1 e 1 , q 1 ) • • • (q n-1 , v n e n , q n )(q n , v n+1 , q n+1 ), where v 1 , . . . , v n+1 ∈ E * uo . Let π = (q 0 , v 1 e 1 , q 1 ) • • • (q n-1 , v n e n , q n ). Then π satisfies Eq.

(3). Hence, q 0 belongs to I (σ o ).

2

As in a logical DES [START_REF] Shu | I-detectability of discrete-event systems[END_REF], in this paper, a possible trajectory of a WA G is represented by an infinite sequence of (label, weight) pairs that G may generate. The set of all possible trajectories of G defines the ω-language L ω (G). For any sequence σ ∈ L ω (G), we denote by P re(σ) the set of all its prefixes.

Definition 7 (Strong (Weak) I-Detectability) Given a UWA G = (Q, E, α, µ), where the set of initial states Q i ⊆ Q is not empty, i.e., Q i = ∅, G is strongly (resp.
weakly) I-detectable with respect to projection P if, for all (resp. some) trajectories, the set of possible initial states shrinks to a singleton after a finite number of observations, i.e.,

(∃n ∈ N)(∀σ (resp. ∃σ) ∈ L ω (G))(∀σ ∈ P re(σ)) |P (σ )| > n ⇒ |I(P (σ ))| = 1.
From the above definition, we know that if a UWA is strongly I-detectable, then it is weakly I-detectable.

Initial-State Opacity

The generated weighted language of a UWA G = (Q, E, α, µ) from an initial state q i ∈ Q i is defined as

L(G, q i ) = {y ∈ (E × D) * | ∃q ∈ Q, ∃ω ∈ E * , ∃π ∈ q i ω q : σ(π) = y}. (4)
In addition, given a subset X ⊆ Q i , we define the weighted language generated from X as L(G, X) = qi∈X L(G, q i ). Due to unambiguity, if

Q i = Q 1 i ∪ Q 2 i , then the generated weighted language of G can be represented by L(G) = qi∈Qi L(G, q i ) = L(G, Q 1 i ) ∪ L(G, Q 2 i ), and P (L(G)) = P (L(G, Q 1 i )) ∪ P (L(G, Q 2 i )).
Given a UWA G = (Q, E, α, µ), an intruder can only observe the projection of the generated language, i.e., P (L(G)). Assume that the intruder has full knowledge of the structure of G, and that the secret is described by an arbitrary subset of Q.

Definition 8 (Initial-state opacity) Given a UWA G = (Q, E, α, µ), projection P : E * → E * o , and a set of secret states Q s ⊆ Q, G is initial-state opaque with respect to Q s and P , if for all q ∈ Q i ∩ Q s and for all σ ∈ L(G, q), there exist q ∈ Q i \ Q s and σ ∈ L(G, q ) such that P (σ ) = P (σ), that is,

P (L(G, Q i ∩ Q s )) ⊆ P (L(G, Q i \ Q s )), or equivalently, P (L(G)) = P (L(G, Q i \ Q s )).
Example 3 Consider again the UWA G in Fig. 1. We assume that the set of secret states is Q s = {3, 4}, and u is the only unobservable label. (1) Detectability analysis: Initially, G can be in state 1 and/or state 4. Once an observation is obtained, i.e., (b, 7), (b, 9) or (d, 4), we can uniquely determine the initial state, from which the system starts. Therefore, G is strongly I-detectable and weakly I-detectable. If the weight of transition (5, u, 6) is now changed from 3 to 1, then G is not strongly I-detectable. It can be found that for infinite weighted sequence (b, 7)(c, 9)(c, 11)(c, 13)

• • • , we have I((b, 7)(c, 9)(c, 11)(c, 13) • • • ) = {1, 4}.
(2) Opacity analysis: Consider the secret initial state 4 and σ = (u, 2)(d, 4) ∈ L(G, 4). No weighted sequence can be generated by G from non-secret initial state 1 such that its projection is equal to P (σ). Therefore G is not initialstate opaque with respect to secret Q s . 2

4 Verification of I-Detectability

I-observer for Strong and Weak I-Detectability

We show that the strong and weak I-detectability of a UWA G can be checked using an I-observer of an augmented automaton obtained from G. The construction of the I-observer extends the approach in Shu & Lin (2013) by using a weighted alphabet that captures the quantitative information associated to observations.

The I-observer G aug obs = (Q aug obs , E aug obs , δ aug obs , q aug i,obs ) of a UWA G = (Q, E, t, Q i , ) is constructed by Algorithm 1, which contains two steps.

Step 1 extends G to its augmented version G aug = (Q aug , E, t aug , Q aug i , aug ), in which each state q aug ∈ Q aug is a pair (q c , q i ), representing that state q c can be reached in G from initial state q i . Step 2 computes the observer G aug obs for the augmented automaton G aug by the approach that is first presented in our recent work [START_REF] Lai | Verification of detectability for unambiguous weighted automata[END_REF]. Observer G aug obs is a deterministic finite state automaton over a weighted alphabet E aug obs ⊆ E o × D. That is, G aug obs has only one initial state, and from a given state no two transitions of G aug obs are labeled by the same weighted label (a, t a ) ∈ E aug obs .

Algorithm 1 Construction of the I-observer of a UWA

Input: A UWA G = (Q, E, t, Q i , ). Output: An I-observer G aug obs = (Q aug
obs , E aug obs , δ aug obs , q aug i,obs ) of G. 1: Construct an augmented automaton G aug of G as follows:

• Q aug i = {(q, q) | q ∈ Q i } is the set of initial states; • aug : Q aug i → D
is the function of initial weights defined as: aug ((q, q)) = (q);

• Q aug 1 = {(q, q i ) | q ∈ Q, q i ∈ Q i }; • t aug : Q aug 1 × E × Q aug 1
→ D is transition function, where t aug ((q c , q i ), a, (q c , q i )) is defined as: t aug ((q c , q i ), a, (q c , q i )) = t(q c , a, q c ) if q i = q i ; otherwise, t aug ((q c , q i ), a, (q c , q i )) = ε;

• Let G aug = (Q aug , E, t aug , Q aug i , aug ) = Ac(Q aug 1 , E, t aug , Q aug i , aug ) 2 . 2:
Compute the observer G aug obs of G aug , i.e., the Iobserver of G as follows:

• q aug i,obs = Q aug i ; • E aug obs consists of all weighted labels (a, τ ) ∈ E o × (D \ {ε}) for which ∃q ∈ Q 1 ∪ Q aug i , ∃q ∈ Q aug , ∃v ∈ E * uo , s.t. t aug (q, va, q ) = τ , where Q 1 = {q ∈ Q aug | ∃q ∈ Q aug , ∃a ∈ E o : t aug (q , a, q) =
ε} is the set of states in G aug that have at least one input transition marked by an observable label; • δ aug obs : 2 Q aug ×E aug obs → 2 Q aug is the state transition function. δ aug obs (q aug obs , (a, τ )) is defined as:

δ aug obs (q aug obs , (a, τ )) = {q ∈ Q aug | ∃q ∈ q aug obs , ∃v ∈ E * uo : t aug (q, va, q ) = τ }.
• Let G aug obs = (Q aug obs , E aug obs , δ aug obs , q aug i,obs ) = Ac(2 Q aug , E aug obs , δ aug obs , q aug i,obs ).

Example 4 Consider again the UWA G in Fig. 1.

2 In this paper, we denote by the Ac(G) the automaton obtained by removing all the states that are not accessible as well as transitions associated with such states in G.

By applying Algorithm 1, we obtain the augmented automaton G aug and I-observer G aug obs visualized in Figs. 2 and3, respectively. 2 Let (E aug obs ) * be the set of all the finite strings over weighted alphabet E aug obs including (λ, e). The language generated by the I-observer is L(G aug obs ) = {ω ∈ (E aug obs ) * | ∃q aug obs ∈ Q aug obs : δ aug obs (q aug i,obs , ω) = q aug obs }, which is a subset of (E aug obs ) * , i.e., L(G aug obs ) ⊆ (E aug obs ) * . For any observed weighted sequence

/ 1 u / 6 b / 2 u / 2 c / 2 c / 4 b / 2 c / 2 d ( ) 1,1 ( ) 2,1 ( ) 3,1 ( ) 4, 4
σ o = (a 1 , τ 1 )(a 2 , τ 2 ) • • • (a n , τ n ) ∈ P (L(G)), we define σ elem o = (a 1 , τ 1 )(a 2 , τ 2 -τ 1 ) • • • (a n , τ n -τ n-1
) to denote its equivalent notation in (E aug obs ) * . Note that τ k -τ k-1 represents the weight for the elementary transition according to a k for k = 2, 3, . . . , n. We denote by P (L(G))

elem the equivalent notation of P (L(G)), that is,

P (L(G)) elem = {σ ∈ (E o × D) * | ∃σ o ∈ P (L(G)) : σ elem o = σ}.
The following lemma states the relationship between the languages generated by G and G aug obs . The proof follows immediately from the construction process of G aug obs , and hence here it is omitted.

Lemma 2 The language L(G aug obs ) generated by Iobserver G aug obs coincides with P (L(G)) elem , that is,

L(G aug obs ) = P (L(G)) elem .
Lemma 3 The possible initial states of UWA G after observing a non-empty weighted sequence

σ o = (a 1 , τ 1 )(a 2 , τ 2 ) • • • (a k , τ k ) ∈ P (L(G)
) is given by

I(σ o ) = {q i ∈ Q i | ∃q c ∈ Q : (q c , q i ) ∈ δ aug obs (q aug i,obs , (a 1 , τ 1 )(a 2 , τ 2 -τ 1 ) • • • (a k , τ k -τ k-1 ))}.
Proof:

I(σ o ) = I (σ o ) (by Lemma 1) = {q i ∈ Q i | ∃q c ∈ Q, ∃σ ∈ L(G) : P (σ) = σ o , q i σ q c , lab f (σ) = lab f (σ o )} = {q i ∈ Q i | ∃q c ∈ Q, ∃σ ∈ L(G aug ) : P (σ) = σ o , (q i , q i ) σ (q c , q i ), lab f (σ) = lab f (σ o )} (by step 1 in Algorithm 1) = {q i ∈ Q i | ∃q c ∈ Q : (q c , q i ) ∈ δ aug obs (q aug i,obs , (a 1 , τ 1 )(a 2 , τ 2 -τ 1 ) • • • (a k , τ k -τ k-1
))} (by Proposition 1 in [START_REF] Lai | Verification of detectability for unambiguous weighted automata[END_REF]) 2

Now we introduce the criteria for checking strong Idetectability and weak I-detectability for UWA G based on its I-observer G aug obs . We denote by S ci the set of all elementary circuits of G aug obs as:

S ci = {(q aug obs , s) ∈ Q aug obs × (E aug obs ) * | δ aug obs (q aug obs , s) = q aug obs ∧ |s| ≥ 1 ∧ [∀s ∈ P re(s) s.t. s = s ∧ |s | ≥ 1 : δ aug obs (q aug obs , s ) = q aug obs ]}.
Besides, we use q aug obs ∈ (q aug obs , s) to represent that node q aug obs is visited by circuit (q aug obs , s). Let Q aug I,obs be a subset of Q aug obs consisting of states with the same second element (initial state of G), that is,

Q aug I,obs = {q aug obs ∈ Q aug obs | ∃q ∈ Q i , ∀(q c , q i ) ∈ q aug obs : q i = q}.
Theorem 1 (Criterion for Checking Strong I-Detectability) A UWA G is strongly I-detectable with respect to projection P iff any node belonging to an elementary circuit of G aug obs is in Q aug I,obs . Formally, ∀q aug obs ∈ Q aug obs , if ∃(q aug obs , s) ∈ S ci such that q aug obs ∈ (q aug obs , s), then q aug obs ∈ Q aug I,obs .

Proof: (If) By Algorithm 1, we know that once G reaches a state belonging to Q aug I,obs , it will stay in Q aug I,obs forever. Suppose ∀q aug obs ∈ Q aug obs , if ∃(q aug obs , s) ∈ S ci such that q aug obs ∈ (q aug obs , s), then q aug obs ∈ Q aug I,obs . Then, Iobserver Q aug obs will eventually reach a state in Q aug I,obs after a finite number of observations for all trajectories of G. Hence, by combining Lemma 3, the set of possible initial states is a singleton for all possible continuations, that is, G is strongly I-detectable. (Only If) Assume ∃q aug obs ∈ Q aug obs , ∃(q aug obs , s) ∈ S ci such that q aug obs ∈ (q aug obs , s) and q aug obs ∈ Q aug obs \ Q aug I,obs . Then a possible evolution of G can iterate indefinitely circuit (q aug obs , s). That is, we cannot determine the initial state of G for such a trajectory forever. Hence, G is not strongly I-detectable.

2

Theorem 2 (Criterion for Checking Weak I-Detectability) A UWA G is weakly I-detectable with respect to P iff Q aug I,obs = ∅.

Proof: (If) If Q aug I,obs is not empty, then I-observer Q aug obs can reach a state in Q aug I,obs after a finite number of observations for some trajectories of G. In this case, we can determine the only possible initial state of G for trajectories corresponding to continuations of these observations. Hence, G is weakly I-detectable. (Only If) If Q aug I,obs = ∅, then the set of possible initial states does not shrink to a singleton whatever the trajectory is, i.e., G is not weakly I-detectable.

2

Example 5 Consider the I-observer G aug obs in Fig. 3 of UWA G in Fig. 1. There are four states in G aug obs , i.e., {(1, 1), (4, 4)}, {(3, 1)}, {(3, 4)}, {(7, 4)}, and we have

Q aug I,obs = {{(3, 1)} , {(3, 4)} , {(7, 4)}}.
It can be checked that all elementary circuits have all their nodes in Q aug I,obs . Therefore, by Theorem 1, G is strongly I-detectable (hence weakly I-detectable).

2

Remark 2 For checking strong I-detectability of UWA G = (Q, E, α, µ), we have to find out all states that belong to an elementary circuit in its I-observer G aug obs = (Q aug obs , E aug obs , δ aug obs , q aug i,obs ). This can be done by finding all the strongly connected components in G aug obs , whose complexity is linear in the number of states and transitions of G We present an approach that has lower complexity (compared with the I-observer-based method), but is still exponential, to verify the strong I-detectability of a UWA G. We show that this verification can be done by means of diagnosability analysis of a new UWA G obtained from G. To perform the diagnosability analysis, we extend the algorithm in [START_REF] Jiang | A polynomial algorithm for testing diagnosability of discrete-event systems[END_REF] to build an automaton denoted by G d from G with a weighted alphabet.

Definition 9 Given a UWA G = (Q, E, t, Q i , )
, where

Q i = {q 1 , q 2 , . . . , q n } is the set of n initial states, a modified UWA G = (Q , E , t , Q i , ) obtained from G is defined by adding a new state q w / ∈ Q and |Q i | new labels f 1 , f 2 , . . . , f |Qi| ∩ E = ∅ as follows. • Q = Q ∪ {q w } is the finite set of states; • E = E ∪ f 1 , f 2 , . . . , f |Qi| ; • t : Q × E × Q → D
is the transition function defined as: t (q w , f i , q i ) = e for i = 1, 2, . . . , |Q i |; t (q, a, q ) = t (q, a, q ) if q, q ∈ Q and a ∈ E; otherwise, t (q, a, q ) = ε;

• Q i = {q w } is the set of initial states; • : Q i → D is the initial weights function: (q w ) = e.
In simple words, G is a UWA obtained from G by adding a new state q w to become the unique initial state, and by introducing new state transitions (q w , f i , q i ) with identity weights for i = 1, 2, . . . , |Q i |.

Example 6 Consider the UWA G in Fig. 1. The modified automaton G is depicted in Fig. 4.

2 1 / u e f / 1 2 / u e f / 2 3 / u 4 / b 6 / b 2 / d 2 / c 2 / c
Let f i , i = 1, 2, . . . , |Q i |, be unobservable labels belonging to |Q i | different types of faults, i.e., E f = E f 1 ∪ E f 2 ∪ • • • ∪ E f |Qi| = {f 1 } ∪ {f 2 } ∪ • • • ∪ f |Qi| ,
and suppose that all labels in E represent normal behavior. We denote by

E uo = E uo ∪ E f the set of unobservable labels in G . For any weighted sequence s = (e 1 , τ 1 )(e 2 , τ 2 ) • • • (e n , τ n ) ∈ (E × D) * or any logical sequence s = e 1 e 2 • • • e n ∈ (E )
* , we denote by lab i (s) = e 1 the first label in s. As in logical automata [START_REF] Sampath | Diagnosability of discrete-event systems[END_REF], the diagnosability of UWA G = (Q , E , t , Q i , ) can be defined as follows.

Definition 10 (Diagnosability) The UWA G = (Q , E , t , Q i , ) is diagnosable with respect to projection P : (E ) * → E * o and |Q i | different types of faults E f = f 1 , f 2 , . . . , f |Qi| if the following holds:

(∀f i ∈ E f )(∃n i ∈ N)(∀v = (f i , e)t ∈ L(G ), |t| ≥ n i ) ⇒ (∀ω ∈ L(G ), P (ω) = P (v))(lab i (ω) = f i ).
In simple words, diagnosability requires that when a fault label f i occurs in G , after a finite number of observations, one can detect its occurrence. The following lemma follows immediately from the construction process of G given in Def. 9, and the relationship between diagnosability and strong I-detectability.

Lemma 4 A UWA G is strongly I-detectable iff the modified UWA G given in Def. 9 is diagnosable with respect to P :

(E ) * → E * o and |Q i | different types of faults E f = f 1 , f 2 , . . . , f |Qi| .
Inspired by the work in [START_REF] Jiang | A polynomial algorithm for testing diagnosability of discrete-event systems[END_REF], we present an approach for determining if G is diagnosable with respect to a single fault type

E f i = {f i }, i ∈ {1, 2, . . . , |Q i |}, based on the construction of a finite state automaton G d = (Q d , E d , δ d , q d w ) over a weighted alphabet E d ⊆ E o × D, as defined in Algorithm 2. Algorithm 2 Construction of G d for E fi = {f i } Input: UWA G = (Q , E , t , Q i , ). Output: NFA G d = (Q d , E d , δ d , q d w )
. 1: Construct an NFA G o from G as follows:

• q o w = (q w , N );

• Q o = {(q, f ) | q ∈ Q 1 ∪ {q w } , f ∈ {N, F i }} where Q 1 = {q ∈ Q | ∃q ∈ Q , ∃a ∈ E o : t (q , a, q) = ε}
is the set of states in G that have at least one input transition marked by an observable label;

• E d consists of all weighted labels (a, τ ) ∈ E o × (D \ {ε}) for which ∃q ∈ Q 1 ∪ {q w }, ∃q ∈ Q , ∃v ∈ (E uo ) * , s.t. t (q, va, q ) = τ ; • δ o ⊆ Q o × E d × Q o is the set of state transitions. ((q, f ), (a, τ ), (q , f )) ∈ δ o iff ∃v ∈ (E uo ) * s.t. 1) t (q, va, q ) = τ ; 2) (f = f ) ∧ (|v| = 0 ∨ (|v| ≥ 1 ∧ lab i (v) = f i )) or (f = F i ) ∧ (|v| ≥ 1 ∧ lab i (v) = f i ). • Let G o = (Q o , E d , δ o , q o w ) = Ac(Q o , E d , δ o , q o w ). 2: Compute automaton G d = G o
G o , the parallel composition of G o with itself, as follows:

• q d w = (q o w , q o w ) = ((q w , N ), (q w , N )); Let S cd be the set of all elementary circuits of G d . Besides, for q d ∈ Q d and (q d , s) ∈ S cd , we use q d ∈ (q d , s) to represent that node q d is visited by circuit (q d , s).

• Q d = (q 1 o , q 2 o ) | q 1 o , q 2 o ∈ Q o ; • δ d ⊆ Q d × E d × Q d is the set of state transitions. ((p 1 o , p 2 o ), (a, τ ), (q 1 o , q 2 o )) ∈ δ d iff (p 1 o , (a, τ ), q 1 o ) ∈ δ o and (p 2 o , (a, τ ), q 2 o ) ∈ δ o . • Let G d = (Q d , E d , δ d , q d w ) = Ac(Q d , E d , δ d , q d w ).
) , 4 d ( ) , 2 c ( ) ,7 b ( ) ,9 b ( ) , 2 c ( ) , 2 c ( ) ( ) ( ) 0, , 0, N N ( ) ( ) ( ) 1 1 3, , 3, F F ( ) ( ) ( ) 7, , 7, N N ( ) ( ) ( ) 3
Theorem 3 UWA G is diagnosable with respect to fault type

E f i = {f i }, i ∈ {1, 2, . . . , |Q i |}, iff for every node ((q 1 d , f 1 ), (q 2 d , f 2 )) ∈ Q d that belongs to an elementary circuit in G d , we have f 1 = f 2 . Formally, ∀((q 1 i , f 1 ), (q 2 i , f 2 )) ∈ Q d , if ∃(q d , s) ∈ S cd such that ((q 1 i , f 1 ), (q 2 i , f 2 )) ∈ (q d , s), then f 1 = f 2 . Proof: (If)
We suppose that G is not diagnosable with respect to f i , i ∈ {1, 2, . . . , |Q i |}. According to Def. 10 on diagnosability, this implies that there exist two infinite weighted sequences in G such that they have equal projection, one contains f i and the other one does not contain f i . Since G d is finite, by considering the relationship between G d and G , there must exist an elementary circuit (q d , s), where the labels of all the first elements of all nodes are N , and the labels of all the second elements of all nodes are F i . Hence, ∃((

q 1 i , f 1 ), (q 2 i , f 2 )) ∈ Q d , ∃(q d , s) ∈ S cd such that ((q 1 i , f 1 ), (q 2 i , f 2 )) ∈ (q d , s) and f 1 = f 2 . (Only If) We assume that ∃((q 1 i , f 1 ), (q 2 i , f 2 )) ∈ Q d , ∃(q d , s) ∈ S cd such that ((q 1 i , f 1 ), (q 2 i , f 2 )) ∈ (q d , s
) and f 1 = f 2 . We can further assume that f 1 = N and f 2 = F i . According to the construction of G d , the above assumptions imply that for any ((q 1 j , f 1 j ), (q 2 j , f 2 j )) ∈ (q d , s), we have f 1 j = f 1 = N and f 2 j = f 2 = F i . Since elementary circuit (q d , s) can be repeated infinitely, by considering the relationship between G d and G , there are two infinite weighted sequences γ 1 and γ 2 in G with the same projection, γ 1 contains fault f i , and γ 2 does not contain f i . More precisely, γ 1 and γ 2 lead from the initial state q d w to the entry point q d of circuit (q d , s) and repeat this circuit infinitely. Assume that q d is reached by v, i.e., q d ∈ δ d (q d w , v). Then we have P (γ 1 ) elem = P (γ 2 ) elem = sv j , where j is a large enough integer. Therefore, by Def. 10 of diagnosability, G is not diagnosable.

2

Remark 3 Similar to Remark 2, the necessary and sufficient condition in Theorem 3 can be verified by finding all the strongly connected components in In this section, we show that the I-observer constructed in Section 4.1 can be used to check I-opacity of UWAs.

G d = (Q d , E d , δ d , q d ω ), whose complexity is O(|Q d | + |Q d | × |E d | × |Q d |). From Algorithm 2, the number of states in G o is bounded by 2 × (|Q| + 1) = 2|Q| + 2. Due to G d = G o G o , |Q d | is bounded by (2|Q| + 2) 2 = 4|Q| 2 +8|Q|+4.
Lemma 5 Given a UWA G = (Q, E, α, µ), projection P : E * → E * o , and a set of secret states Q s ⊆ Q, G is initial-state opaque with respect to Q s and P iff for all σ o ∈ P (L(G)), I(σ o ) Q s holds.

Proof:

(∀σ o ∈ P (L(G)))I(σ o ) Q s ⇔ (∀σ o ∈ P (L(G)))(∃q ∈ I(σ o ))q / ∈ Q s ⇔ (∀σ o ∈ P (L(G)))(∃q ∈ Q i \ Q s )(∃q ∈ Q) (∃σ ∈ L(G))P (σ ) = σ o ∧ q σ q
(by the definition of I(σ o ) in Equation ( 2)) ⇔ (∀σ o ∈ P (L(G)))(∃q

∈ Q i \ Q s )(∃σ ∈ L(G, q )) P (σ ) = σ o (according to Equation (4)) ⇔ (∀σ o ∈ P (L(G, Q i ∩ Q s )))(∃q ∈ Q i \ Q s )
(∃σ ∈ L(G, q ))P (σ ) = σ o (by P (L(G)) = P (L(G,

Q i ∩ Q s )) ∪ P (L(G, Q i \ Q s )) 3 ) ⇔ (∀q ∈ Q i ∩ Q s )(∀σ ∈ L(G, q))(∃q ∈ Q i \ Q s )
(∃σ ∈ L(G, q ))P (σ ) = P (σ) which is consistent with Def. 8 on I-opacity. 2

For any q aug obs ∈ Q aug obs , we define I(q aug obs ) = {q i ∈ Q i | ∃(q c , q i ) ∈ q aug obs } as the initial states of automaton G that it contains. The following theorem illustrates that the I-observer G aug obs can be used to verify the I-opacity for UWA G.

Theorem 4 Given a UWA G = (Q, E, α, µ), projection P : E * → E * o , and a set of secret states Q s ⊆ Q, G is initial-state opaque with respect to Q s and P iff (∀q aug obs ∈ Q aug obs )I(q aug obs ) Q s where Q aug obs is the set of states in I-observer G aug obs of G. Proof: It follows from Lemmas 3 and 5.

2

Example 9 Consider the I-observer G aug obs in Fig. 3 of UWA G in Fig. 1. It can be checked that G is always nonopaque with respect to P and an arbitrary secret Q s so long as Q s ∩ Q i = ∅. In fact, among the states of G aug obs , we have I({3, 1}) = {1} and I({3, 4}) = I({7, 4}) = {4}. Therefore, for a secret Q s containing initial state 1 and/or 4, there necessarily exist one state q aug obs ∈ Q aug obs such that I(q aug obs ) ⊆ Q s . 2

Conclusion

In this paper, an I-observer is constructed to verify the strong I-detectability, weak I-detectability and I-opacity of UWAs. Besides, a diagnosability-based approach, with a lower complexity, is proposed for checking strong I-detectability. Our future work is to explore new techniques that can handle the I-detectability and Iopacity problem for more general classes of WAs.

Definition 4

 4 d) 5,7 = 2, and α =(e, ε, ε, e, ε, ε, ε). All the non-listed coefficients in µ(u), µ(b), µ(c) and µ(d) are equal to ε, meaning that they do not model possible transitions in the automaton. G is not deterministic since it has two initial states.2 Given an arbitrary path π

Fig. 1 .

 1 Fig. 1. An unambiguous weighted automaton G.

Fig. 2 .

 2 Fig. 2. Augmented automaton G aug of G in Fig. 1.

Fig. 3 .

 3 Fig. 3. I-observer G aug obs of G in Fig. 1.

Fig. 4 .

 4 Fig. 4. Modified automaton G of G in Fig. 1.

Example 7 Fig. 5 .

 75 Fig. 5. Diagram of Go of G with respect to fault f1.

(

Fig. 6 .

 6 Fig. 6. Diagram of G d with respect to fault f1.

  aug obs , i.e., O(|Q aug obs | + |Q aug obs | × |E aug obs | × |Q aug obs |). By Algorithm 1, |Q aug obs | is bounded by 2 |Q| 2 -1. Hence, the complexity of verifying strong and weak Idetectability using I-observer is O(2 2|Q| 2 ×|E aug obs |), where |E aug obs | is usually exponential in |Q|.

	2
	4.2 Diagnosability-based Approach for Strong I-
	Detectability

  Algorithm 2 should be applied |Q i | times for testing the diagnosability of G with respect to |Q i | different fault types. Therefore, by combining Lemma 4, the complexity of verifying strong I-detectability based on diagnosability for UWA G is O(|Q| 4 × |E d | × |Q i |), where |E d | is exponential in |Q| in general. Note that under the assumption that all the unobservable events of UWA G are represented by a unique symbol, |E d | is bounded by |Q| 3 × |E o |, and the above exponential complexity is reduced to polynomial. 2 Example 8 Consider G d presented in Fig. 6. According to Theorem 3, we know that the corresponding UWA G in Fig. 4 is diagnosable with respect to fault f 1 . Similarly, it can be checked that G is diagnosable with respect to fault f 2 by constructing another G d on f 2 (here G d is omitted). Therefore, according to Lemma 4, the original system G in Fig. 1 is strongly I-detectable, which is consistent with the result in Example 5. 2 5 Verification of Initial-State Opacity

More precisely, (⇒) holds true simply because P (L(G, Qi∩ Qs)) ⊆ P (L(G)). In order to prove the correctness of (⇐), we have to prove that for all σo ∈ P (L(G, Qi \ Qs)), there exist q ∈ Qi \ Qs and σ ∈ L(G, q ) such that P (σ ) = σo, which is self-evident.

This work is partially supported by the National Key R&D Program of China under Grant 2018YFB1700104, and National Natural Science Foundation of China under Grant 61873442.