Malika Kheireddine
email: malika.kheireddine@kaust.edu.sa

David Antoine

Vincenzo Vellucci

Bernard Gentili

Ras

M D Antoine

H Chami

F Claustre

A D'ortenzio

G Morel

Becu

F Louis

J Ras

E Roussier

A J Scott

D Tailliez

S B Hooker

P Guevel

J F Deste

C Dempsey

D Adams

Diel changes of the optical backscattering coefficient of oceanic particulate matter determined from diel changes in apparent optical properties: a case study in the Mediterranean Sea (BOUSSOLE site)

published or not. The documents may come     Diel changes of the optical backscattering coefficient of oceanic particulate matter determined from diel changes in apparent optical properties: a case study in the Mediterranean Sea (BOUSSOLE site)

INTRODUCTION

Recording from space the changes of the ocean reflectance from dawn to dusk becomes a realistic target with the advent of ocean color sensors on geostationary satellites. The Korean Geostationary Ocean Color Imager (GOCI, [START_REF] Choi | GOCI, the world's first geostationary ocean color observation satellite, for the monitoring of temporal variability in coastal water turbidity[END_REF]) has opened the way since 2010, followed by its successor GOCI-II in 2020 [START_REF] Ahn | Missions and user requirements of the 2nd Geostationary Ocean Color Imager (GOCI-II)[END_REF]3], and other instruments are announced for launches in the 2020 decade (e.g., the NASA Geosynchronous Littoral Imaging and Monitoring Radiometer [START_REF]NASA Targets Coastal Ecosystems with New Space Sensor[END_REF]). Such missions can typically record the ocean reflectance hourly throughout the course of a day, either targeting specific areas or embracing the full Earth disk [5]. Taking full advantage of these high-frequency top-of-atmosphere observations requires numerical inversion techniques capable of retrieving geophysical properties of interest with a consistently high accuracy under the changing 45 observation conditions encountered during a day [6].

Research Article

way is opened for determination of net community production at large scale from satellite observations. In open ocean Case I waters (as per [17]), b bp typically increases from sunrise to about 2-4 h before sunset, and then decreases until next sunrise [18]. The amplitude of this cycle is determined by changes in concentrations, size structure, and refractive index of the different populations of particles, and of phytoplankton in particular. Although the term diel changes normally refers to changes involving a 24-hour period that includes a day and the adjoining night, it will be used hereafter for the sake of conciseness to refer to the changes of b bp that occur during daylight only.

The questions to be addressed are [START_REF] Choi | GOCI, the world's first geostationary ocean color observation satellite, for the monitoring of temporal variability in coastal water turbidity[END_REF] whether the b bp diel cycle observed in situ generates a measurable diel variability in reflectance and, assuming the answer is positive, (2) whether the diel b bp cycle can be retrieved from existing numerical inversion algorithms of the reflectance spectrum also measured in situ at different times of the day and, (3) whether ( 2) is still feasible when the reflectance is derived from satellite ocean color radiometry. This paper only addresses points (1) and [START_REF] Ahn | Missions and user requirements of the 2nd Geostationary Ocean Color Imager (GOCI-II)[END_REF].

It is stressed here that correctly retrieving the amplitude of the diel cycle is somehow little dependent on the absolute accuracy of b bp retrievals, provided that the associated uncertainties do not change significantly for retrievals carried out from reflectances measured at different times of a day (so for various Sun elevations). These questions were addressed by using high-frequency (every 15 min) field observations of b bp ; the irradiance reflectance, R; the diffuse attenuation coefficient for downward irradiance, K d ; and the remote-sensing reflectance, R rs , performed at the BOUée pour l'acquiSition d'une Série Optique à Long termE site in the Mediterranean Sea (BOUSSOLE [19,20]). Contrary to what was done by Kheireddine and Antoine [18] using the same dataset, average seasonal patterns are not aimed at here. A few case studies have been selected to answer question [START_REF] Choi | GOCI, the world's first geostationary ocean color observation satellite, for the monitoring of temporal variability in coastal water turbidity[END_REF] above, and to illustrate the first challenge in answering question [START_REF] Ahn | Missions and user requirements of the 2nd Geostationary Ocean Color Imager (GOCI-II)[END_REF].

It must be kept in mind that no fully analytical bio-optical model exists that would allow making fully explicit the dependence of reflectance on IOPs plus the environment conditions (the Sun zenith angle in particular). Therefore, our goal here is precisely to assess whether existing semi-analytical algorithms still can be used for deriving diel variability of b bp from diel variability of reflectance. We do not claim that we have selected the most appropriate models to perform these retrievals. Instead, we have selected among the most popular methods in order to study their relevance to this specific objective. and Dana [START_REF] Maffione | Instruments and methods for measuring the backward-scattering coefficient of ocean waters[END_REF]) using the beam attenuation coefficient mea-155 sured in parallel (see below) and the total absorption coefficient 156 derived from inversion of K d and R (Eqs. 12 and 13 in [START_REF] Morel | Bio-optical properties of high chlorophyll case 1 waters and of yellow-substance-dominated case 2 waters[END_REF]).

DATA AND METHODS

157

The b bp (λ) coefficient is derived from the corrected β(λ)(140)

158
as follows [START_REF] Boss | Relationship of light scattering at an angle in the backward direction to the backscattering coefficient[END_REF]:

159 b bp (λ) = 2π χ p (β(λ)(140) -β w (λ)(140)) m -1 , (1)
where χ p = 1.13 [START_REF] Dana | Determining the backward scattering coefficient with fixed-angle backscattering sensors? Revisited[END_REF] and where β w (λ)(140), the contribution 160 of pure seawater scattering at 140 (3)

This method assumes that non-water absorption, by colored dissolved organic matter in particular, is negligible at 650 nm [START_REF] Bricaud | Absorption by dissolved organic matter of the sea (yellow substance) in the UV and visible domains[END_REF]. The instrument bodies are covered with copper tape.

Source and detector windows are equipped with copper rings and are cleaned about every 2 weeks by divers using soft brushes.

These measures have proven efficient in preventing biofouling in most cases. Possible remaining corrupted data are identified from the comparison of data collected before and after cleaning operations. They are eliminated and not used here.

The same instruments are deployed on the monthly casts, and their measurements are used to correct for possible calibration bias or drifts in the buoy transmissometer data.

Particulate Backscattering Ratio

The backscattering ratio, bbp (λ), is derived from b bp (λ) and c p (650) as follows:

bbp (λ) = b bp (λ) c p (650) λ 650 n (adimensional), (4)
R = E u (0 -) E d (0 -) (adimensional), (5)
where E d (0 -) is the downward irradiance just below the sea sur-239 face. It is calculated as E s reduced by transmission across the air-240 water interface, i.e., E s times 0.97 [START_REF] Austin | The remote sensing of spectral radiance from below the ocean surface[END_REF]. The upward irradiance 241 just beneath the surface, E u (0 -), is obtained from

242 E u (0 -) = E u (z 4 )e K u Z 4 µW cm -2 nm -1 , (6)
where z 4 is the depth of the shallowest buoy arm (nominally 243 at 4 m) and K u is the diffuse attenuation coefficient for the 244 upward irradiance. K u is computed from the measurements of 245 E u collected at the two buoy measurement depths:

246

K u = - 1 z l n E u (z 9 ) E u (z 4 ) (m -1 ), (7)
where z 9 is the depth of the deepest buoy arm (nominally at 9 m) 247 and z = z 9 -z 4 .

248

The remote-sensing reflectance, R rs , is computed as 249

R rs = L w E s (sr -1 ), (8)
where E s is the above-surface downward solar irradiance and the 250 water-leaving radiance, L w , is computed as

251 L w = L u (0 -) 1 -ρ n 2 (µW cm -2 nm -1 sr -1 ), ( 9
)
where ρ is the water-air reflection coefficient (set to 0.043), and 252 n is the refractive index of seawater (1.34).

253

L u (0 -) is obtained as 254 L u 0 - = L u (Z 4 ) e K Lu Z 4 (µW cm -2 nm -1 ), ( 10
)
Research Article where 255

K Lu = - 1 z ln L u (z 9 ) L u (z 4 ) m -1 . ( 11
)
The diffuse attenuation coefficient for the downward irradi-256 ance in the upper layers is computed as 257

K d = - 1 Z 9 ln E d (Z 9 ) E d (0 -) m -1 . ( 12
)
The effect of buoy and instrument self-shading on underwa- 265

K d ,n = K d * µ d (θ s , [Chl]) m -1 , ( 13
)
where µ d is the average cosine of the downward irradiance and This approach combines two equations [START_REF] Gordon | Computed relationships between inherent and apparent optical-properties of a flat homogeneous ocean[END_REF][START_REF] Gordon | Can the Lambert-Beer law be applied to the diffuse attenuation coefficient of ocean water[END_REF][START_REF] Morel | Radiation transport within oceanic (case 1) water[END_REF] (wavelength 283 omitted):

284 K d = 1.0395 a + b b µ d (θ S , [Chl]) m -1 , (14)
and 285 R = f (θ S , [Chl]) b b a + b b , ( 15
)
where µ d is as described above [Eq. ( 13)] and f is a factor that 286 conceals dependencies of R on the Sun zenith angle and the 287 IOPs [START_REF] Morel | Diffuse reflectance of oceanic waters-its dependence on sun angle as influenced by the molecular scattering contribution[END_REF]. A precomputed lookup table is derived from the 288 same radiative transfer simulations as for K d .

