Qiang Song
email: qiang.song@imt-atlantique.fr

Yoran Eli Pigeon

Kevin Heggarty

Faceted Fresnel DOEs creating the perception of a floating 3D virtual object under divergent illumination

Keywords: Micro optics, Diffractive optics, Iterative algorithm, Binary Optics, Photo-Lithography

An approach for the optimization and fabrication of a phase-only faceted Fresnel type diffractive optical element (FDOE) creating 3D virtual object is proposed. The FDOE is a transmissive Fresnel type DOE array, which produces the perception of a customized floating 3D virtual object behind the FDOE when illuminated with a divergent monochromatic Light Emitter Diode (LED) source. Each DOE unit of the FDOE is optimized by a modified iterative Fourier transform algorithm (M-IFTA). Every unit of the FDOE locally deflects the incident light to the same position to form a designated view in the target plane. The FDOE is fabricated using our home-built parallel writing photo-lithography machine. Numerical simulations and optical experiments are performed to verify the proposed design method. This work may find important applications in the advanced design of optical security hologram and anti-counterfeiting component.

Introduction

Security holograms are one of the most important optical security technologies, which are widely used in anticounterfeiting protection of bank notes , passport and ID cards [START_REF] Van Renesse | Optical Document Security[END_REF][START_REF] Yue | Nanometric holograms based on a topological insulator material[END_REF][START_REF] Ferri | Visualization of 3D information with digital holography using laser printers[END_REF][START_REF] Lee | Optical configuration and color-representation range of a variable-pitch dot matrix holographic printer[END_REF][START_REF] Kevin | Holographic colour prints for enhanced optical security by combined phase and amplitude control[END_REF][START_REF] Andrulevicius | Hologram Origination Combining Rainbow and Dot-matrix Holograms[END_REF][START_REF] Vasconcellos | Printable Surface Holograms via Laser Ablation[END_REF][START_REF] Firsov | Fabrication of digital rainbow holograms and 3-D imaging using SEM based e-beam lithography[END_REF][START_REF] Skeren | Synthetic diffractive elements for security applications realized on an enhanced integral dot-matrix system[END_REF][START_REF] Lee | Pixelgram: an application of electron-beam lithography for the security printing industry[END_REF][START_REF] Huang | Large-format grating image hologram based on e-beam lithography[END_REF][START_REF] Girnyk | Stereographic and animated rainbow diffractive images in optical security[END_REF]. There are two main approaches for producing security holograms: one is based on the original holography technique, using light-sensitive materials to record the optical interference fringes with an optical interferometer system. This requirement makes it difficult to produce nonphysical or virtual objects, and to address cost-effective mass production. The other main approach uses surface relief diffractive optical elements (DOEs), whose advantage compared to optically recorded holograms is that they do not require the interferometer system. This enables mathematical models of diffraction processes to be used with computer programming, which is more suitable for making optical security holograms. Security holograms which can produce new visual perception are still desperately needed for modern anti-counterfeiting protection. What fundamentally new visual perception effect can we expect?

A DOE generating perceived floating 3D virtual object under divergent illumination might be a good candidate, because 3D floating virtual object can provide observers a perception close to real life. In addition, a LED source is universally available in daily life, and safety for the human eye. Several related research literatures on this aspect have appeared in recent years. For instance, Anton Goncharsky and his co-authors proposed a subwavelength grating array method using electron-beam lithography (EBL) technology that can generate visual 2D and 3D images on the surface of the DOE plane [START_REF] Goncharsky | Diffractive optical element with asymmetric microrelief for creating visual security features[END_REF][START_REF] Goncharsky | Diffractive optical element for creating visual 3D images[END_REF][START_REF] Goncharsky | High-resolution full-parallax computergenerated holographic stereogram created by e-beam technology[END_REF]. This can give an observer an impressive 3D vision effect, but, not a floating visual effect. In addition, the diffraction efficiency of the gratings is generally low, and the higher orders can produce crosstalk. Multilevel gratings or blazed gratings are an effective way to improve the diffraction efficiency, but the fabrication process is very complicated with EBL. Later, the same team used a similar method to realize 360 degree 3D image with grating arrays fabricated on a cylindrical surface [START_REF] Goncharsky | Cylindrical computer-generated hologram for displaying 3D images[END_REF]. Furthermore, a recent study presented a frequency mosaic method using maskless lithography which can produce 3D visual perception [START_REF] Yang | Computer generated full-parallax synthetic hologram based on frequency mosaic[END_REF]. Nevertherless, the diffraction efficiency is low and faraway LED illumination is needed to reproduce the virtual object on the surface of the hologram.

In this work we propose a new way to produce a floating 3D virtual image with divergent LED source, which we name faceted pure phase Fresnel DOEs. To the best of our knowledge, this is the first time it is being reported to the optics community. LED sources are used because they are less dangerous for visual perception than lasers, and are also more convenient, as no beam-expanding optics are required, since it is a divergent light source. On the surface of the FDOE, each DOE facet diffracts incident light to generate different virtual viewing angle images toward to the observer position; it therefore belongs to the off-axis type DOE. To obtain high image quality and diffraction efficiency , a modified iterative Fourier transform algorithm (IFTA) is used to optimize each DOE facet.

After design, the FDOE is fabricated with our home-built gray scale direct parallel writing photolithography machine.

Gray-scale maskless photolithography is an effective technique for fabricating high quality complicated micro relief structures, which can help increase DOE's diffraction efficiency. When the FDOE is illuminated by a divergent LED source, a floating 3D virtual image is formed behind the FDOE; an observer can perceive different views of the virtual 3D object hanging in space. The configuration is described in detail in Section 2. Section 3 is devoted to illustrating the design procedure. Fabrication and experiment results are presented in Section 4, and results are discussed in Section 5. Finally, we conclude our work in Section 6.

Principle and configuration

Fig.1. The configuration of the LED source and observer

The essential idea of the method proposed in this work is sketched in Fig. 1. Fig. 1 describes the layout of observations which comprises the LED source and FDOE. The LED source is behind the FDOE. To produce the 3D desired visual effect, we propose an approach based on faceted off-axis Fresnel DOEs. In such an approach the overall diffractive structure is separated into a matrix of individual Fresnel DOEs. Each facet diffracts to generate a view of the same basic 3D object, but seen from a slightly different angle. By carefully optimizing the angles for each facet, when an observer moves their eye from side to side at the FDOE plane, they observe the different views, generating the slightly different angular views of the 3D object. The overall effect is that of observing a real 3D object. To make sure the observer has the impression that each view comes from the same object in the same position, the different DOE facets must generate their reconstructed images with the correct off-axis offset. The lateral displacement of each reconstructed image should be calculated carefully. The basic calculation principle is shown below in Fig. 2.

The global coordinates of the FDOE plane are denoted by (x0,y0), and the local coordinates of the m-th DOE are denoted by (xm, ym) . The distance between the FDOE plane and the virtual image plane is denoted by -z, where z is a positive value, the distance is negative. The complex amplitude distribution of (m, n)-th DOE on FDOE plane can be described as followings: 𝐻 𝑚,𝑛 (𝑥 𝑚 , 𝑦 𝑛 ) = 𝐴 𝑠𝑜𝑢𝑟𝑐𝑒 exp [𝑗𝜑 𝑚,𝑛 ]

(1)

Fig.2. Principle of proposed FDOE calculation process

where 𝐴 𝑠𝑜𝑢𝑟𝑐𝑒 and 𝜑 𝑚 are the amplitude of the incident light and, the m-th DOE's phase distribution, respectively.

According to Fresnel diffraction theory, the diffracted field distribution of the (0,0)-th DOE in the virtual image plane can be calculated by

In eq.( 2), k=2π/λ is the wave number in free space, λ is the design wavelength, and x0, y0, x, y, are the coordinates axis in Fig. 2, respectively. In order to generate the reconstructed virtual image of the (m, n)-th DOE toward to the same position as that of the (0,0)-th DOE , eq.( 2) should be modified. See in Fig. 2, xm and yn can be expressed as

𝑥 𝑚 = 𝑥 0 + 𝑚𝑃 𝐷𝑂𝐸 (3)
𝑦 𝑛 = 𝑦 0 + 𝑛𝑃 𝐷𝑂𝐸 [START_REF] Lee | Optical configuration and color-representation range of a variable-pitch dot matrix holographic printer[END_REF] In eq.(3) and eq.( 4), 𝑃 𝐷𝑂𝐸 represents the facet size of the FDOE, m and n are the displacement along the x0 axis and the y0 axis, respectively. According to Fresnel diffraction theory, the complex amplitude distribution Um,n(x, y) in the virtual image plane contributed by the wavefront from the (m, n)-th DOE can be calculated as