289

Combining Eqs. ( 14) and (15) leads to

290 b b = 0.962 µ d (θ s , [Chl])K d (R/f (θ s, [Chl])) (m -1 ), ( 16
)
and 291 a = 0.962 µ d (θ s , [Chl])K d 1 -R/f (θ s, [Chl]) (m -1 ). ( 17
)
The second inversion algorithm that was used is version 6 of the quasi-analytical algorithm (QAA, [START_REF] Lee | Deriving inherent optical properties from water color: a multiband quasi-analytical algorithm for optically deep waters[END_REF]), which uses R rs as input quantity. Briefly, the total absorption spectrum is first derived from R rs , after selecting a reference wavelength (555 or 667 nm). The calculation of absorption and backscattering are propagated to other wavelengths, assuming that the particulate backscattering varies spectrally in a power-law function [START_REF] Gordon | Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery: A Review[END_REF].

The non-water absorption and backscattering coefficients are derived without any assumption about the spectral shapes of the absorption coefficients of each constituent.

The third algorithm is the GSM semi-analytical bio-optical model [START_REF] Garver | Inherent optical property inversion of ocean color spectra and its biogeochemical interpretation: 1. Time series from the Sargasso Sea[END_REF][START_REF] Maritorena | Optimal tuning of a semi-analytical model for global applications[END_REF], which also uses R rs as input quantity. GSM is based on the quadratic relationship between R rs and the absorption and backscattering coefficients from [START_REF] Gordon | Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery: A Review[END_REF], which assumes known spectral shapes for specific absorption coefficients of Whatever the algorithm, the particulate backscattering coefficient is obtained by removing the water contribution from the total b b , in the same way than for the in situ measurements (i.e., from [START_REF] Zhang | Estimating scattering of pure water from density fluctuation of the refractive index[END_REF][START_REF] Zhang | Scattering by pure seawater: effect of salinity[END_REF]).

C. Radiative Transfer Computations

The consistency between diel changes of the measured AOPs (1) Absorption and scattering coefficients of pure seawater, as well as the molecular phase function are known and prescribed in the code [START_REF] Pope | Absorption spectrum (380-700 nm) of 814 pure water. 2. Integrating cavity measurements[END_REF][START_REF] Morel | Optical properties of pure water and pure sea water[END_REF]. They are kept constant throughout the day.

(2) The non-water absorption coefficient of the water body is unknown (not measured). Therefore, it was computed following [START_REF] Morel | Bio-optical properties of oceanic waters: a reappraisal[END_REF], which uses the irradiance and the diffuse attenuation for downward irradiance to determine the total non-water absorption from [Chl] (Eqs. (8) to ( 14) in [START_REF] Morel | Bio-optical properties of oceanic waters: a reappraisal[END_REF]).

RMSE = 10 ∧    N i=1 [log 10 (m i )-log 10 (O i )] 2 N    , ( 18
) BIAS = 10 ∧ N i=1 ( log 10 (M i )-log 10 (O i ) ) N , (19)
MAE = 10 ∧ N i=1 |(log10(Mi )-log 10 (O i ))| N , ( 20
)
where N is the number of measurements, M is the retrieved value, and O the measured in situ values, respectively.

RESULTS

A The diel variability in R rs , R, and K d measurements was calculated as the relative amplitude for each parameter from sunrise to their observed maximum during the day (Fig. 2). Overall, the R values range between 15% to 55%, with a changing slope around noon whatever the wavelength (Figs. 2A-2H).

In general, the relative changes tend to increase with increasing , 4D, 4E, and4F), whereas the MOR algorithm 450 showed only a good comparison at 555 nm (Figs. 4A and4B) 451 (Table 2). Finally, the comparison between the relative amplitude calculated from modeled and from measured b bp data has been extended beyond the four case studies here selected, by applying the inversion to a larger dataset (see Section 2.D and Table 1).

The comparison shows that it is possible to obtain satisfactory values of ˜ b bp at 555 nm (Fig. 7B), whereas the performances 483 at 443 nm are poorer (Fig. 7A), although this is largely due to 484 the lower performance of MOR at 443 nm (Table 3). We also and7D). Unsurprisingly, we 488 observed a better accuracy in the green (Fig. 7D) than in the blue 489 (Fig. 7C) (Table 3). However, significant differences were observed between the b bp values retrieved from AOPs and those obtained from the in situ measurements, questioning whether the performance of the inversion algorithms are good enough to allow deriving rather small daily changes of b bp from AOPs measurements taken at different times of a day. These inaccuracies are not specific to the four case studies selected here, as if their selection would have been inappropriate. When the inversion algorithms were applied to a larger dataset, the retrieved versus measured b bp values show a significant dispersion (Fig. 4, Table 2). In particular, MOR performs better in the green (i.e., 555 nm) (r 2 = 0.67, Fig. 4B) with no noticeable degradation for low b bp values, as compared to the blue (r 2 = 0.37 instead of 0.67; Fig. 4A).

These results are questioning the possibility to obtain accurate

Research Article