𝑈 𝑚,𝑛 (x, y) = exp(-𝑗𝑘𝑧) -𝑗𝜆𝑧 ∬ 𝐻 𝑚,𝑛 (𝑥 𝑚 , 𝑦 𝑛 ) 𝑒 𝑗𝑘 -2𝑧
[(𝑥-𝑥 𝑚 ) 2 +(𝑦-𝑦 𝑛 ) 2 𝑑𝑥 𝑚 𝑑𝑦 𝑛 [START_REF] Kevin | Holographic colour prints for enhanced optical security by combined phase and amplitude control[END_REF] Substituting eq.(3) and eq.( 4) into eq.( 5),

𝑈 𝑚,𝑛 (x, y) = 𝑒 -𝑗𝑘𝑧 -𝑗𝜆𝑧 e [ jk -2z (x 2 +y 2 )] 𝑒 𝑗𝑘(𝑚𝑃 𝐷𝑂𝐸 𝑥+𝑛𝑃 𝐷𝑂𝐸 𝑦) 𝑧 × ∬ 𝐻 𝑚,𝑛 (𝑥 𝑚 , 𝑦 𝑛 )𝑒 𝑗𝑘 -2𝑧 [(𝑥 0 +𝑚𝑃 𝐷𝑂𝐸 ) 2 +(𝑦 0 +𝑛𝑃 𝐷𝑂𝐸 ) 2 ] 𝑒 𝑗𝑘 -2𝑧 (𝑥𝑥 0 +𝑦𝑦 0 ) 𝑑𝑥 0 𝑑𝑦 0 (6)
Where, 𝑓 𝑥 ,𝑓 𝑦 are defined as

𝑓 𝑥 = 𝑥 𝜆𝑧 𝑓 𝑦 = 𝑦 𝜆𝑧 (7)
With the 2D Fourier transform, eq.( 6) can be simplified as

𝑈 𝑚,𝑛 (x, y) = 𝑒 -𝑗𝑘𝑧 -𝑗𝜆𝑧 𝑒 jk -2z (x 2 +y 2 ) 𝑒 𝑗𝑘(𝑚𝑃 𝐷𝑂𝐸 𝑥+𝑛𝑃 𝐷𝑂𝐸 𝑦) 𝑧 • ℱ 2𝐷 [𝐻 𝑚,𝑛 𝑒 𝑗𝑘 -2𝑧 [(𝑥 0 +𝑚𝑃 𝐷𝑂𝐸 ) 2 +(𝑦 0 +𝑛𝑃 𝐷𝑂𝐸 ) 2 ] ] (8)
According to eq.( 8), the diffracted field of the off-axis Fresnel DOE can still be expressed in a 2D Fourier transform format. Therefore, it can be computed quickly with the fast Fourier transform method. Sampling constraints for the Fresnel diffraction calculation should also be considered. According to the Nyquist sampling rule, the absolute value of the propagation distance z should satisfy the condition as below [START_REF] Nguyen | Computational efficient scalar non-paraxial modeling of optical wave propagation in the far-field[END_REF]:

z > 𝑁∆𝑝 2 λ (9)
Here, in eq.( 9), N is the number of sampling points and ∆𝑝 is the DOE pitch. Generally, it is a phase retrieval problem to get the DOE phase. In this work, a modified iterative Fourier transform algorithm is implemented to design each off-axis Fresnel type DOE, it belongs to the IFTA family [START_REF] Hao | Wide angle 2D bean splitter design based on vector diffraction theory[END_REF][START_REF] Wang | A hybrid Gerchberg-Saxton-like algorithm for DOE and CGH calculation[END_REF][START_REF] Wyrowski | Diffractive optical elements: iterative calculation of quantized, blazed phase structures[END_REF][START_REF] Nguyen | Iterative scalar nonparaxial algorithm for the design of Fourier phase elements[END_REF]. The proposed method optimizes the image quality and the fabrication performance simultaneously incorporate with a soft quantization iterative algorithm. The fabricated phase levels are quantized, so the surface relief of the DOE should be divided into several levels. In theory, the higher the number of phase levels, the better the reconstructed image quality and diffraction efficiency. But, the efficiency increase for greater numbers of phase level is often less than the decrease from fabrication errors. Thus, eight-levels is a good practical compromise. The details of the design procedures are provided in Section 3.

Design procedure and numerical simulations 3.1 Design method

The flowchart of the proposed design approach is shown in Fig. 3. The stack of different viewing angle target images is obtained from a 3D model. Each target image is input to the DOE optimization algorithm. This iterative algorithm can be divided into three stages, the so called three-part IFTA. The first stage is the standard Gerchberg-Saxton (GS) algorithm, it is implemented to get the rough phase solution, the GS operation process can be considered as the transformation of signals between the input DOE plane and diffracted image plane.

Fig.3. Flowchart of the proposed iterative optimization process of FDOE

The GS iteration process can be described as four basic steps: (1) the process starts from the diffracted image plane.

An initial random phase is generated in the range 0-2π, and then multiplied by an amplitude matrix with the current viewing angle target image to form the initial complex amplitude. [START_REF] Yue | Nanometric holograms based on a topological insulator material[END_REF] This field is propagated to the FDOE plane using inverse Fresnel diffraction field propagation. [START_REF] Ferri | Visualization of 3D information with digital holography using laser printers[END_REF] Note that in what follows, we assume that the DOE facets are illuminated by plane waves. Illuminating such a DOE with spherical wave will simply shift the position of the observed image along the optical axis; see discussion part for more details. The amplitude of the wave-front 𝑂(𝑥 0 , 𝑦 0 ) is replaced by the intensity distribution of the source amplitude, while the wavefront phase is retained, to form the updated complex amplitude 𝐻(𝑥 0 , 𝑦 0 ) at the DOE plane. ( 4) Propagating the result back to the diffracted image plane, the phase obtained from the Fresnel diffraction propagation is combined with the corresponding view image amplitude. The process then repeats from step 1. Based on the above steps, the phase is optimized for 20 iterations, then the algorithm goes on to the second stage, where the iteration process is similar to the GS, but the constraint factor strategy on the diffracted image plane is different to the conventional GS algorithm, which is given as

|𝑈 𝑘 | = { 2𝛽 𝑘 |𝐹 𝑡𝑎𝑟𝑔𝑒𝑡 | -|𝐹 𝑘 |,𝐹 𝑘 ∈ 𝑆 |𝐹 𝑘 |,𝑒𝑙𝑠𝑒 (10)
Where, 𝐹 𝐾 is the reconstruction result of the k-th iteration, 𝐹 𝑡𝑎𝑟𝑔𝑒𝑡 is the target amplitude, 𝛽 𝑘 is a weighting factor , and S is the signal area, see Fig. 4. For a more detailed depiction of the modified amplitude constraint, we refer the interested readers to our previous work [START_REF] Wang | A hybrid Gerchberg-Saxton-like algorithm for DOE and CGH calculation[END_REF]. After that, the optimization process goes into the third stage; this stage is called the iterative soft quantization method. The soft quantization is a stepwise iterative quantization method. In this part, the iterative process is similar to that of the second stage, but another loop is introduced into it. This process of iteration can be divided into 10 cycles, each cycle includes 30 iterations. The soft quantization operation can be demonstrated as [START_REF] Wyrowski | Diffractive optical elements: iterative calculation of quantized, blazed phase structures[END_REF] :

First we define 𝜑 1 as

𝜑 1 = 𝜑 𝑝ℎ𝑎𝑠𝑒 𝑘 2𝜋 𝐿𝑒𝑣𝑒𝑙𝑠 (11) 𝜑′ 𝑝ℎ𝑎𝑠𝑒 𝑘 = 𝜑 1 , if - 𝜀 𝑃 2 ≤ 𝜑 1 -𝑟𝑜𝑢𝑛𝑑(𝜑 1 ) < 𝜀 𝑃 2 (1 < 𝑝 < 10)(12)
In eq.( 10), 𝜑 𝑝ℎ𝑎𝑠𝑒 𝑘 is the un-quantized phase of the k-th iteration, Levels is number of phase levels, 𝜑′ 𝑝ℎ𝑎𝑠𝑒 𝑘 represents the quantized phase in eq.( 12). The quantity 𝜀 𝑃 increases with the index p, p is the current loop, the maximum loop is 10 in this work. p is an empirical parameter. Here, the collection of the values in this work is 𝜀 1 = 0.15, 𝜀 2 = 0.3, 𝜀 3 = 0.45, 𝜀 4 = 0.6, 𝜀 5 = 0.7, 𝜀 6 = 0.8, 𝜀 7 = 0.85, 𝜀 8 = 0.9, 𝜀 9 = 0.95, 𝜀 10 = 1,

The reason why we choose 30 iterations in this process are discussed in next sub-section where the convergence behavior is described. So, in this way, the phase can be quantized into eight levels with high a quality reconstructed image, and adapted to the fabrication constraints. The iterative calculation stops when the error between the intensity distribution on the reconstructed plane and target image is small enough, or the maximum number of iterations is reached. The root mean square error (RMSE) can be defined as

RMSE = ∬(𝑈 𝑡𝑎𝑟𝑔𝑒𝑡 -𝛾|𝑈 𝐾 |) 2 ∬ 𝑈 𝑡𝑎𝑟𝑔𝑒𝑡 2 ( 14
)
where 𝛾is a scale factor, which can be written as

𝛾 = ( ∬ 𝑈 𝑡𝑎𝑟𝑔𝑒𝑡 2 ∬|𝑈 𝐾 | 2 )( 15
)
When calculations of all the DOE facets are finished, the set of DOE phases are assembled to get the final FDOE.

Numerical simulation

For the numerical simulations, a 3D "cube" object was used to generate the different viewing angle images. In this study, we use horizontal views. To examine the convergence behavior of the proposed optimization algorithm, one view of these image is selected as the target image, which is shown in Fig. 5(a), it is the center view of the 3D object.

The RMSE value is plotted against the iteration numbers in Fig. 5(b). In this numerical simulation work, the wavelength is 525nm, the DOE pixel pitch is 0.75μm, the number of sampling points is 3600 in each DOE facet, and the distance between the FDOE plane and the virtual image is set as -120mm. Acoording to eq.( 14 [START_REF] Goncharsky | Cylindrical computer-generated hologram for displaying 3D images[END_REF] Where I is the intensity distribution of the target plane, M and N are the number of sampling points in the X and the Y direction , respectively.

For our application, the aim is to generate horizontal parallax 3D views. 11 off-axis Fresnel type DOE elements are used. The size of the FDOE is 29.7mm×2.7mm, it can be regarded as the system eye-box. However, this is too small to be convenient in the vertical direction, as an observer will no longer see the virtual object if their eye moves beyond this small area. The pupil size of human eye is between 6mm~8mm, so the vertical direction should be expanded to satisfy this size. To address this problem, we propose a method to increase the vertical eye-box size.

The assembly is shown in Fig. 7.

Fabrication and experiment results

With the algorithm mentioned above, a FDOE is designed with 11×3 facets. The smallest pixel size of our photolithography machine is 0.75μm, but here we use 1.5μm as the design pixel size. It helps minimize pixel deformation in the off-axis configuration which tends to decrease the DOE fringe size. Each DOE facet is1800×1800 pixels with pixel size 1.5μm, the total FDOE size is 29.7mm×8.1mm.The FDOE was fabricated with our home-built parallel direct-writing photolithography machine, an updated version of that demonstrated in [START_REF] Kessels | Versatile stepper based maskless microlithography using a liquid crystal display for direct write of binary and multilevel microstructures[END_REF][START_REF] Kessels | Optical proximity correction for a versatile LCD based direct write maskless photoplotter[END_REF][START_REF] Kessels | Direct write of optical diffractive elements and planar waveguides with a digital micromirror device based UV photoplotter[END_REF]. In Fig. 8, the light source is LED based with a center wavelength of 435nm. The divergent light is first homogenized using group of lenses and a diffuser. Then, the collimated and uniform light field is modulated with the designed phase pattern loaded onto a 1920×1080 pixel Spatial Light Modulator(SLM) (Epson L3D135-55g00) and de-magnified onto the surface of the photoresist on the precision moving platform through the objective lens group. The magnification ratio of the objective lens group is 5. The projected area of SLM is small, about 4mm×3.2mm, so the FDOE is made by stitching different SLM images by moving the stage in the x and y directions with an accuracy of about 250nm. The fabrication process is now described briefly. The glass substrate is cleaned with ultrasonic water cleaning, and then heated for 15 minutes to dry the glass on a hot plate. The next step is the spin coating of a photoresist onto the substrate. In this work, the photoresist is the S1813 from Micro Resist Technology GmbH. The spin coating speed is chosen to give uniform photoresist layer. Etch depth is chosen to give the maximum phase difference of each cell is 2π. The etch depth is dependent on the wavelength λ of the incident light and the refractive index of the resist.

According to the thin element approximation (TEA), the groove depth of the surface relief structure can be written where number of phase levels, N=8, 𝑛 𝑚𝑎𝑡𝑒𝑟𝑖𝑎𝑙 is the refractive index of photoresist, its value is 1.66 at a wavelength of 532nm. According to the formula, the maximum depth d is about 691nm. The spin coating speed is 3200rpm, giving a total photoresist thickness about 1.3μm; this depth is enough for etching the DOE. The substrate is baked for 2 minutes at temperature of 120 degree to extract the solvent, and then cooled down to room temperature. Then, exposure is implemented using our home-built parallel direct-writing photo-plotter with a lookup table (LUT) to establish a linear relation between the addressed gray phase level and the etch depth. The exposed pattern is mapped from the SLM screen to photoresist plane. This step needs about 15 minutes. The exposure time and LUT can be controlled to ensure the optimum etching depth. The last step is to develop the exposured photoresist in a chemical solution for 2 minutes. A fragment of the FDOE and final fabricated FDOE are shown in Fig. 9 (a The schematic of experimental set up used to capture the output views is shown in Fig. 10. This optical system comprises an LED source, a FDOE and a camera. A green LED source behind the FDOE is used as the illumination source for reproducing the expected 3D virtual image. The camera is mounted on the support, which can be rotated to gain different viewing angles equivalent to those can observer sees moving his head to view the apparent 3D object. The distance between LED and FDOE is about 155mm, and position of the virtual image is about 70mm below the FDOE. The position of virtual image can be calculated according to the imaging formula, which is expressed in eq. ( 18).

d 𝑖𝑚𝑎𝑔𝑒 =

𝑓 𝑓𝑟𝑒𝑠𝑛𝑒𝑙 𝑑 𝑠 𝑑 𝑠 -𝑓 𝑓𝑟𝑒𝑠𝑛𝑒𝑙 [START_REF] Nguyen | Computational efficient scalar non-paraxial modeling of optical wave propagation in the far-field[END_REF] In eq.( 18), 𝑓 𝑓𝑟𝑒𝑠𝑛𝑒𝑙 is the designed focal length of the Fresnel type DOE whose value is -120mm, 𝑑 𝑠 is the distance between the source and the FDOE plane. With the LED light illumination, the light diffracted by the FDOE propagates to human eye and the 3D virtual image is perceived in the air between the DOE plane and the LED source.

Fig. 11 shows one perspective image of the 3D target "cube". When the hologram is optically replayed, a floating virtual image appears hanging in space behind the FDOE, a finger is near to the virtual image. An observer also can see a realistic 3D image by looking at the virtual object from different angle. Fig. 12

Discussion

This paper has presented the design process for a diffractive structure, based on an array of computer generated offaxis Fresnel phase holograms. When illuminated by a readily available LED source, the overall structure gives an observer the perception of a floating 3D object in space behind the hologram. The design model has been described and confirmed through both numerical simulation and experimental results. Although, this work focused on presenting the design methodology to generate perceived 3D virtual vision in the horizontal direction, the fullparallax vision can also be generated with the same approach. As illustrated in Fig. 7, we just let the Hm1, Hm2 and shift the image to the desired position. In Fig. 13 (a), when the LED source is on the axis of the DOE1, the virtual image is still formed on the axis because there is no blazed grating component in the DOE1. The position of the generated virtual image can be calculated with the geometrical imaging formula of eq. ( 18). The DOE2 and DOE3 contain the blazed grating component, so, their corresponding virtual image can be deflected to the same position with DOE1. When the LED source is moved with respect to the original position (see Fig. 13 (b)), the virtual image will move in the same direction, this process is similar to a lens imaging system, because the Fresnel lens component works as a lens. In this situation, the corresponding images generated by DOE2 and DOE3 just move in the same way.

But, one thing must be borne in mind, the LED source is a divergent source, and the intensity is cosine with the divergent angle. The bigger divergent angle, the weaker the intensity. So, the shifting of the LED source cannot be too far.

Another important factor is the spatial coherence of the LED source, which has been discussed in detail in reference [START_REF] Deng | Coherence properties of different light sources and their effect on the image sharpness and speckle of holographic displays[END_REF]. According to the results, the size of the LED cannot be too big; otherwise, the reconstructed image will be not sharp due to the low degree of spatial coherence. In this work, the LED size is 1mm×1mm, but as the distance from the LED source to FDOE plane is more than 150mm, the LED source can be still regarded as 'point source'.

Another essential DOE parameter to be evaluated is the diffraction efficiency: the image must be sufficiently bright to be visible to the observer. Unfortunately measuring the diffraction efficiency of this type of DOE is not straightforward because the image is virtual and the divergent illumination. To obtain an estimate of the diffraction efficiency we used the setup shown in Fig. 14. The condenser lens is used to convert the virtual image into a real image which is projected onto the photodetector. A diaphragm is placed between the LED source and the DOE so as to illuminate only the active DOE area. Without the DOE, the total collected light power on the detector gave an estimate of the incident light power. With the DOE present the light power in the desired image area gave an estimate of the light power in the signal. In this way we obtained a DOE diffraction efficiency estimate of 46%. The DOE zeroth order was estimated in a similar way at 20%. The remaining light losses due mainly to Fresnel reflection (~10%) and higher diffraction orders. With this preliminary un-optimized DOE, illuminated with divergent nonmonochromatic (bandwidth ~33nm) light such light losses are not unexpected. We are currently improving our fabrication facilities to decrease the DOE pixel size which will increase diffraction efficiency by suppressing the higher diffraction orders [START_REF] Schwarzer | Non-pixelated" design of computer-generated diffractive elements for increased diffraction efficiency[END_REF] and optimizing the setup and algorithm to make the DOE easier to fabricate. Based on values obtained with similar DOEs, we expect to be able to reduce zero order power to a few percent. However, even with the current diffraction efficiency of 46%, the observed image is already clearly visible in daylight with a readily available LED.

Conclusion

This paper has presented a computer generated horizontal moving parallax synthetic DOE based on facet Fresnel type DOE array. The design model for the FDOE has been demonstrated through both numerical simulation and experimental results. In this work, we limit viewing of the FDOE to a horizontal displacement of the observer's eye exclusively. A visually convincing 3D virtual object can be observed when the FDOE is illuminated with a divergent LED source. The proposed optimization method can also be implemented to calculate a full-parallax DOE, which has been described in Section 5. As the FDOE is a surface relief diffractive structure, the FDOE can be replicated

using roll to roll nanoimprint techniques for mass production. Therefore, the FDOE developed in this paper can find applications in anti-counterfeiting protecting of bank notes and ID cards. The approach may potentially also be used in practical holographic 3D display application.

Fig. 4 .

 4 Fig.4. The definition of signal area in the diffracted image plane.

Fig. 5 .

 5 Fig.5. Target image and RMSE curves +20(the second stage)+500(the third stage)). In Fig.5 (b), the RMSE is small enough when the total number of iteration is 340. To further compare the convergence behavior between GS and the proposed algorithm in this work, the RMSEs of two algorithm against versus the iteration number are plotted in Fig.5(c), respectively. The results suggest that the proposed algorithm can converge to a smaller RMSE than GS algorithm when the soft quantization process is implemented. The GPU is used to accelerate the calculation speed; the iteration optimization time is about 120 seconds on a computer equipped with Intel Core™ I7-8750H CPU@2.2GHz, GPU GTX1060 and the computation platform is MATLAB software. The final optimized DOE is shown in Fig.6 (a) , the inset is a cross section of the phase value distribution in the red circle. It is clearly indicating an 8 phase-level structure. The simulated reconstructed virtual image is shown in Fig.6 (b), the reconstructed image is similar to the target image in Fig.5 (a). But, the reconstructed image is a bit darker than the target image, because the calculated diffraction efficiency is about 82%, part of the energy being lost as noise. The diffraction efficiency is defined as 𝜂 = ∑ ∑ 𝐼(𝐼∈𝑆) 𝑁 𝑀 ∑ ∑ 𝐼 𝑁 𝑀

Fig. 6 .

 6 Fig.6. Optimal phase map and the reconstructed result

Fig. 7 .

 7 Fig.7. Assemble of the FDOE In Fig.7, the elements Hm0 , Hm1 and Hm3 generate the same virtual image Um at the same position. When people move their eye in vertical direction, they can still see the same virtual image at specific view. In this way, the eyebox can be expanded in vertical direction, and the final size of the FDOE is 29.7mm×8.1mm. Due to the use of divergent illumination and DOE fabrication limitations, the LED source is visible to the observer. If the FDOE are designed with an on axis target, this visibility of the LED source perturbs the desired image. To avoid this, an offaxis configuration is used.

Fig. 8 .

 8 Fig.8. Diagram of photo-lithography system

  ) and Fig.9 (b), respectively.

Fig. 9 .

 9 Fig.9. Fragment part of FDOE and real sample.(a) Fragment of FDOE;(b) Real sample on photoresist with glass substrate. The photo of Fig.9 (a) is captured by an optical microscope (Reichert-Jung Polyvar Met Microscope), with a 100x objective.

  Fig.10. The different views all correspond to a single 3D object located in the same spatial position. The field of view (FOV) is about 20 degree. The brightness of the reconstructed images vary slightly because the exposure time of camera is slightly different for the different image capture angles. In practice, the visual effect perceived by the observer is of a 3D virtual image with more uniform brightness.

Fig. 10 .

 10 Fig.10. Schematic diagram of capturing the 3D floating virtual object reconstruction. Inset: spectrum of the LED source, which is obtained by Ocean Optics HR4000CG-UV-NIR, its center wavelength is 525nm.

Fig. 11 .Fig. 12 .

 1112 Fig.11. One perspective of the 3D cube floating in the air

Fig. 13 .

 13 Fig.13. Schematic of imaging principle, (a) LED source is on the axis, (b) LED source is off axis .

Fig. 14 .

 14 Fig.14. Setup for diffraction parameters measurement

Acknowledgments

The authors acknowledge the support by the European Union's Horizon 2020 research and innovation program through PHENOmenon under grant agreement No.780278 and Brest Metropole.

Funding

European Union's Horizon 2020 (780278)


Qiang Song
email: qiang.song@imt-atlantique.fr

Yoran Eli Pigeon

Xavier Theillier

Kevin Heggarty

An optimization approach of computer generated hologram (CGH) for divergent light shaping

An optimization method of pure phase CGH for beam shaping with divergent incoherent light is proposed. The sample experiment fabricated with our home built photo-plotter shows a good agreement with the numerical simulation results.

Introduction

The researches of divergent light illumination are mainly concentrated in refractive or reflective optical element (ROE) for a long time [START_REF] Oliker | Optical design of freeform two-mirror beam-shaping systems[END_REF][START_REF] Song | Free-form illumination of a refractive surface using multiple-faceted refractors[END_REF]. Although these ROE based method can get an asymmetric pattern with a freeform lens or facet surface structure, the large scale size and thickness limit their applications.

Computer generated hologram (CGH) is a good choice to solve this problem up to now. However, there is rarely research work relevant with the CGH design illuminated by divergent incoherent light. Here in this paper, the strategy of divergent source model, complex degree of coherence and a modified magnitude constraint are introduced into the standard iterative Fourier transform algorithm(IFTA) process [START_REF] Wyrowski | Diffractive optical elements: iterative calculation of quantized, blazed phase structures[END_REF]. With this proposed method, a good performance for incoherent light shaping is obtained.

Methodology and Results

The classical IFTA is one of the most popular algorithms, which has been widely used in phase retrieval problem and calculation of CGHs. Numerically, the IFTA is a forward and inverse Fourier transformation process of a complex function between the real(object) space and Fourier(image) space, in which the constraints are applied to the computed wave-front by a prior knowledge. Generally the standard IFTA is valid only in paraxial region and often used for plane wave shaping. So the model should be modified when the divergent light source is used. The proposed strategy is to use the scalar non-paraxial diffraction formula as the propagating function between two domains [START_REF] Nguyen | Iterative scalar nonparaxial algorithm for the design of Fourier phase elements[END_REF], and introduce the divergent light source model in the object space and the modified constraint in the Fourier space, as shown in Fig. 1. The constraint in the kth iteration is implemented by modifying |Ek| in the way as

ℛ[𝐸 𝑘 ] = { ℇ 𝑘 [𝜇 𝑘 (𝜔 1 𝐸 𝑡𝑎𝑟𝑔𝑒𝑡 -𝜔 2 𝐸 𝑘 )] Ω ∈ 𝑠𝑖𝑔𝑛𝑎𝑙 𝑤𝑖𝑛𝑑𝑜𝑤 𝐸 𝑘 Ω ∉ 𝑠𝑖𝑔𝑛𝑎𝑙 𝑤𝑖𝑛𝑑𝑜𝑤 (1)
ℇ 𝑘 = ∑ 𝐸 𝑡𝑎𝑟𝑔𝑒𝑡 𝛺 𝑠 ∑ 𝐸 𝑘 𝛺 𝑠 (2) ℛ[𝜇 𝑘 ] = 𝐸 𝑡𝑎𝑟𝑔𝑒𝑡 𝐸 𝑘 𝜇 𝑘 (3)
In eq.( 1), the subscript k indicates the number of iteration, 𝐸 𝑡𝑎𝑟𝑔𝑒𝑡 is the target amplitude in Fourier domain, 𝜔 1 and 𝜔 2 are empirical parameters. A binary image ′π′ is used to demonstrate the proposed algorithm . In this design, the distance between LED source and CGH is 50mm, and the effective area of the CGH is 12mm×12mm with pixel size 0.75µm , the collective divergent angle is 14°.

The fabricated phase levels are quantized, so the surface relief of the CGH should be divided into several levels. In theory, the higher the number of phase levels, the better the reconstructed image quality. However, the fabrication performance must also be considered. Empirically, eight levels is a good practical compromise. The optimized CGH was fabricated with our home-built photo-plotter [START_REF] Kessels | Direct write of optical diffractive elements and planar waveguides with a digital micromirror device based UV photoplotter[END_REF],

which has 16000×16000 pixels . When the fabricated CGH is illuminated by using a divergent LED (wavelength 525nm), see Fig.

Conclusion

In this work, the modified IFTA method is proposed by introducing the strategy of divergent light source model and scaler non-paraxial diffraction model. The optical experiment result shows high relevance to the numerical simulation result, which demonstrates that the proposed method has good performance in divergent light shaping. The proposed design method can also find applications in diffractive optical elements design for anti-counterfeiting protecting, as well as computational imaging.

Fig. 1 .

 1 Fig.1. Flow chart of the proposed optimization method.

  Fig.2. (a) Simulation result (b) Experiment setup (c) Intensity distribution on the screen as a result of illuminating the CGH with a green LED source.


Derdei M Bichara

Aboudramane Guiro

Abderrahman Iggidr

Diene Ngom

State and Parameter Estimation for a Class of Schistosomiasis Models

Keywords: Mathematics Subject Classification: 92D25, 92D30, 93B05, 93B07, 93C40 Schistosomiasis, MacDonald's model, observer, state estimation, parameter estimation, control strategies. 1

We develop a general framework to estimate the proportion of infected snails and snail-human transmission parameter of a class of models that describes the evolution of schistosomiasis. To do so, we consider simultaneously the dynamics of schistosomiasis, captured by the homogeneous version of the classical MacDonald's model, and the measurable output: the number of female schistosomes per single host. The proposed method consists of designing an auxiliary dynamical system, called observer, whose solutions converge exponentially to those of the system capturing the schistosomiasis model. Moreover, we derive an estimation of the snail-human transmission rate, an unknown but key parameter in the dynamics of schistosomiasis. These estimations are central in two of the strategies of controlling schistosomiasis, namely the use of molluscicides and mass drug administration. To further investigate control strategies on a larger scale, we consider a heterogeneous model which consists of an arbitrary number of human groups or patches and an arbitrary number of fresh-water sources, natural habitats of snails. Provided that the data of infected humans' worm burden in each patch or group is available, we provide a method of estimating the proportion of infected snails in each snail natural habitat, thereby providing a map on where to implement control strategy to mitigate or eliminate Schistosomiasis.

Introduction

Schistosomiasis, also known as bilharzia, is a major parasitic disease that causes over 200 million in morbidity, and over 700 million people are at risk across the world [START_REF]Neglected tropical diseases, Centers for Disease Control and Prevention[END_REF]54]. The disease is prevalent mostly in tropical and subtropical areas and falls under the umbrella of neglected tropical diseases and is topped only by malaria in terms of incidence and economic burden in endemic regions. Consequently, understanding the mechanisms that drive the transmission Schistosomiasis and providing tools to control it should be a worldwide priority.

Schistosomiasis is the end-product of a complex biological cycle that involves two hosts -humans and snails; two free-living transmission pathogen, namely, cercariae and miracidia; and a variety of environments. Indeed, the overall transmission goes as follows: Schistosoma eggs are excreted from humans into the environment via urine or faeces (depending on the species of schistosomes). The eggs hatch into miracidia after a contact with fresh water. These miracidia infect water snails by penetrating the latter through their soft tissue, which results in amplification of the numbers of parasites by asexual reproduction. The parasite leaves the snail intermediate host as cercariae which can penetrate the skin of humans in contact with water, thus completing the life cycle. Different intervention strategies have been implemented in an attempt to eradicate schistosomiasis or alleviate its morbidity in many endemic countries [START_REF] Inobaya | Prevention and control of schistosomiasis: a current perspective[END_REF]. Indeed, mass drug administration campaign was widely used with varying degrees of success. However, this method has its own set of challenges due factors such as the development of drug resistance [START_REF] Inobaya | Prevention and control of schistosomiasis: a current perspective[END_REF], the difficulties of sustaining coverage of a critical mass of the target population, and the scarcity of resources to undertake nationwide programs [START_REF] Gray | Schistosomiasis elimination: lessons from the past guide the future[END_REF][START_REF] Utzinger | Sustainable schistosomiasis control-the way forward[END_REF]. The second control strategy consists of controlling the intermediate host snails using chemical molluscicides. However, there are two main challenges to this control strategy: assessing the number of infected snails in a given fresh-water source, and assessing which snail habitats to target. The goal of this paper is to address these two issues.

The use of mathematical models in understanding the dynamics of the interactions between host and vectors goes back to Ross' 1911 paper [START_REF] Ross | The prevention of malaria[END_REF]. Indeed, in [START_REF] Ross | The prevention of malaria[END_REF], Ross provides a model capturing the overall dynamics of the malaria, a disease caused by the parasite Plasmodium falciparum and transmitted between human hosts and female Anopheles mosquitoes when the latter bites the former. As mentioned above, schistosomiasis has a somehow similar but more complex transmission cycle and its mathematical model was first derived by MacDonald [START_REF] Macdonald | The dynamics of helminth infections, with special reference to schistosomes[END_REF] and later developed in different directions (see [START_REF] Anderson | Helminth infections of humans: mathematical models, population dynamics, and control[END_REF][START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF][START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF][START_REF] Woolhouse | On the application of mathematical models of schistosome transmission dynamics. i. natural transmission[END_REF][START_REF]On the application of mathematical models of schistosome transmission dynamics. ii. control[END_REF] and the references therein). However, estimating the state variables and/or parameters of the model is a perennial problem in epidemiology in general and in understanding the dynamics of schistosomiasis in particular. For instance, some parameters in schistosomiasis models, such as the life span of snails, could be obtained experimentally (2-3 years) [START_REF] Feng | Estimation of parameters governing the transmission dynamics of schistosomes[END_REF][START_REF] Van Cleave | Studies on the life cycle of the snail, viviparus contectoides[END_REF][START_REF] Woolhouse | On the application of mathematical models of schistosome transmission dynamics. i. natural transmission[END_REF] or the life span of the human hosts. Some other parameters, however, are unknown and very difficult to estimate. This the case for the snail-human transmission rate. Indeed, this parameter represents a simplified com-pound value of a complex biological cascade of events leading to the infection [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF]. And so, estimating the value of this snail-human transmission term is crucial in a control strategy that consists of blocking the transmission of schistosomiasis. Furthermore, a facet of this strategy is the chemical control strategy, in which chemical molluscicides are used to target infected snails in a given environment. This is also a difficult task as it requires one to distinguish between uninfected and infected snails in natural environments. Although some ad hoc estimations of the snail-human transmission rate exist for particular cases [START_REF] Feng | Estimation of parameters governing the transmission dynamics of schistosomes[END_REF][START_REF] Woolhouse | On estimating the basic reproduction number for schistosoma haematobium[END_REF], a comprehensive and general method of estimating the snail-human transmission is lacking. Moreover, to the best of our knowledge, there is no method in the literature that estimates the infected snail in water source.

In this paper, we use elements of control theory to estimate the proportion of infected snails in a given freshwater source and the snail-human transmission rate. To this end, we consider the dynamics of schistosomiasis along with the measurable output. For schistosomiasis' dynamical model, we consider the classical MacDonald's model [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF][START_REF] Macdonald | The dynamics of helminth infections, with special reference to schistosomes[END_REF] and we assume that the measurable output is the state variable that represents the number of female schistosomes. Indeed, for schistosomiasis, it is widely reported that it is possible to measure the number of parasites per hosts [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF][START_REF] Woolhouse | On the application of mathematical models of schistosome transmission dynamics. i. natural transmission[END_REF][START_REF]On the application of mathematical models of schistosome transmission dynamics. ii. control[END_REF], using faeces or urine.

Therefore, the first focal point of this paper is to develop a general framework to estimate the proportion of infected snails and the snail-human transmission parameter using the available data and the temporal evolution of schistosomiasis. More precisely, the method consists of constructing an observer. An observer is essentially an auxiliary dynamical system that uses the available data on some state variables and/or parameters and delivers dynamical estimates of the unmeasured state variables. The solutions of this auxiliary dynamical system must converge (as fast as possible) towards the solutions of the original system. The notion of observers originates in the theory of systems [START_REF] Luenberger | An introduction to observers[END_REF]. The use of this powerful tool of control theory is recent in epidemiology [START_REF] Abdelhedi | Extended observer to estimate the spreading of contagious disease[END_REF][START_REF] Abdelhedi | Observer design for a class of nonlinear piecewise systems. Application to an epidemic model with treatment[END_REF][START_REF] Alvarez-Ramirez | Feedback control of the chemotherapy of HIV[END_REF][START_REF] Aronna | Interval observer for uncertain time-varying SIR-SI model of vector-borne disease[END_REF][START_REF] Bichara | Étude de modèles épidémiologiques: Stabilité, observation et estimation de paramètres[END_REF][START_REF] Bichara | On the estimation of sequestered infected erythrocytes in Plasmodium falciparum malaria patients[END_REF][START_REF] Bliman | Interval observers for SIR epidemic models subject to uncertain seasonality[END_REF][START_REF] Bliman | A class of nonlinear adaptive observers for sir epidemic model[END_REF][START_REF] Degue | Interval estimation of sequestered infected erythrocytes in malaria patients[END_REF][START_REF] Diaby | Observer design for a schistosomiasis model[END_REF][START_REF] Iggidr | State estimators for some epidemiological systems[END_REF][START_REF] Tami | Parameters and states estimation for dengue epidemic model[END_REF][START_REF] Tendeng | Etude de modèles de transmission de la Schistosomiase: Analyse mathématique, reconstruction des variables d'état et estimation des paramètres[END_REF][START_REF] Velasco-Hernández | Mathematical approaches for emerging and reemerging infectious diseases: An introduction[END_REF]. In estimating the unmeasured proportion of infected snails and the unknown snail-human transmission rate, our framework allows a better understanding of the schistosomiasis dynamics and better implementation of snail control strategies. After the preliminaries on the observers (Section 2), the bulk of this part is described in Section 3.

The second focal point of the paper is to address schistosomiasis transmission on a larger scale. Indeed, schistosomiasis, also known as a disease of poverty, is endemic in developing countries where resources are very limited and snail habitats (ponds, marshland, fresh-water sources, ...) are widespread. And so, determining which snail habitats to target for biological control or which human age-group(s) or patch(es) to target for mass drug administration, is a major problem. To deal with this question, we consider a heterogeneous schistosomiasis model, derived originally by [START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF] (see also [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF][START_REF] Mari | Heterogeneity in schistosomiasis transmission dynamics[END_REF]). The model consists of human hosts, structured into l groups or patches, interacting with snails, structured into m natural habitats (ponds, fresh-water, etc). Assuming that we could measure the worm burden in hosts in the l host group -this is a routine technique available in most countries -our method provides a tool to estimate the proportion of infected snails in each snail habitat, thereby providing a map to optimize resource allocation in controlling schistosomiasis. This part is treated in Section 4.

Preliminaries

In this section we briefly introduce the notion of observer. Let us consider a phenomenon whose evolution is captured by a dynamical system:

ẋ = f (x, θ) (1)
where x is the vector representing state variables and θ the vector of parameters of the considered model. Suppose that only some or a combination of state variables are available as measurements.

In control theory nomenclature, this quantity is called measurable output and is denoted by

y = h(x)
where h is a smooth function of the state x.

Hence the input-output system that accounts for the dynamics of the system along with the available data is given by: ẋ

= f (x, θ), y = h(x). (2)
An observer is an auxiliary dynamical system that takes into account the available output into the dynamics of the system and that converges asymptotically to the original system [START_REF] Abdelhedi | Extended observer to estimate the spreading of contagious disease[END_REF]. Particularly, if the vector of parameters θ is known, an observer of System (1) has the following form: ż(t) = g(ẑ(t), y(t), θ),

x = k(ẑ(t), y(t), θ). (3)
The challenge in designing an observer consists of finding the two functions g and k in such a way that the solutions of System (3) converge exponentially to that of System (2), irrespective of initial conditions. That is, the functions g and k must satisfy, for any initial condition x(t 0 ):

||x(t) -x(t)|| ≤ e -λt ||x(t 0 ) -x(t 0 )||,
for some λ ∈ IR + . When System (1) is linear, the problem of finding an observer is completely solved by the Luenberger observer [START_REF] Luenberger | An introduction to observers[END_REF], and the possibility of achieving arbitrary fast estimation is equivalent to the so-called observability property (see below Definition 2.3 and Theorem 2.1). If the considered system is non-linear, typically, there is no "off the shelf" method to design an observer (see [START_REF] Besanc ¸on | Nonlinear observers and applications[END_REF] and the references therein for an overview on the observer theory, and [START_REF] Iggidr | Encyclopedia of life support systems (eolss)[END_REF] for a short survey of observability and observers with applications to various life-support-systems).

Throughout the rest of the paper, we will use some terminologies and notions that we recall in the rest of section. When the vector of parameters θ is not fully known, it may still be possible to design an observer for System [START_REF] Abdelhedi | Observer design for a class of nonlinear piecewise systems. Application to an epidemic model with treatment[END_REF]. These type of observers are defined as follows:

Definition 2.1. An adaptive observer for System (2) is an observer for joint estimation of the state x(t) and the parameter θ.

An (exponential) adaptive observer for System (2) is a dynamical system of the form

     ż(t) = g 1 (ẑ(t), θ(t), y(t)) θ(t) = g 2 (ẑ(t), θ(t), y(t)) x(t) = k(ẑ(t), θ(t), y(t)) (4)
such that the solutions of ( 4) and those of (2) satisfy x(t) -x(t) and θ(t) -θ converge (exponentially fast) to zero as time t goes to infinity.

The following two definitions and theorem will be used to prove some of the results in Section 3 and Section 4.

Definition 2.2. A square matrix A is said to be Hurwitz (or stable) if all its eigenvalues have negative real parts.

Definition 2.3 (Kalman observability rank condition [START_REF] Wonham | Linear multivariable control: a geometric approach[END_REF]). Let A be a n × n matrix and C be a q × n matrix. The pair (A, C) is said to be observable or to satisfy the Kalman observability rank condition if the matrix:

O (C,A) =          C CA CA 2 . . . CA n-1          is of rank n.
Theorem 2.1 (Pole-Shifting Theorem [50, page 61]). If the pair (A, C) satisfies the Kalman observability rank condition then it is possible to find a matrix K in such a way that the eigenvalues of the matrix A -KC are all with negative real part.

For linear systems of the form:

ẋ(t) = Ax(t) y(t) = Cx(t), (5)
Definition 2.3 provides a characterization of the observability of the corresponding linear system (5). That is, given the measurable output y, it is possible to determine the behavior of the entire system x from the system's output, for all t ∈ [0, T ]. Theorem 2.1 guarantees the existence of a matrix K such that the spectrum of A -KC belongs to the open left half-plane of l C, for an observable pair (A, C). Therefore, if the linear system (5) satisfies the Kalman observability rank condition then an exponential observer is simply given by:

ẋ(t) = Ax(t) -K(C x(t) -y(t))
since the error e(t) = x(t) -x(t) satisfies ė(t) = (A -KC)e(t) and the matrix K can be chosen in such a way that the matrix A -KC is Hurwitz. In the following section, we design an observer for a schistosomiasis model in order to estimate the number of infected snails and transmission rate from snails to human hosts.

Observer design for a schistosomiasis model

Let w be the number of female schistosomes (worms) per single host and s the proportion of infected snails. The female schistosomes per single host decays at a per capita rate γ and are replenished at a rate a. The latter process is proportional to the proportion of infected snails s. The proportion of susceptible snails 1 -s could be infected through an indirect contact with schistosomes that are excreted from hosts w, and naturally die at a rate µ. Finally, the overall dynamics of the schistosomiasis is captured by the following model:

   ẇ = -γw + as ṡ = b(1 -s)w -µs (6)
Model ( 6) is a version of MacDonald's model and is proposed first by Barbour [9] in investigating the dynamics of bilharzia -See also [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF] for a more detailed formulation of the model, and [START_REF] Anderson | Helminth infections of humans: mathematical models, population dynamics, and control[END_REF][START_REF] Basáñez | A research agenda for helminth diseases of humans: modelling for control and elimination[END_REF] for a literature review. It is easy to show that the compact set

Ω = {(w, s) ∈ IR 2 | 0 ≤ w ≤ a γ ; 0 ≤ s ≤ 1} (7)
is a positively invariant and absorbing set for System [START_REF] Antsaklis | State Feedback and State Observers[END_REF], making the solutions biologically substantiated.

It has been documented that the number of female schistosomes (worms) per single host could be measured using urine or faeces samples [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF]. Particularly, authors in [START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF] provided data of the numbers of schistosome eggs per gram of feces per host (Fig 1 in [START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF]). This is a widespread diagnosis technique for parasitological infections in general and schistosomiasis in particular. We suppose henceforth that the state variable w in Model ( 6) is measured. Thus, the measurable output is y(t) = w. By denoting x = (w, s) T , the dynamical system that describes the evolution of schistosomiasis along with the available data could be written in a compact form as follows:

   ẋ = M x + U (y) x y = C x, (8)
where M = -γ a b -µ , U (y) = 0 0 0 -by , and C = (1, 0). As many models representing a complex system, some of the parameters of Model (8) are unknown. For instance, the death rates of schistosomes γ, and infected snails µ, respectively, could be measured by biological considerations (see [START_REF] Colley | Immunology of human schistosomiasis[END_REF]). However, transmission parameters are usually unknown, and the existing proposed values are ad hoc.

In the next two subsections, we design observer to estimate the state variable s in two cases: when all the parameters are known, and when only one of the transmission parameters is known. Actually, the proposed observers provide tools to estimate the state variable s without using the above-mentioned transmission parameters a. Moreover, we use this observer to estimate this parameter.

State variable estimate when all parameters are known

The following theorem provides an observer (state-estimator) using an output-injection for the system (8).

Theorem 3.1. An exponential estimator for the MacDonald's Bilharzia model [START_REF] Antsaklis | State Feedback and State Observers[END_REF] with measurable output y(t) = w(t) is given by

ẋ = M x + U (y) x -K (C x -y) (9)
where

K T = (k 1 , a + b) for any positive constant k 1 .
Proof.

It suffices to show that the trajectories of System (9) converge exponentially to those of System (8) for any initial conditions of System (9) . To do so, we show that the equation describing the evolution of the error e(t) = x -x converges exponentially to zero. Indeed,

ė = ẋ - ẋ = M x + U (y) x -K (C x -y) -M x -U (y) x = (M -KC)e + U (y)e ( 10
)
Let us consider the Lyapunov function L = 1 2 e T e. The derivative of L along the trajectories of ( 10) is:

L = 1 2 ėT e + 1 2 e T ė = 1 2 [(M -KC)e + U (y)e] T e + 1 2 e T [(M -KC)e + U (y)e] = 1 2 e T (M -KC) T + e T U T (y) e + 1 2 e T [(M -KC)e + U (y)e] = 1 2 e T (M -KC) T + (M -KC) e + 1 2 e T U T (y) + U (y) e (11)
Using the expressions of M , K, C and U given in [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF], Equation ( 11) leads to:

L(e) = e T -k 1 -γ 0 0 -µ -by e ≤ 0,
since y ≥ 0. Moreover, L is a quadratic negative-definite function. It follows that, by Lyapunov's theorem all trajectories of (10) approach the origin with an exponential convergence rate as t goes to infinity. We deduce therefore that System ( 9) is an exponential observer of System (8).

Remark 3.1.

In Theorem 3.1 we chose K T = (k 1 , k 2 ) with k 2 = a + b. However K could be any vector that makes M -KC Hurwitz.

We carry out numerical simulations to showcase how the overall behavior of the solutions of the original system (8) compared to the observer [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF]. Unless otherwise stated, we use the following parameters: γ = 0.05 day -1 , µ = 0.04 day -1 , a = 2 female worm × day -1 and b = 0.01 female worm -1 × day -1 . The initial conditions (w(0), s(0)) = (3, 0.3) and ( ŵ(0), ŝ(0)) = (1, 0.1) are used for the systems ( 8) and ( 9), respectively. Figure 1 displays the dynamics of infected snails s (red), given by Model ( 6), and its estimated value, ŝ (blue), obtained using the observer [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF]. As anticipated, the curves of ŝ and that of s(t) merge rather quickly (convergence-time ≈ 7 days). Moreover, the merger occurs way before the two systems reach the equilibrium value (dashed green). It is worthwhile to note that the trajectories given by the observer (9) converge to those of Model ( 6), for any initial conditions of the observer system. This avoids the often difficult task of estimating the initial conditions to integrate epidemic systems.

Figure 1: Coupled dynamics of the proportion of infected snails s(t) (dashed red) as given by System (8) and its estimate ŝ(t) (blue) as given by the observer (9).

Adaptive observer: A tool to estimate both state variables and parameters

In Section 3.1, we designed an observer to estimate a state variable, namely the proportion of infected snails s, given that all the parameters of the model are known. Indeed, the observer (9) depends on all the parameters of Model ( 6). However, unfortunately the values of some parameters in epidemic models are usually not known. A perennial problem in epidemic models is the estimation of the transmission parameters. For the schistosomiasis model, it is the parameter a which represents the snail-host infection rate that is difficult to estimate [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF]. In this section, we design an observer that allows to estimate simultaneously the state variable s and the snail-host transmission rate a given that w is measured. Such observer is called adaptive observer (see for instance [START_REF] Besancon | High gain observer based state and parameter estimation in nonlinear systems[END_REF] and [START_REF] Besanc ¸on | Nonlinear observers and applications[END_REF], chap. 7). Let a, the unknown parameter. Using the linear transformation (z 1 , z 2 ) T = (w, -γw + as) T = T x, where

T = 1 0 -γ a , Model (6)
can be written, with the new coordinates as

   ż(t) = (A + B(y))z(t) + ψ(y(t))a, y = C 0 z, (12)
where

A =   0 1 -γµ -(µ + γ)   , B(y) = -by   0 0 γ 1   , ψ(y) =   0 by   , and C 0 = (1, 0).
Let λ be any positive real number. Define the following functions of the output y := y(t):

                 k 1 (y) = 3(λ + 1) -γ -µ -b y, k 2 (y) = γ 2 + µ 2 + γµ + b 2 y 2 + (2µ + γ) by + 3λ(λ + 2) + 2 -3(λ + 1)(by + µ + γ), = -3(µ + γ) + γ 2 + µ 2 + γµ + b 2 y 2 + (2µ + γ -3 -3λ) by + 3λ(λ + 2) + 2 -3λ(µ + γ), Λ(y) = λ(λ + 1)(λ + 2)
by .

(13) The expression of Λ(y) in ( 13) is well defined since y(t) > 0 for all t ≥ 0 if y(0) > 0. Indeed, y is the measurable output and satisfies the equation ẏ = -γy + as. Moreover, we assume that the ecosystem is in an endemic situation which means that R 0 = ab γµ > 1 and implies [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF] that

y(t) → ȳ = a γ 1 - γµ ab > 0.
The following theorem gives an observer that estimates both the states variables and the transmission parameter a.

Theorem 3.2.

An exponential adaptive observer for System ( 12) is given by the following system:

   ż(t) = (A + B(y))ẑ(t) + ψ(y(t))â(t) -K(y(t))(C 0 ẑ(t) -y(t)), ȧ(t) = -Λ(y(t))(C 0 ẑ(t) -y(t)) (14)
where K(y(t)) = (k 1 (y(t)), k 2 (y(t))) T and Λ(y(t)) are defined by relations [START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF].

The proof of Theorem 3.2 uses a classical result concerning a class of time-varying linear systems, which we recall here: [START_REF] Coppel | Stability and asymptotic behavior of differential equations[END_REF], Chap. IV, Theorem 2). Let D be a constant matrix whose characteristic roots λ 1 , . . . , λ n are all simple, and let ξ i be a characteristic vector of D belonging to the characteristic root

Theorem 3.3 ([
λ i (i = 1, . . . , n). If R(t) is a continuous matrix defined for t ≥ t 0 such that ∞ t 0 |R(t)| dt < ∞,
then the equation ẋ = (D+R(t))x has a fundamental system of solutions

x 1 (t), x 2 (t), . . . , x n (t) satisfying for t → ∞ x k (t) ∼ e λ k t ξ k , (k = 1, . . . , n).
The symbol ∼ in Theorem 3.3 denotes the asymptotic equivalence.

Proof of Theorem 3.2.

Recall that the goal is to show that the error between System (12) and System ( 14) converges exponentially to zero. The estimation error e(t) = (e 1 (t), e 2 (t), e a (t)

) T = (ẑ 1 (t) -z 1 (t), ẑ2 (t) -z 2 (t), â(t) -a(t)) T satisfies ė(t) =    -k 1 (y(t)) 1 0 -k 2 (y(t)) -γµ -bγy(t) -γ -µ -by(t) by(t) -Λ(y(t)) 0 0    e(t) := M (y(t)) e(t). (15)
The time-varying matrix M (y) could be written as

M (y) = Ã(y) -K(y) C, with Ã(y) =    0 1 0 -γµ -bγy -γ -µ -by by 0 0 0    , C = 1 0 0 , and K(y) =    k 1 (y) k 2 (y) Λ(y)    .
For any fixed y > 0, it is not difficult to show that the pair ( Ã(y), C(y)) is observable. That is, it satisfies the Kalman observability rank condition (see Definition 2.3). Therefore, by the pole-shifting theorem (see Theorem 2.1), it is possible to find a y-dependent gain K(y) such that sp(M (y)) = {-λ-2, -λ-1, -λ}, where sp(M (y)) denotes the spectrum of M (y) and λ is any positive real number. This gain K(y) is given by the relations [START_REF] Bethony | Exposure to schistosoma mansoni infection in a rural area in brazil. ii: household risk factors[END_REF] and can be computed using for instance Ackermann's formula (see [START_REF] Antsaklis | State Feedback and State Observers[END_REF] page 382):

K(y) = ( Ã(y) + λI 3 ).( Ã(y) + (λ + 1)I 3 ). ( Ã(y) + (λ + 2)I 3 ).O -1 (y).    0 0 1    ,
where I 3 is the 3 × 3 identity matrix and O(y) is the observability matrix associated to the pair ( Ã(y), C(y)) and it is given by O

(y) =    C(y) C(y) Ã(y) C(y) Ã(y) 2    .
The spectrum of the matrix M (y(t)) is fixed, and lies entirely in the left hand side of the complex plane as the eigenvalues of M (y(t)) are -λ -2, -λ -1 and -λ, with λ > 0. We must emphasize that, this, of course, is not sufficient to conclude that the solutions of ė = M (y(t))e converge to zero at an exponential rate since System (15) is a non autonomous system. Now, we will show that System (15) can be written in the form of u = (D + R(t))u, where D is a constant matrix. Under some conditions on R(t), the behavior of the solutions for these types of non autonomous systems could be obtained using Theorem 3.3 [START_REF] Coppel | Stability and asymptotic behavior of differential equations[END_REF], for example.

Let D = diag(-λ, -(λ + 1), -(λ + 2)) and P (t) the invertible time-dependent matrix satisfying

D = P -1 (t) M (y(t)) P (t).
A straightforward computation of the eigenvectors of M (y(t)), associated with the corresponding eigen-values leads to determine P (t) and thus P -1 (t). We successively obtain

P (t)= -by(t)          - 1 (λ + 1)(λ + 2) - 1 λ(λ + 2) - 1 λ(λ + 1) γ + µ -2λ -3 + by(t) (λ + 1)(λ + 2) γ + µ -2λ -2 + by(t) λ(λ + 2) γ + µ -2λ -1 + by(t) λ(λ + 1) - 1 by(t) - 1 by(t) - 1 by(t)         
, and

P -1 (t) =           - λ(λ + 1)(λ + 2)(γ + µ -λ + by(t)) 2by(t) - λ(λ + 1)(λ + 2) 2by(t) (λ + 1)(λ + 2) 2 λ(λ + 1)(λ + 2)(γ + µ -λ -1 + by(t)) by(t) λ(λ + 1)(λ + 2) by(t) -λ(λ + 2) - λ(λ + 1)(λ + 2)(γ + µ -λ -2 + by(t)) 2by(t) - λ(λ + 1)(λ + 2) 2by(t) λ(λ + 1) 2           .
We perform the change of variable u = P (t) -1 e, and the new variable u(t) satisfies the time-varing linear equation

u = (D + R(t)) u, (16)
where

R(t) = dP -1 (
     . Therefore, it follows that R(t) ∞ = y (t) y(t) |(3(λ + 1)(λ -by(t)) + 3λ + 1)| .
The solutions of of System (6) evolve in the compact set Ω defined in [START_REF] Aronna | Interval observer for uncertain time-varying SIR-SI model of vector-borne disease[END_REF], so y(t) (= w(t)) is a positive and bounded function of t. We choose λ large enough in such a way that λ -by(t) ≥ 0. Thus,

R(t) ∞ = y (t) y(t) (3(λ + 1)(λ -by(t)) + 3λ + 1) = y (t) y(t) 3λ 2 + 6λ + 1 -3b(λ + 1)y(t) . R(t) ∞ = (3λ 2 + 6λ + 1) y (t) y(t) -3b(λ + 1) |y (t)| .
Since System ( 6) is a planar cooperative system there exists T ≥ 0 such that y(t) is nonincreasing or nondecreasing on t ≥ T (see [START_REF] Smith | Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems[END_REF], Chapter 3, Theorem 2.2, page 43, [START_REF] Hirsch | Handbook of Differential Equations[END_REF], Theorem 3.21, or [START_REF] Hirsch | Systems of differential equations which are competitive or cooperative. I: Limit sets[END_REF], Proposition 2.5). This implies that y (t) changes sign at most once (for t = T ). Using Chasles' identity, we have

∞ 0 R(t) ∞ dt = T 0 R(t) ∞ dt + ∞ T R(t) ∞ dt, and ∞ T R(t) ∞ dt = ± ∞ T y (t) y(t) 3λ 2 + 6λ + 1 -3b(λ + 1)y(t) dt. ( 17
)
The sign in ( 17) is positive when y (t) ≥ 0 and negative otherwise. Thus, we obtain, using the fact that y(t) tends to ȳ = a γ 1 -γµ ab as t tends to +∞

∞ T R(t) ∞ dt = ± 3λ 2 + 6λ + 1 ∞ T y (t) y(t) dt -3b(λ + 1) ∞ T y (t)dt = ± 3λ 2 + 6λ + 1 log ȳ y(T ) -3b(λ + 1) (ȳ -y(T )) < ∞.
Hence, The observer in Theorem 3.2 gives an estimation ẑ of z and, â of a. Hence, using the linear and invertible transformation z = T x, we obtain an estimation of x. Indeed, in the original coordinates x = ( ŵ, ŝ), the adaptive observer ( 14) is given by:

∞ 0 R(t) ∞ dt < ∞,
           ẋ(t)=   -γ â(t) b -µ -by(t)   x(t) -       1 0 γ â(t) 1 â(t)    K(y(t)) - ŝ(t) â(t) 0 Λ(y(t))    (C 0 x(t) -y(t)), ȧ(t)=-Λ(y(t))(C 0 x(t) -y(t)) (18)
It is worthwhile noting that the estimators ( 14) and ( 18) are completely independent of a. Figure 3 provides the dynamics of z 2 (t) (red), solution the shifted system [START_REF] Besancon | High gain observer based state and parameter estimation in nonlinear systems[END_REF] and that of its estimate ẑ2 (t) (dashed blue), obtained using the corresponding observer [START_REF] Bichara | Étude de modèles épidémiologiques: Stabilité, observation et estimation de paramètres[END_REF]. The dynamics of the original variable s(t) and its estimate ŝ(t) are represented in Figure 4. The estimate ŝ(t) is obtained using the adaptive observer, in original coordinates. That is, System [START_REF]Neglected tropical diseases, Centers for Disease Control and Prevention[END_REF]. Also, as expected from Theorem 3.2, ẑ2 (t) (Figure 3) and ŝ(t) (Figure 4) converge exponentially to z 2 (t) and s(t) respectively. Moreover, the efficiency of the designed observer is showcased in Figure 4, as it can be noticed that the estimated state ŝ(t) converges toward the state s(t) much earlier than they both reach the equilibrium value (dashed green line).

Figure 3: Dynamics of the state variable z 2 (t) (dashed red) given by system [START_REF] Besancon | High gain observer based state and parameter estimation in nonlinear systems[END_REF] and its corresponding estimate ẑ2 (t) (blue) delivered by the observer [START_REF] Bichara | Étude de modèles épidémiologiques: Stabilité, observation et estimation de paramètres[END_REF]. The initial conditions are such that (z 1 (0), z 2 (0)) = (3, 0.3) and (ẑ 1 (0), ẑ2 (0)) = (1, 0.1).

Figure 4: Dynamics of proportion of infected snails s(t) (dashed red) and its estimate ŝ(t) (blue) delivered by the observer [START_REF]Neglected tropical diseases, Centers for Disease Control and Prevention[END_REF]. The estimation ŝ(t), of s(t) is obtained without using the unknown parameter a. The initial conditions are such that (w(0), s(0)) = (3, 0.3) and (ẑ 1 (0), ẑ2 (0)) = (1, 0.1). The green dashed line corresponds to the equilibrium value of s(t), using Model (6).

Figure 5: Dynamics of proportion of infected snails s(t) (dotted red curve) and its estimate ŝ(t) (blue) when the measurements y(t) are corrupted by noise. The estimation ŝ(t), of s(t), is obtained without using the unknown parameter a, using the adaptive observer [START_REF]Neglected tropical diseases, Centers for Disease Control and Prevention[END_REF].

Remark 3.2. The observer ( 9) and the adaptive observer (18) both estimate the proportion of infected snails s(t), and for our model and baseline parameter values, numerical simulations indicate that the estimations given by both observers converge to the state variables of Model [START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF] with approximately the same time period of around 5 days (see Figure 1 and Figure 4). However, the adaptive observer (18) allows to estimate not only the unmeasured state variables but also the unknown snail-man transmission parameter a. Moreover, numerical simulations suggest that the adaptive observer (18) is less sensitive to measurement noise as it can be seen by comparing Figure 5 and Figure 2.

Remark 3.3. The state and parameters estimation problem for the MacDonald's model [START_REF] Antsaklis | State Feedback and State Observers[END_REF] has been addressed in [START_REF] Tendeng | Etude de modèles de transmission de la Schistosomiase: Analyse mathématique, reconstruction des variables d'état et estimation des paramètres[END_REF] using high-gain observer (we refer to [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF] for a comprehensive introduction of a highgain observer). However, the use of high-gain observer requires the construction of a prolongation of the vector field that has to be globally Lipschitz not only on the set Ω but on the whole IR 2 . It is also necessary to extend the diffeomorphism (defining the coordinates change) to the whole IR 2 and this extended diffeomorphism must also be globally Lipschitz [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF]. The explicit construction of these extensions is not an easy task. The high-gain observer may fail to converge if these Lipschitz extensions are not used as it has been shown in [START_REF] Guiro | On the stock estimation for some fishery systems[END_REF].

In Theorem 3.1, we provided an estimation of the proportion of infected snails (s(t)) when all parameters are supposedly known. And, in Theorem 3.2, we designed an adaptive observer to estimate simultaneously both the state variable s(t) and a key unknown parameter -the snail-to-human transmission rate. However, these estimations correspond to our particular set of data, and thus may not be used in the context of the actual disease control. Indeed, these estimations are based on Macdonald's model, a crude representation of the dynamics of Schistosomiasis. To address parts of the coarseness Figure 6: Comparison between the estimation of the unknown parameter a, as limit of the solution â(t) (blue) of the estimator ( 14) or [START_REF]Neglected tropical diseases, Centers for Disease Control and Prevention[END_REF] and the value of a, taking in the simulations of Model [START_REF] Antsaklis | State Feedback and State Observers[END_REF]. It can be observed that the estimation convergence is quite fast. of Macdonald's model, we provide in the following subsection an estimation when the infection process is represented by a function f (w, s), instead of an exponential distribution.

State Estimate when the snail-host infection process is not known

In estimating the proportion of infected snails, we used MacDonald's model [START_REF] Antsaklis | State Feedback and State Observers[END_REF], for which the snailhuman infection term is linear, namely "a s". Here, we propose a tool that allows to estimate the proportion of infected snails when the infection process is captured by any function f (w, s) whose analytic expression is not necessarily known. Indeed, in this case, the modified MacDonald's model could be written as

         ẇ = -γw + f (w, s) ṡ = b(1 -s)w -µs y = w, (19)
or equivalently,

ẋ = Gx + Ef (x) + g(x), y = Cx (20a) (20b)
where is an exponential observer for system [START_REF] Colley | Immunology of human schistosomiasis[END_REF], that is, x(t) converges exponentially towards x(t).

x = (w, s) T , G =   -γ 0 b -µ   , E =   1 0   , C = (1,

Proof.

The observer construction is inspired from [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF]. Let us decouple the System (19) from the unknown input f (w, s) = f (x). Indeed, from Equation (20b), it follows that

ẏ = C ẋ = C (Gx + Ef (x) + g(x)) = CGx + CEf (x) + Cg(x).
Since CE=1, we deduce that f (x) = ẏ -CGx -Cg(x). Thus, Equation (20a) will be decoupled from the unknown input f (x). Indeed, we obtain:

ẋ = Ḡx + E ẏ + (I 2 -EC)g(x).
Now it suffices to investigate the dynamics of the error e(t) := x(t) -x(t) between Systems ( 20) and [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF]. We have:

ė = ẋ - ẋ = ( Ḡ -KC)e + ĝ -(I 2 -EC)g(x) = ( Ḡ -KC)e -by 0 0 0 1 e = -k 1 0 b -k 2 -µ -by e. (22)
We choose k 2 = b, we then have

d dt ( e(t) 2 2 ) = 2(-k 1 e 2 1 -(µ + b y(t))e 2 2 ). Since y(t) ≥ 0, it follows that d dt ( e(t) 2 2 ) ≤ -2(k 1 e 2 1 + µe 2 2 ) ≤ -2 min(k 1 , µ) e(t) 2 2 .
Hence, e(t) 2 2 ≤ e -2 min(k 1 ,µ) t e(0) 2 2 which shows that the estimation error e(t) converges exponentially towards zero. This ends the proof. Remark 3.4. System ( 21) is called observer with unknown input [START_REF] Bichara | On the estimation of sequestered infected erythrocytes in Plasmodium falciparum malaria patients[END_REF][START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF][START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF][START_REF] Hou | Disturbance decoupled observer design: A unified viewpoint[END_REF].

Observer design for a patchy Schistosomiasis model

In this section, we propose to design an observer for Schistosomiasis model that accounts for individuals and ponds heterogeneity. Indeed, we revisit a dynamical model proposed in [START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF] for which individuals are structured into l patches or groups that interact with the snails that are dispatched in m freshwater contact sites or pond locations. Using the same setup as the homogeneous case and by denoting w = (w 1 , w 2 , . . . , w l ) the vector corresponding to the number of female schistosomes (worms) per host (for l patches or groups) and s = (s 1 , s 2 , . . . , s m ) the vector corresponding to the proportion of infected snails (for m freshwater sites or pond locations). Hence, following [START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF], the dynamics of Schistosomiasis that results from the interactions of m groups of individuals into l freshwater sources is given by:

   ẇ = -diag(γ)w + As ṡ = diag(1 -s)Bw -diag(µ)s (23)
where A = (a ij ) ∈ M l,m (IR) is a nonnegative l×m matrix, which captures the snail to man transmission, and B = (b ij ) ∈ M m,l (IR) a nonnegative m × l matrix, which captures the man to snail transmission. More precisely, a ij is the transmission from snails in Pond j to humans of Patch i and b ij is the transmission from humans of Patch i to snails of Pond j. As for the homogeneous case, γ ∈ IR l + is the vector of death rates of female schistosome in the human host populations, µ ∈ IR m + is the vector of death rates of snails in the water site. By denoting 1, the vector (1, • • • , 1) T of IR m , it could be shown easily (see [START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF] for instance) that the set

D = {(w, s) ∈ IR l+m + | 0 ≤ w ≤ diag(γ) -1 A 1; 0 ≤ s ≤ 1}
is an absorbing and positively invariant set for system [START_REF] Diaby | Observer design for a schistosomiasis model[END_REF], thereby making the solutions of (23) biologically grounded. A detailed formulation of Model [START_REF] Diaby | Observer design for a schistosomiasis model[END_REF] and its dynamics can be found in [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF][START_REF] Mari | Heterogeneity in schistosomiasis transmission dynamics[END_REF][START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF] for example. We assume that the measurable output is the vector w. That is, we measure the worms in humans host in each group. Noting that diag(1 -s)Bw -diag(µ)s = B w -diag(Bw) s -diag(µ)s, and since the measurable output is y = w, Model [START_REF] Diaby | Observer design for a schistosomiasis model[END_REF] along with the measured data could be written as

   ẋ = M x + U (y) x y = C x, (24)
where

x = (w, s) T , M = -diag(γ) A B -diag(µ) , U (y) = 0 0 0 -diag(By) , and
C =        1 0 • • • 0 0 1 • • • 0 . . . . . . 0 0 • • • 1 l columns 0 • • • 0 0 • • • 0 . . . . . . 0 • • • 0 m columns        .
The following theorem gives an exponential observer for System [START_REF] Feng | Estimation of parameters governing the transmission dynamics of schistosomes[END_REF].

Theorem 4.1.

An exponential state estimator for System (24) is given by

ẋ = M x + U (y) x -K (C x -y) (25)
where Thus, it is possible to find K 1 and K 2 that make V is negative definite. A simple choice is K 2 = A T + B and K 1 > 0.

K = K 1 K 2 with K 2 = A T + B,
The observer obtained in Theorem 4.1 enables the estimation of the proportion of infected snails in each freshwater site. This is particularly useful in providing a tool to better control Schistosomiasis, such as biological control of snails [START_REF] Swartz | Infection with schistosome parasites in snails leads to increased predation by prawns: implications for human schistosomiasis control[END_REF] or in allocating resources (such as the use of chemical molluscicides [START_REF] King | Historical perspective: snail control to prevent schistosomiasis[END_REF]) at appropriate snails' freshwater habitats according to their infected snails importance.

Conclusion and discussions

In this paper, we investigated tools and methods to estimate state variables and parameters of a class of schistosomiasis models using elements of control theory. Particularly, we address two fundamental problems in modeling and control of schistosomiasis transmission. Namely, we propose a method to estimate the proportion of infected snails in a given snail habitat and the snail-host transmission parameter. We revisited the classical MacDonald's model [START_REF] Anderson | Helminth infections of humans: mathematical models, population dynamics, and control[END_REF][START_REF] Barbour | Modeling the transmission of schistosomiasis: an introductory view[END_REF][START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF] for schistosomiasis and used the fact that a key state variable could be measured, namely the worms burden in humans. Our method consists of using elements of system theory in designing an auxiliary dynamical system, called observer, that converges exponentially to the original model.

First, we estimated the proportion of infected snails, assuming that all parameters are known (Theorem 3.1) and when the snail-host transmission parameter is unknown (Theorem 3.2) in a homogeneous system. Moreover, using the latter observer, we estimated the unknown parameter. We then compared the asymptotic value of the obtained parameter with that a constant valued used in the dynamics of Schistosomiasis model (Figure 6).

Second, we considered a heterogeneous schistosomiasis model with an arbitrary number of age-group in the structure of human hosts and where the snails are distributed into m water sources, considered in [START_REF] Barbour | Macdonald's model and the transmission of bilharzia[END_REF][START_REF] Mari | Heterogeneity in schistosomiasis transmission dynamics[END_REF][START_REF] Riveau | A Transmission Model of Bilharzia A Mathematical Analysis of an Heterogeneous Model[END_REF]). We derived a method to estimate the proportion of infected snails in each snail habitat (Theorem 4.1). This result, we hope, could help in better designing effective public health policies, particularly in allocating spare resources where they are needed the most, with respect to host groups or snail habitats.

In both the homogeneous and heterogeneous cases, we estimated the snail-human transmission rate. However, the human-snail transmission b is also unknown. This parameter become the WAIFW (Who Acquires Infection From Whom [START_REF] Anderson | Infectious Diseases of Humans. Dynamics and Control[END_REF]) matrix B for the heterogeneous case. This is very challenging as the cross term is nonlinear and is the subject of an ongoing investigation. In this paper, the observers are designed with the assumption that the measurable output is a function continuous in time. However, in real life, one can hope at best to have a limited number of samples of measurements at discrete times. It would be interesting to design observer for epidemic models -that are continuous in timefor which the measurable output is discrete in time.

Figure 2 :

 2 Figure 2: Coupled dynamics of the proportion of infected snails s(t) (dotted red curve) as given by System (8) and its estimate ŝ(t) (blue) as given by the observer (9) when the measurements y(t) are corrupted by noise.

Figure 7 :

 7 Figure 7: Comparison between the estimation of the unknown parameter a, as limit of the solution â(t) (blue curve) of the estimator (14) when the measurements y(t) are corrupted by noise and the value of a, taking in the simulations of Model (6). It can be observed that the estimation convergence is quite fast.

Theorem 3 . 4 .

 34 0), and g(w, s) Let z(t) = (z 1 (t), z 2 (t)) T , Ḡ = (I 2 -EC)G, K = (k 1 , k 2 ), with k 1 > 0, k 2 = b > 0 and ĝ(t) = (0, -byz 2 ) T . The following system ż = ( Ḡ -KC)z + (K + ( Ḡ -KC)E)y + ĝ, x = z + Ey[START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF]

  and K 1 any l × l square matrix whose components are all positive.Proof. Let e = x -x, the error between System (25) and System[START_REF] Feng | Estimation of parameters governing the transmission dynamics of schistosomes[END_REF]. The evolution of the errorequation is ė = (M + U (y) -K C) e. Considering the function V = 1 2e T e, we obtainV = -e T 1 (K T 1 + diag(γ))e 1 -e T 2 (diag(µ) + diag(By))e 2 + e T 2 (A T + B -K 2 )e 1 .

Acknowledgements

We are grateful to two anonymous referees and the handling editor for valuable comments and suggestions that led to an improvement of this paper.