A. Field Measurements and Sampling 1 . 2 .

 12 Description of the BOUSSOLE Site and Related Field Operations Essential information about the site characteristics, the measurement platforms and the instrumentation is provided in [19-21]. Only few elements relevant to the present work are reminded. The BOUSSOLE site is located in deep waters (2440 m) of the Ligurian Sea, one of the sub-basins of the north-114 western Mediterranean Sea. An important seasonality exists in 115 physical conditions, with deep mixed layers in winter (∼400 m, 116 sometimes down to the bottom), and a marked stratification 117 during summer. This seasonal dynamics drives the seasonal 118 changes in trophic levels from oligotrophy in summer to bloom 119 conditions in spring. There is, accordingly, a large range of opti-120 cal properties at this site. A buoy has been permanently deployed 121 at the BOUSSOLE site since September 2003 and operates 122 in a quasi-continuous mode, with data acquisition sequences 123 of one minute every 15 min night and day. This platform was 124 specifically designed to measure radiometric quantities in the 125 water column (nominally at 4 m and 9 m), and above the surface 126 (4.5 m above sea level), from which apparent optical properties 127 (AOPs) are derived, typically the irradiance reflectance, R, the 128 remote-sensing reflectance, R rs , and the diffuse attenuation 129 coefficient for downward irradiance, K d . IOPs, including b bp , 130 are collected at the same depths in the water column. Two sister 131 buoys equipped with the same set of instruments are used, with 132 rotation performed about every 6 months. The site is visited 133 monthly for buoy servicing, during which 0-400 m CTD casts 134 are performed for acquisition of hydrological data (conduc-135 tivity, temperature, and pressure), IOPs, and water sampling for 136 subsequent phytoplankton pigment analyses and particulate 137 absorption measurements [22]. 138 Particulate Backscattering Coefficient 139 The total angular scattering coefficient at 140 • , β(λ)(140), 140 was measured at BOUSSOLE using a Hydro-Optics, Biology, 141 and Instrumentation Laboratories (HOBI Labs) Hydroscat-4 142 backscattering meter installed at 9 m and equipped with filters at 143 443, 488, 555, and 620 nm. The instrument faceplate is made 144 of copper, which has proven to fully prevent biofouling. The 145 instrument operates at 1 Hz, so that about 60 measurements are 146 collected during each of the 1 min data collection sequences. 147 The median of these 60 measurements is used to derive a rep-148 resentative value for β(λ)(140). Dark current measurements are 149 performed on site by divers with a black neoprene cap covering 150 the instrument windows. Their average for a given 6 month 151 deployment are subtracted from the raw β(λ)(140) measure-152 ments. The β(λ)(140) values are also corrected for attenuation 153 along the measurement path (the σ (λ) correction of Maffione 154

  Discrete sampling is performed during the BOUSSOLE cruises between the surface and 200 m. Seawater samples are collected from Niskin bottles and filtered through 25 mm diameter Whatman GF/F filters (0.7 µm retention capacity), stored in liquid nitrogen during the cruise, and then transferred at -80 • C in the laboratory until algae pigment contents are measured using high-performance liquid chromatography, following [31]. The total chlorophyll a concentration ([Chl]) is computed as the sum of the concentrations of Chl a , chloro-216 phyllide a , and divinyl Chl a . The data obtained from samples 217 at 5 m and 10 m are used here. 218 6. Radiometric Quantities and Apparent Optical Properties 219 All radiometric measurements described hereafter were col-220 lected using Satlantic OCI-200 Series radiometers. They are 221 placed at the extremity of horizontal arms on the buoy for meas-222 uring the downward and the upward irradiance, E d and E u , 223 and the upwelling nadir radiance, L u , at 4 and 9 m, and also on 224 the top of the buoy (4.5 m a.s.l.) to measure the above-surface 225 downward irradiance, E s . These radiometers are equipped with 226 seven cosine collectors at different spectral bands (412, 443, 227 490, 510, 555, 665, and 680 nm). The underwater instruments 228 containers are covered with copper tape to protect them against 229 biofouling, and divers clean the optical surfaces about every 230 2 weeks. These radiometers operate at 6 Hz so that about 360 231 measurements are collected during each of the 1 min data collec-232 tion sequences. The median value of each 1 min record is kept 233 as representative of the measurement. An average dark current 234 value is computed from data collected during nighttime and is 235 subtracted from all measurements. 236 From these measurements, the irradiance reflectance is deter-237 mined as (wavelength dependence omitted in what follows) 238

  258ter radiometers was determined with the SimulO 3D backward 259 Monte Carlo code[START_REF] Leymarie | ProVal: a new autonomous profiling float for high quality radiometric measurements[END_REF] and corrected through lookup tables 260 generated for individual instrument geometries[START_REF] Białek | Monte Carlo-based quantification of uncertainties in determining ocean remote sensing reflectance from underwater fixed-depth radiometry measurements[END_REF], whereas 261 extrapolation errors on L u (0 -) and E u (0 -) are corrected 262 trough lookup tables [20] generated with HydroLight code [35]. 263 Because K d significantly depends on the Sun zenith angle, θ s , 264 its normalized values are also computed as follows:

  266 is taken from precomputed lookup tables generated through 267 HydroLight radiative transfer computations including Raman 268 emission and a [Chl]-varying volume scattering function (VSF) 269 for particles [36]. The full dependence of K d on the Sun zenith 270 angle and IOPs is therefore taken into account. 271 There is no need for such a correction for R rs , which by defini-272 tion [Eq. (8)] already includes this normalization. 273 The photosynthetically available radiation (PAR) was also 274 measured above the surface (4.5 m a.s.l.) by a Satlantic PAR 275 sensor (400-700 nm, µmol quanta m -2 s -1 ). 276 B. Inversion Algorithms 277 Three algorithms have been used to derive IOPs from numerical 278 inversion of AOPs, in an attempt to evaluate the uncertainties 279 related to algorithm formulation. The first one was developed 280 by Morel et al. [24] (MOR in the following) and allows deriving 281 total absorption, a , and total backscattering, b b , from R and K d .282

  and measured b bp was verified through radiative transfer computations, before AOPs are introduced into inversion algorithms.In this optical closure exercise, radiative transfer simulations are parameterized with IOPs that follow the observed diel patterns as measured in situ at BOUSSOLE for two wavelengths (443 and 555 nm). The values of R, K d , and R rs obtained from these simulations can then be compared with the observations. The HydroLight code[START_REF] Gordon | Computed relationships between inherent and apparent optical-properties of a flat homogeneous ocean[END_REF] has been used, which solves the radiative transfer equation through the invariant imbedding technique and includes Raman scattering. Input parameters to HydroLight were set as follows.

( 4 )

 4 It was kept constant whatever the time of the day given that only one [Chl] value per day was available, implicitly assuming no diel changes in the non-water absorption. The high-frequency chlorophyll fluorescence measurements collected at BOUSSOLE cannot be used to infer possible 350 daily changes in [Chl] because of non-photochemical 351 quenching of fluorescence.352(3) For particulate scattering, the measured particle attenua-353 tion coefficient at 650 nm, c p (650), is used as a proxy for 354 b p (650). A spectral dependence is applied to get b p at 443 355 and 555 nm, following Eq. 14 in[START_REF] Morel | Bio-optical properties of oceanic waters: a reappraisal[END_REF], where the spectral 356 dependence is expressed as a power law whose exponent 357 varies from -1 to 0 when [Chl] increases from ∼0.02 to 358 ∼2 mg m -3 . 359 The scattering phase function of the particulate matter 360 was derived from the [Chl]-dependent parameterization 361 in [36]. This parameterization was run for each time step 362 using the [Chl] value that led to a VSF with the same bbp (λ) 363 as the measured value. By doing so, the diel cycle of b bp can 364 be taken into account. 365 (5) Boundary conditions included the Sun zenith angle (com-366 puted as a function of the day, time, and latitude), the sky 367 radiance conditions (clear-sky model [49] with an aerosol 368 optical thickness set to a fixed value of 0.2 at 550 nm), and 369 the sea state (wind speed of 5 m s -1 ). 370 D. Selected BOUSSOLE Datasets 371 Four situations with increasing total chlorophyll a concentra-372 tions (0.13, 0.26, 0.84, and 2.23 mg m -3 ) have been selected to 373 show how b bp , R, R rs , and K d evolve along the course of a day. 374 The two cases with the lowest concentrations are from 375 summer (July-August), and the other two from the spring 376 phytoplankton bloom (April). This selection follows the expec-377 tation that the magnitude of b bp will increase, and the shape of 378 its diel variability will change when waters become greener [18]. 379 The four cases are for clear skies.

380A

  larger dataset (N = 1290) has also been used to perform 381 the closure experiment and then to quantify the performance 382 of the inversions. This larger dataset includes clear-sky data 383 collected from 9 am to 5 pm local time at all seasons and for years 384 2006 to 2015. The selection criteria essentially include stable 385 physical conditions, as in [9].

  . Examples of Diel Variations of b bp , R, R rs , and K d The concurrent diel changes of the measured b bp , R, R rs , and K d for the four selected cases are shown in Fig. 1 for λ = 443 nm and 555 nm. The b bp diel cycles (Figs. 1A-1D) are characterized by a general increase as the day progresses, with minimal b bp around sunrise and maximum observed 3-6 h before sunset. The relative daily changes with respect to the value at sunrise observed in b bp ( ˜ b bp ) vary between 15% and 60%. These observations are consistent with [18]. Diel variability is also observed for R and K d at 555 nm and 443 nm whatever the trophic situation (Figs. 1E-1T), with a decrease from a maximum at sunrise to a minimum at solar noon and an increase toward another maximum at sunset. These diel variations are actually largely driven by changes in solar elevation. Little diel variations are observed for K d ,n , which can be considered a proxy of absorption. This observation supports the hypothesis we made about a constant non-water absorption coefficient in the radiative transfer calculations. Other studies (e.g., Schallenberg et al. [50]) have also shown an absence of diel variations in [Chl] in other environments.

[

  Chl]. The K d are generally below 15% (Figs. 2I-2P). The K d ,n are much smaller and do not display a significant diel cycle. B. Optical Closure As mentioned previously, optical closure was performed on a larger dataset in order to verify the consistency of the measured AOPs and b bp . The comparison indicates good agreement (Fig. 3, Table 2) confirming consistency of the dataset. Furthermore, the diel changes of the modeled R and K d compare well with those of the field determinations (red versus black curves in Fig. 1), except perhaps when the assumption of no diel changes in [Chl] is invalid, i.e., for the highest [Chl] values (Figs. 2H, 2L, 2P, and 2T). This consistency among field measurements of b bp and AOPs indicates that it can be

Fig. 1 .

 1 Fig. 1. Examples of daily changes of b bp (A-D, sampling dates and daily course of PAR, orange line, are displayed too), R and R rs (E-L), and K d (M-T) at 443 and 555 nm. The continuous black curves are for the measured quantities (except R rs , dotted curve), the dashed black curves for the normalized values (K d ,n ), and the red curves for the R and K d values obtained from radiative transfer computations. 4/C

452

  Nevertheless, degraded performances appear for the QAA 453 and GSM algorithms in the low range for λ = 555 nm, with a 454 significant overestimation of b bp values < 0.001 m -1 (Figs. 4D 455 and 4F), whereas MOR still works reasonably well (Fig. 4B; 456 significant dispersion but lower bias). 457 These results also show that it is difficult to get accurate 458 estimations of b bp in the blue by using only blue bands for the 459 inversion (MOR, Fig. 4A), whereas better results are obtained 460 when using methods that constrain the b bp derivation using more spectral bands (QAA and GSM; Figs. 4C and 4E). The comparison between the b bp values obtained from inversion of the measured R and K d or R rs and their measured values for the four case studies is shown in Figs. 5A-5D for λ = 443 nm and Figs. 5E-5H for λ = 555 nm. Again, the best results are obtained at 555 nm. Although the retrieved values generally reproduce the observed daytime increase of b bp , significant differences are observed in the amplitude of b bp (e.g., Fig. 5G) and in the shape of its diel changes (e.g., Fig. 5A). The results for the same four examples are also displayed in relative amplitudes of the diel change in b bp from sunrise, ˜ b bp , as calculated either from direct measurements or from AOPs inversion (Figs. 6A-6D for λ = 443 nm and Figs. 6E-6H for λ = 555 nm). The retrieved are higher by a factor of 1.4 ± 0.5% in the green and of 1.7 ± 1.1% in the blue with respect to measured values.

Fig. 2 .

 2 Fig. 2. Examples of daily relative changes (%) of R (black line; A-D at 443 nm, E-H at 555 nm), R rs (gray line; A-D at 443 nm, E-H at 555 nm), and K d at 443 nm (I-L) and K d at 555 nm (M-P). The dashed black curves in panels I-P are for the normalized values (K d ,n ).

  485 looked at the comparison between the absolute diel changes 486 of b bp , [ b bp ], calculated either from modeled or measured 487 data at 443 and 555 nm (Figs. 7C

490 4 .

 4 DISCUSSION AND CONCLUSIONS 491 Our case study has shown that the amplitude of the b bp diel cycle 492 observed in situ is large enough to generate a measurable diel 493 variability in reflectance. Although this result could be antic-494 ipated because the reflectance is at first order proportional to 495 b b , it still had to be verified quantitatively. Reflectance is, how-496 ever, also inversely proportional to absorption. Diel increases 497 in absorption would have to parallel those in b bp to make their 498 ratio, and then the reflectance, little variable over the day. This 499 is not what we observe, suggesting that possible changes in 500 absorption are not large enough to counterbalance those in b bp 501 in shaping the diel variability of reflectance. If absorption would evolve in the opposite way to b bp , which is not a priori expected, the result would be to enhance the diel variability in reflectance. The closure calculations have confirmed that the diel changes in AOPs could be well reproduced by accounting only for diel changes in b bp , ignoring possible diel changes in absorption.

Fig. 3 .

 3 Fig. 3. R rs (top row), R (middle), and K d (bottom) values computed by HydroLight compared to those measured in situ at 443 nm (A, C, E) and at 555 nm (B, D, F).

  estimations of b bp at 443 nm using Eq. (16) in oligotrophic waters.It could be, however, hypothesized that these inaccuracies would not change when inverting AOPs taken at different times of a day, so that they only marginally affect the retrieval of the quantity of interest here, namely the daily change in b bp and ultimately the daily change in POC, POC. This hypothesis seems confirmed in the green for MOR, where ˜ b bp and POC are retrieved with good accuracy (Figs.7B and 8B) but, again, not in the blue, where a large dispersion is observed (Figs.7A and 8A). Therefore, using only the blue part of the spectrum is likely inappropriate to perform inversion of R and K d in view of deriving ˜ b bp . The QAA and GSM algorithms provide similarly good results in the green (Figs. 7B and 8D, 8F,

Fig. 4 .Fig. 5 .

 45 Fig. 4. b bp retrieved from inversion of AOPs as a function of the measured b bp , for λ = 555 nm (A, C, E) and λ = 443 nm (B, D, F). Results are shown for the three inversion algorithms, as indicated on each panel. The dashed line is the 1:1 line.

Fig. 6 .

 6 Fig. 6. As in Fig. 5 but for the relative change in b bp from sunrise.

Fig. 7 .

 7 Fig. 7. A, B, ˜ b bp calculated from b bp values obtained from inversion of AOPs using the three algorithms, as a function of the same quantity derived from the field measurements. A is for λ = 555 nm, B is for λ = 443 nm, and the color coding is indicated in panel D. C, D, as in A, B but for the absolute diel change of b bp . The dashed line is the 1:1 line. 4/C

  Fig. 8. Diel change of POC calculated from the diel change in b bp , where b bp was retrieved from inversion of AOPs, as a function of the same quantity derived from the field measurements. A, C, and E are for λ = 443 nm and the MOR, QAA, and GSM algorithms, respectively. B, D, and F are for λ = 555 nm and the same three algorithms. The dashed line indicates the 1:1 line.

Fig. 9 .

 9 Fig. 9. R rs changes from its morning value to when it reaches its daily maximum, for λ = 443 nm A and λ = 555 nm B. The hatched histogram is for the simulated data and the grayed histogram for the measurements.

Funding.

  Agence Nationale de la Recherche; European Space Agency; Centre National d'Etudes Spatiales; Centre National de la Recherche Scientifique; Institut national des sciences de l'Univers; Sorbonne Université, Institut de la Mer de Villefranche (IMEV); Conseil Régional Provence-Alpes-Côte d'Azur. 3 Acknowledgment. The authors are grateful to the BOUSSOLE technical staff for their work in performing buoy deployments and maintenance, monthly cruises, and laboratory analyses. Our thanks go in particular to Emilie Diamond and Melek Golbol for the monthly cruises, and to Joséphine Ras, Céline Dimier, and Mustapha Ouhssain for the phytoplankton pigment analyses (SAPIGH analytical platform of the Laboratoire d'Océanographie de Villefranche, CNRS, France). This study is a contribution to the ANR Bio-optics and Carbon Experiment (BIOCAREX) and BOUSSOLE projects.

Finally, we are

  grateful to the external reviewers for providing constructive comments to the paper. Disclosures. The authors declare no conflicts of interest. NASA Press Release written by Felicia Chou, https://www.nasa.gov/ 666 press-release/nasa-targets-coastal-ecosystems-with-new-space-Ocean-Color Coordinating Group (IOCCG), "Ocean-669 colour observations from a geostationary orbit," in Reports of the 670 International Ocean-Colour Coordinating Group IOCCG No12, D.

Table 1

 1

	386	provides environmental
	387	information for this dataset.
	388	E. Statistical Indicators
	389	In order to assess the performances of the inversion algorithms,
	390	the following metrics were calculated, in log10 space, between
	391	retrieved and measured optical properties: the square of the
	392	Pearson correlation coefficient (R 2 ), the root mean square error
	393	[RMSE, Eq. (18)], the BIAS [i.e., the mean error, Eq. (19)], and
	394	the mean absolute error [MAE, Eq. (20)]:

Table 1 .

 1 Seasonal Repartition of the Selected Dataset

		Winter	Spring	Summer	Fall
		(DJF)	(MAM)	(JJA)	(SON) Total
	Number of days	9	20	28	15	72
	Number of in situ	155	361	504	270	1290
	data points
	[Chl] range	0.1-1.6 0.1-3 0.1-0.6 0.1-0.8 0.1-3
	(mg m -3 )

Table 2 .

 2 Statistics

		N	R 2	BIAS	MAE	RMSE
	HydroLight R r s versus in situ R r s
	443	1290	0.78	1.16	1.22	1.27
	555	1290	0.79	1.04	1.10	1.13
	HydroLight R versus in situ R
	443	1290	0.72	1.12	1.19	1.24
	555	1290	0.79	0.97	1.13	1.16
	HydroLight K d versus in situ K d
	443	1290	0.71	1.19	01.20	1.26
	555	1290	0.81	1.09	1.16	1.20
	b bp retrieved versus in situ b bp
	MOR
	443	1290	0.37	1.35	1.45	1.53
	555	1290	0.67	0.86	1.25	1.33
	QAA
	443	967	0.67	1.33	1.34	1.41
	555	967	0.65	1.40	1.60	1.72
	GSM
	443	1003	0.49	0.88	1.34	1.32
	555	1003	0.48	1.36	1.41	1.53

for the Comparisons Presented in Figs. 3 and 4: Number of Points (N), Correlation Coefficient (R 2 ), BIAS, Mean Absolute Error (MAE), and Root Mean Square Error (RMSE)

Table 3

 3

	) and better estimates in the blue
	region (Figs. 7A, 7C and Figs. 8C, 8E, and Table 3) as compared
	to MOR.
	As far as the R -K d algorithm (i.e., MOR) is concerned,
	the conclusions can be tempered, however, by reminding that
	the algorithm was not operated in ideal conditions, i.e., a single

[Chl] value has been used for a given day to derive µ d and f from lookup tables, whereas the actual [Chl] value at a given time of the day should in principle have been used to perform

Table 3 .

 3 Number of Points (N), Correlation Coefficient (R 2 ), BIAS, Mean Absolute Error (MAE), and Root Mean Square Error (RMSE), Values Corresponding to Data Displayed in Figs. 7 and 8

	N	R 2	BIAS	MAE	RMSE
	[

b bp ] retrieved versus in situ [ b bp ]

	MOR
	443	1290	0.37	1.12	2.28	2.83
	555	1290	0.78	0.86	1.61	1.78
	QAA
	443	967	0.71	1.08	1.53	1.75
	555	967	0.76	1.17	1.44	1.64
	GSM
	443	1003	0.73	1.17	1.44	1.65
	555	1003	0.79	0.91	1.49	1.75
	P OC retrieved versus in situ P OC
	MOR
	443	1290	0.37	1.13	1.85	2.07
	555	1290	0.78	1.03	1.39	1.43
	QAA
	443	967	0.71	1.07	1.38	1.46
	555	967	0.76	0.92	1.22	1.31
	GSM
	443	1003	0.71	1.19	1.42	1.51
	555	1003	0.77	0.87	1.36	1.47

Vol. 61, No. 18 / 20 June 2022 / Applied Optics

Antoine, ed. (2012).

K. Ruddick, G. Neukermans, Q. Vanhellemont, and D. Jolivet,

"Challenges and opportunities for geostationary ocean colour remote

sensing of regional seas: a review of recent results," Remote Sens.

Environ. 146,63-76 (2014).

D. A. Siegel, T. D. Dickey, L. Washburn, M. K. Hamilton, and B. G.

Mitchell, "Optical determination of particulate abundance and pro-

duction variations in the oligotrophic ocean," Deep-Sea Res. 36,

211-222 (1989).

H. Claustre, A. Morel, M. Babin, C. Cailliau, D. Marie, J.-C. Marty,

and D. Vaulot, "Variability in particle attenuation and stimulated

fluorescence in the tropical and equatorial Pacific: scales, patterns

and some biogeochemical implications," J. Geophys. Res. 104,

3401-3422 (1999).

P. Gernez, D. Antoine, and Y. Huot, "Diel cycles of the particulate

beam attenuation coefficient under varying trophic conditions in

the northwestern Mediterranean Sea: observations and modeling,"

Limnol. Oceanogr. 56, 17-36 (2011).

H. Loisel, V. Vantrepotte, K. Norkvist, X. Mériaux, M. Kheireddine, J.

"Penetration of UV-visible solar radiation in the global oceans:

insights from ocean color remote sensing," J. Geophys. Res. Oceans

118, 4241-4255 (2013).

Data availability. The BOUSSOLE dataset is accessible at http:// www.obs-vlfr.fr/Boussole/html/boussole_data/login_form.php.


Jean Auriol

Federico Bribiesca Argomedo

Observer design for n + m linear hyperbolic ODE-PDE-ODE systems

In this paper, an observer design is presented for a system of n+m hetero-directional transport partial differential equations (PDEs) coupled on both boundaries of the domain to ordinary differential equations (ODEs). This class of systems can represent, for instance, actuator and load dynamics at the boundaries of a hyperbolic system. The results in this paper provide a constructive way to reconstruct the state of the ODEs and the PDEs using only available measurements on one of the two ODEs. This observer design completes existing fullstate feedback designs for this class of systems and enables, together with previous results, the construction of outputfeedback stabilizers. As a complementary result, this paper includes a simple (constructive) method to obtain a stable leftinverse of a specific transfer matrix required in the observer design, which can also be applied (after a transposition) to the computation of the analogous stable right-inverses required for the control design.

I. INTRODUCTION

Many recent publications have focused on the control of interconnections involving hyperbolic PDEs (particularly ODE-PDE-ODE systems where the ODEs correspond to actuator and load dynamics). Indeed, this class of systems can model, e.g., the propagation of torsional waves in drilling systems [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] or deepwater construction vessels [START_REF] Stensgaard | Subsea hardware installation from a FDPSO[END_REF] (where the top of a cable is attached to a crane on a vessel at the ocean surface and the bottom attached to equipment to be installed at the seafloor).

Interconnections of hyperbolic PDEs and ODEs are not a new problem. Most constructive control results in the literature for these systems are based on the backstepping approach. Indeed, in the seminal paper [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], a reinterpretation of the classical Finite Spectrum Assignment [START_REF] Manitius | Finite spectrum assignment problem for systems with delays[END_REF] was proposed, modeling ODEs with input delays as PDE-ODE interconnections. Subsequently, this result has enabled the design of observers, controllers, or parameter estimation methods for a plethora of interconnected systems: systems with varying delays [START_REF] Bekiaris-Liberis | Nonlinear control under nonconstant delays[END_REF], [START_REF] Bresch-Pietri | Commande robuste de systèmes à retard variable: Contributions théoriques et applications au contrôle moteur[END_REF], cascades of PDEs [START_REF] Auriol | Delay robust state feedback stabilization of an underactuated network of two interconnected PDE systems[END_REF]. Of particular note are results concerning cascaded interconnections of hyperbolic PDE-ODE systems, such as [START_REF] Aamo | Disturbance rejection in 2 x 2 linear hyperbolic systems[END_REF], [START_REF] Hasan | Boundary observer design for hyperbolic PDE-ODE cascade systems[END_REF], [START_REF] Zhou | Boundary stabilization of a coupled wave-ode system with internal anti-damping[END_REF]. More recently, constructive results on non-cascaded PDE-ODE systems have been proposed. In particular, a stabilizing state-feedback control law has been proposed in [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF], [START_REF] Wang | Control of a 2× 2 coupled linear hyperbolic system sandwiched between 2 ODEs[END_REF]. Recently, a backstepping approach has been proposed in [START_REF] Deutscher | A backstepping approach to output regulation for coupled linear wave-ode systems[END_REF] for the output regulation for coupled linear wave-ODE systems. Some results also concern other types of PDE-ODE interconnections, such as parabolic PDEs, e.g. [START_REF] Tang | State and output feedback boundary control for a coupled PDE-ODE system[END_REF], [START_REF] Auriol | A differential-delay estimator for thermoacoustic oscillations in a Rijke tube using indomain pressure measurements[END_REF], [START_REF] De Andrade | Backstepping stabilization of a linearized ODE-PDE rijke tube model[END_REF], however this class of systems is beyond the scope of this article.

Regarding ODE-hyperbolic PDE-ODE systems with full interconnections (not limited to cascades), a stabilizing observer-controller robust to delays has been proposed in the case of a scalar proximal ODE in [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF]. In [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ODE-PDE-ODE systems[END_REF], an outputfeedback controller is designed based on assumptions that guarantee the existence of a Byrnes-Isidori normal form for the proximal ODE, as well as a relative degree one condition. These restrictions are partially avoided in [START_REF] Bou Saba | Backstepping stabilization of 2× 2 linear hyperbolic pdes coupled with potentially unstable actuator and load dynamics[END_REF] for the case of a scalar PDE system as the proximal ODE is simply assumed to be minimum phase for the output that affects the PDE. In [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF], a strictly-proper control law was proposed using less restrictive assumptions on the structure of the ODE components. This approach was later extended in [START_REF] Wang | Delay-compensated control of sandwiched ODE-PDE-ODE hyperbolic systems for oil drilling and disaster relief[END_REF] to encompass a state-observer (thus allowing the design of an output-feedback controller). However, the PDE subsystem is still assumed to be scalar in this work. Some recent results have also been obtained for interconnected PDE systems with non-linear ODEs [START_REF] Irscheid | Observer design for 2× 2 linear hyperbolic PDEs that are bidirectionally coupled with nonlinear ODEs[END_REF].

This paper deals with observer design for systems of n+m linear first-order hyperbolic Partial Differential Equations coupled with Ordinary Differential Equations at both boundaries of a one-dimensional spatial domain. The available measurements are obtained from one of the ODEs. The proposed approach extends the methodology introduced in [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF] to n+m transport equations (instead of 1+1), and combines the backstepping methodology with time-delay approaches. The assumptions made on the structure of the ODE components do not require the system to be of relative degree 1, nor to be in a specific canonical form. As such, existing results in the literature do not cover this class of systems. Our approach relies on the backstepping methodology. Using an invertible integral transformation, we map the system into a simpler target system. Then, we design a state observer for this target system using a time-delay representation and an analysis in the Laplace domain. The proposed design relies on the stable dynamic left-inversion of the (measured) ODE dynamics. We present a simple constructive procedure to find such a stable left-inverse. This procedure requires the same condition on the transfer function (full-column rank for all s ∈ C + ) than presented in Theorem 2 of [START_REF] Moylan | Stable inversion of linear systems[END_REF], yet operates directly on the transfer function instead of requiring switching back and forth between a frequencydomain and time-domain representations. In practice, a statespace realisation of this inverse system will only be found after adding an adequate low-pass filter, as in the case of the control (e.g. [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF]). Such a filtering procedure also guarantees the existence of robustness margins.

The paper is organized as follows: in Section II, we present the problem under consideration, as well as the hypotheses required for the observer design. In Section III, we design the state-observer using the backstepping methodology. In Section IV we present a simple constructive procedure to construct the a stable left-inverse for the ODE dynamics. Some concluding remarks are given in Section V.

Notations We denote the state space

χ = R p × L 2 ([0, 1]; R) n+m × R q ,
where p, n, m, q are positive integers.

For

(X 0 , u, v, X 1 ) ∈ χ, we introduce the corresponding χ-norm ||(X 0 , u, v, X 1 )|| 2 χ = ||X 0 || 2 R p + ||u|| 2 L 2 + ||v|| 2 L 2 + ||X 1 || 2 R q .
The variable s denotes the Laplace variable. To simplify the notations, the Laplace transform of a function f (t) (provided it is well-defined) will be denoted f (s). The space C + corresponds to the complex right half plane: C + = {s ∈ C, Re(s) ≥ 0}. The notation I r stands for the r × r identity matrix (if the dimensions are not ambiguous, the subindex will be omitted).

II. PROBLEM STATEMENT

A. Presentation of the system

In this paper, we consider a n+m linear hetero-directional hyperbolic system coupled through its boundaries with linear ODEs:

Ẋ0 (t) = A 0 X 0 (t) + E 0 v(t, 0), (1)
u(t, 0) = C 0 X 0 (t) + Qv(t, 0), (2)
u t + Λ + u x = Σ ++ (x)u(t, x) + Σ +-(x)v(t, x), (3) v t -Λ -v x = Σ -+ (x)u(t, x) + Σ --(x)v(t, x), (4) v(t, 1) = Ru(t, 1) + C 1 X 1 (t), (5)
Ẋ1 (t) = A 1 X 1 (t) + E 1 u(t, 1), (6)
defined for a.e.

(t, x) ∈ [0, +∞] × [0, 1]. The state of the system is (X 0 (t), u(t, •), v(t, •), X 1 (t)) ∈ χ.
The initial condition is taken as ((X 0 ) 0 , u 0 , v 0 , (X 1 ) 0 ) ∈ χ and we consider weak solutions to (1)-( 6) [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. The system is wellposed in the sense of [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem A.6, page 254]. The matrices Λ + and Λ -are diagonal and represent the transport velocities. We have Λ + = diag (λ i ) and Λ -= diag (µ i ) and we assume that their coefficients satisfy

-µ m < • • • < -µ 1 < 0 < λ 1 < • • • < λ n .
The spatially-varying matrices Σ •• are continuous (each coefficient of the matrix is a continuous function). With no loss of generality, we assume that the matrices Σ ++ and Σ --have zero diagonal elements [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ODE-PDE-ODE systems[END_REF]. The different coupling matrices satisfy

A 0 ∈ R p×p , E 0 ∈ R p×m , C 0 ∈ R n×p , A 1 ∈ R q×q , E 1 ∈ R q×n , C 1 ∈ R m×q , R ∈ R m×n , Q ∈ R n×m .
As we consider an observer design problem in this paper, and since the system under consideration is linear, we do not consider any control input in the system. We measure a part of one of the ODE state (as it is the case for drilling applications [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF]), i.e. the measurement Y (t) is defined by

Y (t) . = CX 1 (t), (7)
where C ∈ R d×q , d ≥ n. Having d ≥ n means that the number of measurements is greater than the dimension of the PDE state u. Without this assumption, it may not be possible to reconstruct the value of the function u(t, 1). As no result currently exist for under-measured hyperbolic PDE systems, it is a reasonable assumption.

B. General objectives

The objective of this paper is to design a state observer for the system (1)-( 6) based on the available measurement Y (t) given by equation [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. Even if the velocities are considered constant here, the proposed methodology could be extended to the case of spatially-varying transport velocities (following the approach given in [START_REF] Hu | Boundary exponential stabilization of 1-dimensional inhomogeneous quasi-linear hyperbolic systems[END_REF]). Note that, due to the symetry of the system, measuring the X 0 state would not change the nature of the problem.

C. Assumptions

Similarly to what has been done in [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF], we make the following assumptions Assumption 1: The pairs (A 1 , C), (A 0 , C 0 ) are detectable (i.e. there exist L 0 ∈ R p×n and L 1 ∈ R q×d such that Ã1 . = A 1 + L 1 C and Ã0 . = A 0 + L 0 C 0 are Hurwitz). Assumption 1 is a classical requirement found in most of the papers dealing with ODE-PDE-ODE systems [START_REF] Gehring | A systematic backstepping design of tracking controllers for ODE-PDE-ODE systems with nonlinear actuator dynamics[END_REF], [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ODE-PDE-ODE systems[END_REF]. It is not overly conservative since without the detectability of (A 0 , C 0 ) it becomes impossible to reconstruct the X 0 -state independently of the PDE or interconnection structure. The detectability condition on (A 1 , C) allows us to for a simpler observer design (no modes of X 1 are reconstructed indirectly through the PDE). It results in a set of conditions that can be easily tested.

Assumption 2: For all s ∈ C + , the matrices

(A 1 , E 1 , C) satisfy rank ( sI -A 1 E 1 C 0 ) = q + n. (8)
This last assumption (that requires d ≥ n) serves multiples purposes. First, it implies that C and E 1 are not identically zero. It is necessary to observe the PDE and the X 1 states from the available measurements. Under Assumption 2, the function

P 1 (s) . = C(sI -Ã1 ) -1 E 1 (9)
does not have any zeros in C + that is common to all its components. Thus, the function P 1 (s) admits a left inverse (denoted P - 1 ) whose entries have no unstable poles (such a left inverse is not proper) [START_REF] Moylan | Stable inversion of linear systems[END_REF]. A possible choice is given by the Moore-Penrose left inverse P - 1 (s) = (P T 1 (s)P 1 (s)) -1 P T 1 (s) (which should be verified to be stable a posteriori). If this is not stable, then a more involved stable inversion procedure is needed, an example of which is given in Section IV. This assumption is directly used in the constructive design of a observer and can be tested in a simple way.

Assumption 3: The system defined for all i ∈ [1, n] by

z(t) = m k=1 n ℓ=1 Q ik R kℓ z(t - 1 µ k - 1 λ ℓ ), (10)
is exponentially stable. Assumption 3 constitutes a reasonable assumption since it prevents system (1)-( 6) from having an asymptotic chain of eigenvalues with non-negative real parts [START_REF] Hale | Introduction to functional differential equations[END_REF], [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF]. It has been shown in [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF] that having an open-loop transfer function with an infinite number of poles on the closed right half-plane implies no (delay-)robustness margins in closed-loop (i.e., the introduction of any arbitrarily small delay in the actuation will destabilize the closed-loop system). Therefore, Assumption 3 is slightly stronger than a necessary condition for delay-robust stabilization and guarantees exponential stability of system system (1)- [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF]. If the delays are rationally independent 1 , Assumption 3 is equivalent to the following condition [START_REF] Hale | Introduction to functional differential equations[END_REF]: sup

θ kℓ ∈[0,2π] n×m Sp ( m k=1 n ℓ=1 Q ik R kℓ exp(jθ kℓ )) < 1,
where Sp denotes the spectral radius and j is the imaginary unit. This condition is simplified if the delays are rationally dependent. Furthermore, since the spectral radius of a matrix is upper-bounded by any norm of the matrix, easy to compute sufficient conditions for this spectral radius condition to hold can be derived using different norms of the matrices involved at the cost of increased conservatism.

III. OBSERVER DESIGN

The objective of this section is to design a state-observer for the system (1)-( 6) using the available measurement [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. In what follows, we denote with a • superscript the estimated states and with a • superscript the error state, i.e. the difference between the real state and the observer state. The objective is that the estimated state converges to the real state or equivalently that the error state, converges to 0, in the sense of the χ-norm. To simplify the computations and the design of the observer, we will work with a target system obtained from ( 1)-( 6) using a backstepping transformation, for which the in-domain coupling terms Σ •• have been moved to the boundary. Then, we design a state-observer that is a copy of this dynamics with some output injection operators.

A. Backstepping transformation

Inspired by [START_REF] Hu | Boundary exponential stabilization of 1-dimensional inhomogeneous quasi-linear hyperbolic systems[END_REF], we use an integral transformation to move the local coupling terms Σ •• to the boundary (in the form of integral terms). Consider the Volterra transformation T , similar to the one introduced in [START_REF] Hu | Boundary exponential stabilization of 1-dimensional inhomogeneous quasi-linear hyperbolic systems[END_REF], [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF]

X 0 (t) = ξ(t) - 1 0 L 1 (y)α(y) + L 2 (y)β(y)dy,(11) u(t, x) = α(t, x) - 1 x L αα (x, y)α(y)dy - 1 x L αβ (x, y)β(y)dy + γ α (x)X 1 (t), (12) v(t, x) = β(t, x) - 1 x L βα (x, y)α(y)dy, - 1 x L ββ (x, y)β(y)dy + γ β (x)X 1 (t), (13) X 1 (t) = X 1 (t), (14)
where the kernels are bounded functions defined either on

T u = {(x, y) ∈ [0, 1] 2 , x ≤ y}, or [0, 1]. This transformation rewrites (X 0 , u, v, X 1 ) = T (ξ, α, β, X 1 ). Denoting Λ = diag(Λ + , -Λ-), Σ = Σ ++ Σ +- Σ -+ Σ -- and L = L αα L αβ L βα L ββ , γ = (γ α , γ β ), we obtain ΛL x + L y Λ = Σ(x)L, Λγ x (x) = Σ(x)γ -γ ⊤ A 1 (15) (L 1 (x)) x Λ + = A 0 L 1 (x) + E 0 L βα (0, x) -L 0 (L αα (0, x) -QL βα (0, x) -C 0 L 1 (x)), ( 16
) (L 2 (x)) x Λ -= -A 0 L 2 (x) -E 0 L ββ (0, x) + L 0 (L αβ (0, x) -QL ββ (0, x) -C 0 L 2 (x)), (17)
with the boundary conditions [START_REF] Hasan | Boundary observer design for hyperbolic PDE-ODE cascade systems[END_REF] and γ α (1) = 0, γ β (1) = C 1 . Finally, we define L αα ij (0, y) for i ≤ j by

ΛL(x, x) -L(x, x)Λ = Σ(x), (18)
L 1 (0) = L 0 (Λ + ) -1 , L 2 (0)Λ -= L 1 (0)Λ + Q + E 0 ,
L αα (0, y) = QL βα (0, y) + C 0 L 1 (y). (20)
To this set of equations, we add arbitrary values for L αα ij (x, 1) (when i > j) and L ββ ij (x, 1) (when i > j) and L ββ ij (0, y) (when i ≤ j). Reinterpreting the ODEs in ( 16)- [START_REF] Gehring | A systematic backstepping design of tracking controllers for ODE-PDE-ODE systems with nonlinear actuator dynamics[END_REF] as PDEs evolving in the triangular domain T u with horizontal characteristic lines (since there is only an evolution along the x axis), it is possible to adjust the results from [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF]Theorem 3.2] to guarantee that the set of PDEs and ODEs ( 15)-( 20) has a unique solution which is piecewise continuous.The boundedness of transformation ( 12)-( 14) is a direct consequence of the structure of the transform (identities, integral operator and matrices) and the regularity of the different kernels. Its invertibility is a consequence of the structure of the transformation which is block triangular with the blocks on the diagonal being either identities (for the ODEs) or invertible Volterra operators (for the PDEs).The invertible backstepping transformation ( 12)-( 14) maps the original system (1)-( 6) to the following target system

ξ(t) = Ã0 ξ(t) + G 3 α(t, 1) + G 4 X 1 (t), (21)
α(t, 0) = Qβ(t, 0) + C 0 ξ(t) + (Qγ β (0) -γ α (0))X 1 (t) + 1 0 F α (y)α(t, y) + F β (y)β(t, y)dy, (22)
α t (t, x) + Λ + α x (t, x) = G 1 (x)α(t, 1), ( 23
) β t (t, x) -Λ -β x (t, x) = G 2 (x)α(t, 1), (24)
β(t, 1) = Rα(t, 1), Ẋ1 (t) = A 1 X 1 (t) + E 1 α(t, 1). ( 25
)
The functions G 1 and G 2 satisfy

G 1 (x) = 1 x L αα (x, y)G 1 (y) + L αβ (x, y)G 2 (y)dy -L αα (x, 1)Λ + + L αβ (x, 1)Λ -R -γ α (x)E 1 , (26)
G 2 (x) = 1 x L βα (x, y)G 1 (y) + L ββ (x, y)G 2 (y)dy -L βα (x, 1)Λ + + L ββ (x, 1)Λ -R -γ β (x)E 1 . (27)
The set of equations ( 26)-( 27) has a unique solution (Volterra equations of the second kind [START_REF] Yoshida | Lectures on differential and integral equations[END_REF]). The matrices G 3 and G 4 are defined by

G 3 (x) = L 2 (x, 1)Λ -R -L 1 (x, 1)Λ + + 1 0 (L 1 (x)G 1 (x)+L 2 (x)G 2 (x))dx, and G 4 (x) = E 0 γ β (0)+ L 0 (Qγ β (0) -γ α (0)).
Finally, the matrix F β the matrix F α are defined by F α (y) = L αα (0, y)-QL βα (0, y)-C 0 L 1 (y), and F β (y) = L αβ (0, y) -QL ββ (0, y) -C 0 L 2 (y). Note that F α is strictly lower triangular due to equation [START_REF] Hu | Boundary exponential stabilization of 1-dimensional inhomogeneous quasi-linear hyperbolic systems[END_REF]. The measurement Y (t) remains unchanged.

B. Observer equations

We can now design an observer for the target system ( 21)- [START_REF] Moylan | Stable inversion of linear systems[END_REF]. The observer state ( ξ, α, β, X1 ) is the solution of a set of equations that is a copy of the original dynamics to which we add dynamical output injection gains. We denote Ỹ (t) = Y (t) -C X1 (t), the difference between the real output and the observer output. The observer equations read as

ξ(t) = Ã0 ξ(t) + G 3 α(t, 1) + G 4 X1 (t) -O 0 ( Ỹ ), (28)
α(t, 0) = Q β(t, 0) + C 0 ξ(t) + (Qγ β (0) -γ α (0)) X1 (t) + 1 0 F α (y)α(t, y) + F β (y) β(t, y)dy -O 1 ( Ỹ ), (29)
αt (t, x) + Λ + αx (t, x) = G 1 (x) α(t, 1) -O α (x, Ỹ ), ( 30
)
βt (t, x) -Λ - βx (t, x) = G 2 (x)α(t, 1) -O β (x, Ỹ ), (31)
β(t, 1) = R α(t, 1), Ẋ1 (t) = A 1 X1 (t) + E 1 α(t, 1)-L 1 C Ỹ , (32)
with any (arbitrary) initial conditions in χ. The stable operators O i still have to be defined. Substracting the observer dynamics to the real one, we obtain the error system

ξ(t) = Ã0 ξ(t) + G 3 α(t, 1) + G 4 X1 (t) + O 0 ( Ỹ ), (33)
α(t, 0) = C 0 ξ(t) + Q β(t, 0) + (Qγ β (0) -γ α (0)) X1 (t) + 1 0 F α (y)α(t, y) + F β (y) β(t, y)dy + O 1 ( Ỹ ), (34)
αt (t, x) + Λ + αx (t, x) = G 1 (x) α(t, 1) + O α (x, Ỹ ) (35) βt (t, x) -Λ - βx (t, x) = G 2 (x)α(t, 1) + O β (x, Ỹ ) (36) β(t, 1) = R α(t, 1), (37) Ẋ1 (t) = Ã1 X1 (t) + E 1 α(t, 1). ( 38
)
The objective is now to tune the different operators O i such that the error system exponentially converges to zero. To do so, it is sufficient to show the convergence of ξ, α(t, 1) and X1 to zero. More precisely, we have the following lemma Lemma 1: If ξ(t), α(t, 1) and X1 (t) exponentially converge to zero, then the state ( ξ, α, β, X1 ) converges to zero in the sense of the χ-norm. This implies the convergence of the observer state to the real state.

Proof: Due to the stability of the observer operators and using the transport structure of ( 35) and (36), the exponential convergence of X1 and α(t, 1) to zero imply the exponential convergence of the states α(t, x) and β(t, x).

C. Design of the operators O i

We now want to define the operators O i such that ξ, α(t, 1) and X1 exponentially converge to zero. The analysis will be done in the Laplace domain. The Laplace transform 2 of equation ( 38) yields

(sI -Ã1 ) X1 (s) = E 1 α(s, 1). ( 39
)
Due to Assumption 1, the matrix (sI -Ã1 ) is invertible on C + . This implies Ỹ (s) = C(sI -Ã1 ) -1 E 1 α(s, 1). Thus, we obtain α(s, 1) = P - 1 (s) Ỹ (s), where P - 1 is a left-inverse of P 1 in ( 9). This in turns implies X1 (s) = (sI -Ã1 ) -1 E 1 P - 1 (s) Ỹ (s). This means that the terms that are functions X1 and α(s, 1) that appear in the error system can directly be compensated using the observer gains. In particular, we can define O 0 as

O 0 ( Ỹ (s)) = -(G 3 P - 1 (s) + G 4 (sI- Ã1 ) -1 E 1 P - 1 (s)) Ỹ (s) (40)
so that equation (33) can be rewritten as (sI -Ã0 ) ξ(s) = 0, which implies the exponential convergence of ξ to zero due to Assumption 1. Similarly, to get rid of the terms G 1 and G 2 , we define the operators O α (x, Ỹ ) and O β (x, Ỹ ) by 2 We omit the effect of the initial condition when taking the Laplace transform since it does not modify the stability analysis [START_REF] Hale | Introduction to functional differential equations[END_REF] such that equations ( 35)-(36) rewrite as transport equations. Indeed, for t > 1 λ1 + 1 µ1 , for every 1 ≤ i ≤ n and every 1 ≤ j ≤ m we now have for every

O α (x, Ỹ ) = -G 1 (x)P - 1 (s) Ỹ (s) (41) O β (x, Ỹ ) = -G 2 (x)P - 1 (s) Ỹ (s), (42)
x ∈ [0, 1] αi (t, x) = αi (t - x λ i , 0) (43) βj (t, x) = n k=1 R jk αk (t - 1 -x µ j , 1). ( 44
)
In what follows, we consider that t > 1 λ1 + 1 µ1 . The design of the operator O 1 is more involved since this operator must compensate almost all the terms that appear in equation (34) (including the integral terms). To design this observer operator, we will omit the term C 0 ξ that appear in equation ( 34) since this term exponentially converges to zero. For all 1 ≤ i ≤ n, we aim to obtain

αi (t, 0) = m k=1 n ℓ=1 Q ik R kℓ αℓ (t - 1 µ k - 1 λ ℓ , 0), (45)
so that α(•, 0) will exponentially converge to zero in virtue of Assumption 3 (see [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] for details). Since the matrix F α is strictly lower triangular, we will recursively compute the different components of O 1 . Consider the first line of equation ( 34) expressed in the Laplace domain. Using equations ( 43)-( 44), we obtain

α1 (s, 0) = ((Qγ β (0) -γ α (0)) X1 ) 1 + (O 1 ( Ỹ )) 1 + m k=1 n ℓ=1 Q 1k R kℓ e -s µ k -s λ ℓ αℓ (s, 0) + 1 0 m k=1 n ℓ=1 F β 1k (ν)R kℓ e - s(1-ν) µ k e -s λ ℓ αℓ (s, 0)dν, which gives (O 1 ( Ỹ )) 1 = -((Qγ β (0) -γ α (0))(sI -Ã1 ) -1 E 1 P - 1 (s) Ỹ (s)) 1 - m k=1 n ℓ=1 1 0 F β 1k (ν)R kℓ e - s(1-ν) µ k e -s
λ ℓ dν(P - 1 (s) Ỹ (s)) ℓ . By induction, let us consider i > 1 and assume that for any j < i, we have managed to design the j th component of our observer operator O 1 such that equation (45) holds. Consequently, using the triangular structure or F α , we obtain for any y ∈ [0, 1]

(F α (y) α(t, y)) i = i-1 j=1 F α ij (y) αj (t, y) (46)
= i-1 j=1 F α ij (y) m k=1 n ℓ=1 Q jk R kℓ αℓ (t - 1 µ k - y λ j , 1),
where the last equality holds since j < i. Consequently, the i th line of equation (34) now reads as

αi (s, 0) = ((Qγ β (0) -γ α (0)) X1 ) i + (O 1 ( Ỹ )) i + m k=1 n ℓ=1 Q ik R kℓ e -s µ k -s λ ℓ αℓ (s, 0) + 1 0 m k=1 n ℓ=1 F β ik (ν)R kℓ e - s(1-ν) µ k αℓ (s, 1)dν + 1 0 i j=1 F α ij (ν) m k=1 n ℓ=1 Q jk R kℓ e -sν
λ j e -s µ k αℓ (s, 1)dν.

To reach the desired target (as given by equation ( 45)) we choose

(O 1 ( Ỹ )) i = -((Qγ β (0) -γ α (0))(sI -Ã1 ) -1 E 1 P - 1 (s) Ỹ (s)) i - i j=1 m k=1 n ℓ=1 1 0 F α ij (ν)Q jk R kℓ e -sν λ j e -s µ k dν (P - 1 (s) Ỹ (s)) ℓ - m k=1 n ℓ=1 1 0 F β 1k (ν)R kℓ e - s(1-ν) µ k dν • (P - 1 (s) Ỹ (s)) ℓ . (47)
We can now write the following theorem Theorem 1: Consider the operators O 0 , O α , O β , O 1 , respectively defined by equations ( 40), ( 41), ( 42) and (47). Define the observer states ( X0 , û, v, X1 ) = T ( ξ, α, β, X1 ), where the transformation T is defined by equation ( 12)-( 14) and where ( ξ, α, β, X1 ) is the solution of the system ( 21)- [START_REF] Moylan | Stable inversion of linear systems[END_REF]. Then the state ( X0 , û, v, X1 ) exponentially converges to (X 0 , u, v, X 1 ) in the sense of the χ-norm.

Proof: With this choice of operators, we have already shown that X1 and ξ exponentially converges to zero. Moreover, the recursive design of O 1 implies that for all

t > 1 λ1 + 1 µ1 , all i ≤ n, αi (t, 0) is solution of αi (t, 0) = m k=1 n ℓ=1 Q ik R kℓ αℓ (t -1 µ k -1 λ ℓ , 0) + O( ξ),
where O is a linear bounded operator. Thus, α(t, 0) exponentially converges to zero (Assumption 3), which in turns imply the exponential convergence of α(t, 1) due to equation (43). Consequently, Lemma 1 implies that the state ( ξ0 , α, β, X1 ) exponentially converges to zero for the χ-norm. Using the invertibility and boundedness of the linear transformation T , we can easily conclude the proof. Using the linearity and invertibility of the different backstepping transformations, it is possible to express the observer system (28)- [START_REF] Zhou | Boundary stabilization of a coupled wave-ode system with internal anti-damping[END_REF] in the original coordinates. This is omitted due to space restrictions.

Remark 1: In a similar way to what has been done in [START_REF] Redaud | Practical output regulation and tracking for linear ODE-hyperbolic PDE-ODE systems[END_REF]Chapter 7] for the control and observer design in the 2 + 2 case, it is possible (thanks to Assumption 3) to low-pass filter the measured output signal Y (t) to use only strictly proper observer operators. Indeed, without filtering, the observer operators O i may not be strictly proper (due to the use of the left inverses) and the observer system may consequently be sensitive to delays in the measurements.

Remark 2:

The proposed observer could be combined with existing state-feedback laws [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ODE-PDE-ODE systems[END_REF], [START_REF] Gehring | A systematic backstepping design of tracking controllers for ODE-PDE-ODE systems with nonlinear actuator dynamics[END_REF] to obtain outputfeedback controllers. Note that we are currently extending the results of [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF] to deal with the case of non-scalar hyperbolic systems using assumptions that are less restrictive than those existing in the literature.

IV. STABLE LEFT-INVERSION ALGORITHM

In the construction of the observer, a stable left-inverse is required for the ODE system. In this section, we present a simple procedure to construct a stable left-inverse for P 1 (s). This procedure is given in order to give a complete, self-contained method in this paper and also to present an alternative, in many cases simpler, to the computation of a Hermite normal form of a matrix.

A. Preliminary Definitions

Recall that Assumption 2 guarantees that P 1 (s) is fullcolumn rank for all s ∈ C + . Furthermore, taking det(sI q -Ã1 ) as a common denominator (of degree q, with all its roots in the complex open left half-plane) we can factor P 1 (s) as P 1 (s) . = 1 det(sIq-Ã1) P num 1 (s), where P num 1 (s) has real polynomial entries of degree at most q-1 (it is a d×n matrix over a Principal Ideal Domain, which is in fact an Euclidean Domain, see e.g. [START_REF] Adkins | Algebra: An approach via module theory[END_REF]). The full-column rank property for s ∈ C + also applies to P num 1 (s). Given a list of real polynomials P = (l 1 (s), l 2 (s), ...l j (s)), we will denote by gcd(l 1 , l 2 , ...l j ) the polynomial greatest common divisor of the elements of P, and by (a 1 (s), a 2 (s), ..., a j (s)) .

= bezout(l 1 , l 2 , ...l j ) a corresponding list of real polynomial coefficients such that j k=1 a k (s)l k (s) = gcd(l 1 , l 2 , ...l j ). In order to construct a left-inverse for P num 1 (s) we will first transform it into an upper-triangular form. One possible upper-triangular form is the Hermite normal form, see for instance Theorem 2.9, in [START_REF] Adkins | Algebra: An approach via module theory[END_REF]Ch. 5]. However, in practice we do not require the uniqueness provided by this normal form, and it might be simpler to find a different uppertriangular form with the right properties. We provide a simple method that allows for one such construction. We begin defining some matrices that we will use to operate on the rows of P num 1 (s) to construct the desired upper-triangular form. The first matrix, a d × d upper-triangular matrix with real polynomial entries (of degree at most q-1) will allow us to replace the i-th row of a matrix by a combination of that row and the following ones and allow us to place a "pivot" element in the diagonal of the transformed matrix:

T i p [ci, ..., c d ](s) . = Ii-1 0 i-1,d-i+1 0 d-i+1,i-1 U [ci, ..., c d ](s) (48)
where U [c i , ..., c d ](s) is any d -i + 1 × d -i + 1 polynomial matrix with full rank for s ∈ C + and having as first row the polynomials [c i , ..., c d ]. Note that a particular choice for this matrix would be the unimodular (invertible) matrix used to construct the Hermite normal form (see for instance the matrix U 1 in the inductive proof of Theorem 2.9 in [3, Ch. 5]), it is also worth mentioning that, unlike the construction of the Hermite normal form, we do not require the elements above the diagonal to belong to a set of residues modulo the element on the diagonal. We believe that for the application considered in this paper, this formulation simplifies the necessary computations since, in many cases, one can simply complete the first line with d -i adequately chosen rows of the identity I d-i+1 , as long as one avoids rank deficiencies in C + (trivial if at least one of the c k polynomials has no roots in C + ). A second, d × d lower-triangular matrix with with real polynomial entries, T i l [p i+1 , ..., p d ](s), will allow us eliminate the elements under the previously constructed "pivot". It is constructed by replacing all the elements below the diagonal on the i'th column of a d × d identity matrix, by the column of polynomials [p i+1 , ..., p d ] T .

B. Construction of a stable left-inverse

A stable left-inverse of P 1 (s) can then be found by a method similar to Gaussian elimination detailed in Algorithm 1. Let us remark that a completely analogous algorithm can be used to find a stable right-inverse in the control case (simply acting on the columns of the transfer matrix, instead of the rows, or transposing the system).

If the reader already knows P (s), the Hermite normal form of P num 1 (s) and associated unimodular matrix T (s) such that P (s) = T (s)P num 1 (s), they can skip directly to step 11 of the algorithm. Remark that we do not require the classical condition of the Hermite normal form of having the elements above the diagonal (in the full column rank case) belonging to a complete set of residues modulo the elements of the diagonal (see [START_REF] Adkins | Algebra: An approach via module theory[END_REF]Ch. 5]), which also simplifies the procedure. (p i+1 , ..., p d ) = -P i+1:d,i (s)/P i,i (s) T (s) ← T i l [p i+1 , ..., p d ](s)T (s) 10: end for ▷ At the end of this loop, we obtain an upper triangular polynomial matrix P (s) with Hurwitz polynomials on the diagonal and zeros below the diagonal. 11: P (s) ← [I n 0 n,d-n ] P (s)

▷ We extract the first n rows of the matrix P (s), which are full rank in C + . 12: T (s) ← P -1 (s) [I n 0 n,d-n ] T (s)

▷ P , at this step a square, triangular matrix with Hurwitz polynomial entries in the diagonal, has a trivial stable inverse, and T (s) is therefore a stable, left inverse of P num 1 (s) 13: P - 1 (s) ← det(sI q -Ã1 )T (s) ▷ P - 1 (s) now contains a stable, left-inverse of P 1 (s)

V. CONCLUDING REMARKS

In this paper, we have designed a dynamic observer for a class of n + m linear hyperbolic ODE-PDE-ODE systems. The proposed approach combines backstepping transformations and frequency-domain design methods. The resulting observer requires the computation of a stable left-inverse, for which we have proposed a simple constructive procedure that can be easier to construct than a Hermite normal form. This observer can be combined with existing state-feedback control laws to obtain output-feedback controllers. In future works we will consider networks of ODE and PDEs with a more complex structure (star-shaped networks).

Algorithm 1 1 :

 1 P (s) ← P num 1 (s) 2: T (s) ← I d 3: for i = [1, 2, ..., n] do 4: (c i , ..., c d ) ← bezout(P i:d,i (s)) 5: P (s) ← T i p [c i , ..., c d ](s)P (s) 6:T (s) ← T i p [c i , ..., c d ](s)T (s)7:

8 :

 8 P (s) ← T i l [p i+1 , ..., p d ](s)P (s)

	9:

Jean Auriol is with Université Paris-Saclay, CNRS, Centrale-Supélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France. jean.auriol@l2s.centralesupelec.fr,

[START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] Federico Bribiesca Argomedo is with Université de Lyon, INSA Lyon, CNRS Ampère, F-69621, Villeurbanne, France

Extending the variable z, it is always possible to rewrite the system in a situation where the delays are rationally independent[START_REF] Hale | Introduction to functional differential equations[END_REF] Chapter 9].