J Calvo-De La Rosa

P Pomarède

P Antonik

F Meraghni

D S Citrin

D Rontani

A Locquet

Determination of the process-induced microstructure of woven glass fabric reinforced polyamide 6.6/6 composite using terahertz pulsed imaging

Terahertz pulsed imaging, combined with spatial and temporal signal and image processing, is performed to visualize the woven fabric in the various plies of glass-fiber-reinforced polymer laminates and to determine quantitative parameters characterizing the microstructure such as fiber-bundle orientation and the distance between yarns. The results are in quantitative agreement with the microstructure features provided by the manufacturer. In addition, the employed signal processing shows an excellent capability to reveal the weave pattern from noisy C-scans -where a visual-based analysis can be problematic -and to assess the process-induced microstructure.

Introduction

The growing need to replace metals with strong and stiff lightweight materials with other desirable mechanical, thermal, and chemical properties has led to the development of a range of composite materials. Among these, woven glass-fiber-reinforced polymer (GFRP) laminates are currently employed in high-impact sectors such as aeronautics [START_REF] Soutis | Fibre Reinforced Composites Aircraft construction[END_REF][START_REF] Wang | Nondestructive imaging of hidden defects in aircraft sandwich composites using terahertz time-domain spectroscopy[END_REF][START_REF] Ospald | Aeronautics composite material inspection with a terahertz time-domain spectroscopy system to cite this version : HAL Id : hal-00942971 Aeronautics composite material inspection with a terahertz time-domain spectroscopy system[END_REF], automotive [START_REF] Becker | Fatigue testing of GFRP materials for the application in automotive leaf springs[END_REF][START_REF] Sinha | Strength and corrosion analysis in alloy steel and E-glass composite wear ring in automotive engine cooling water pump[END_REF], renewable energy [START_REF] Yang | Progress and trends in nondestructive testing and evaluation for wind turbine composite blade[END_REF][START_REF] Jensen | Wind turbine blade recycling : experiences , challenges and possibilities in a circular economy[END_REF], and civil engineering [START_REF] Zhang | Strain rate studies of pultruded glass fibre reinforced polymer material properties : a literature review[END_REF][START_REF] Fang | Connections and structural applications of fi bre reinforced polymer composites for civil infrastructure in aggressive environments[END_REF]. In addition to a favorable density/strength relationship and significant resistance to impact, GFRPs are also low-cost composites [START_REF] Khan | Machinability of glass fibre reinforced plastic (GFRP) composite using alumina-based ceramic cutting tools[END_REF], in contrast to, for example, carbon-fiber-reinforced polymers (CFRP).

GFRP's properties arise from the combination of the low density of the polymer matrix and the high strength provided by the glass fibers. There exist multiple arrangements depending on the fibers' characteristics (length or width, for instance) and how they are spatially distributed [START_REF] Soutis | Fibre Reinforced Composites Aircraft construction[END_REF][START_REF] Callister | Materials science and engineering[END_REF]. Glass fibers can be randomly oriented in the polymeric matrix, leading to a transversely isotropic material. Another approach consists of unidirectionally oriented fiber resulting in an improved mechanical response in the direction of the fibers. This anisotropy makes these materials useful for applications where the stress direction is constant and well known. There are also composites in which fiber bundles (yarns or tows) are woven into fabrics, mixed with a polymer in thin plies, and stacked together, forming a multilayered woven laminate [START_REF] Soutis | Fibre Reinforced Composites Aircraft construction[END_REF]. It improves mechanical properties ranging from orthotropic to quasi-isotropic depending on the chosen stacking sequence. Consequently, this type of GFRP is most commonly used for structural elements in demanding applications. Therefore, accurate knowledge of the composite's microstructure is vital to rely on its mechanical capabilities. This study focuses on such GFRP laminates incorporating woven glass-fiber fabrics in a thermoplastic matrix. As for all composite types, the physical characteristics of the fibers (diameter, length, orientation, concentration, bundle structure and weave pattern) play a crucial role in the determination of the mechanical behavior of the composite [START_REF] Zhang | Strain rate studies of pultruded glass fibre reinforced polymer material properties : a literature review[END_REF][START_REF] Callister | Materials science and engineering[END_REF].

When these materials are used for critical components (such as structural components in aircrafts), it becomes crucial to test them both in production and service [START_REF] Yang | Progress and trends in nondestructive testing and evaluation for wind turbine composite blade[END_REF][START_REF] Karbhari | Non-destructive evaluation (NDE) of polymer matrix composites[END_REF]. Regarding production, a proper quality-control process is essential to confirm that the produced GFRP laminate has the desired microstructure and that there are no significant defects before being deployed as a load-bearing component. For in-service testing, routine evaluation may be required to detect the presence of damage that may lead to the component's failure. There are many different nondestructive-evaluation (NDE) methods that have been tested or used with this purpose [START_REF] Karbhari | Non-destructive evaluation (NDE) of polymer matrix composites[END_REF][START_REF] Miqoi | Detection and evaluation of barely visible impact damage in woven glass fabric reinforced polyamide 6.6/6 composites using ultrasonic imaging , Xray tomography and optical profilometry[END_REF][START_REF] Eckel | Investigation of damage in composites using nondestructive nonlinear acoustic spectroscopy[END_REF]. The most common techniques are ultrasonic inspection, x-ray tomography, thermography, and microwave imaging. Terahertz (THz) pulsed imaging is a promising alternative that provides several attractive features but presents challenges. It has appropriate spatial resolution (due to its wavelength, λ ~ 300 μm in the air at 1 THz) to identify critical structural characteristics and defects. Moreover, since THz electromagnetic radiation penetrates polymers and glass, quantitative structural information at depth can potentially be obtained. In addition, THz radiation is nonionizing, as opposed to x-rays, and thus presents minimal health risks to the operator. Compared to standard approaches such as ultrasonic imaging, THz waves require no coupling fluid, and the materials are relatively transmitting, unlike ultrasonic signals at the high frequencies required to obtain quantitative three-dimensional information at depth.

THz-based techniques have already drawn attention in the last decade for the NDE of GFRP materials. Continuous-wave (CW) THz systems, which work with monochromatic electromagnetic radiation, have already provided interesting results on the damage detection and inspection of various polymeric and composite materials [START_REF] Cristofani | Nondestructive testing potential evaluation of a terahertz frequencymodulated continuous-wave imager for composite materials inspection[END_REF][START_REF] Dai | Improved terahertz nondestructive detection of debonds locating in layered structures based on wavelet transform[END_REF][START_REF] Zhang | Broadband stepped-frequency modulated continuous terahertz wave tomography for non-destructive inspection of polymer materials[END_REF], though usually not providing depth information. Several studies [START_REF] Wang | Nondestructive imaging of hidden defects in aircraft sandwich composites using terahertz time-domain spectroscopy[END_REF][START_REF] Ospald | Aeronautics composite material inspection with a terahertz time-domain spectroscopy system to cite this version : HAL Id : hal-00942971 Aeronautics composite material inspection with a terahertz time-domain spectroscopy system[END_REF][START_REF] Dong | Nondestructive evaluation of forced delamination in glass fiber-reinforced composites by terahertz and ultrasonic waves[END_REF][START_REF] Dong | Polarization-resolved terahertz imaging of intra-and inter-laminar damages in hybrid fi ber-reinforced composite laminate subject to low-velocity impact[END_REF][START_REF] Han | Nondestructive evaluation of GFRP composite including multidelamination using THz spectroscopy and imaging[END_REF] have also successfully used pulsed THz radiation to observe the presence of buried defects in GFRP laminates. For instance, Jördens and co-workers [START_REF] Jördens | Terahertz spectroscopy to study the orientation of glass fibres in reinforced plastics[END_REF] used THz time-domain spectroscopy (THz-TDS) to explore the birefringent properties of unidirectionally oriented GFRP. They showed how the refractive index and absorption coefficient depend on the relative orientation between the THz polarization and the glass fibers. Similar procedures have been later applied by other authors [START_REF] Albrecht | Measuring fibre orientation in sisal fibre-reinforced, injection moulded polypropylenepros and cons of the experimental methods to validate injection moulding simulation[END_REF]. In a recent paper [START_REF] Park | Predicting the stacking sequence of E-glass fiber reinforced polymer (GFRP) epoxy composite using terahertz time-domain spectroscopy (THz-TDS) system[END_REF], Park et al. used the same birefringence principle to predict the primary fibers' orientation and differentiate between three different stacking sequences of woven fibers. This THz-TDS approach leveraging off the fabric-induced birefringence has hitherto been restricted to determining the fibers' orientation, and not used to extract more detailed information concerning the woven fabrics. With a somewhat different type of material, Tang et al. [START_REF] Tang | Quantitative evaluation of fiber structure by using coherent terahertz wave[END_REF] measured the orientation and spacing between fibers in a silk-nylon woven sample based on Fraunhofer diffraction. Similar research has been conducted on CFRP, where defects and damage have been observed [START_REF] Dong | Polarization-resolved terahertz imaging of intra-and inter-laminar damages in hybrid fi ber-reinforced composite laminate subject to low-velocity impact[END_REF][START_REF] Wang | Study of the skin depth and defect detection in carbon fiber composites with Terahertz waves[END_REF][START_REF] Xu | Terahertz imaging and vibrothermography for impact response in carbon fiber reinforced plastics[END_REF][START_REF] Dong | Visualization of subsurface damage in woven carbon fiber-reinforced composites using polarization-sensitive terahertz imaging[END_REF]. However, the THz wave penetration was limited because of the conductivity of the carbon fibers; hence hindering analyses in-depth.

In this work, we observe woven fabrics and quantify their processinduced microstructure in three dimensions in GFRP laminates. Our approach is based on a THz raster scan in reflection by a concatenation of image and signal processing tools to extract specific information about the weave structure, such as the fiber orientation and the gap between two tows, in the various plies. Specifically, the approach provides a quantitative and global measure of the orientation of tows within the weave pattern and their period. These quantities are characterizable and can ascertain the origin of the material's manufacturing or its damage in service. The proposed approach could be easily implemented in standard C-scan inspection.

Materials and methods

Sample

The studied material is a woven twill 2/2 composite called Vizilon SB63G1-T1.5-S3, manufactured by DuPont de Nemours using a thermocompression molding process. It is a co-polyamide 6.6/6 (PA 6.6/6, thermoplastic matrix) reinforced with a layup of three balanced glass fiber plies oriented at ± 45 • (warp and weft yarns) as shown in Fig. 1. The composite, with a thickness of 1.53 mm, has a fiber weight ratio of 63%, corresponding to a fiber volume fraction of 43%. The sample has a surface dimension of 21 mm × 77 mm. A detailed structural characterization of the same material can be found in Ref. [START_REF] Pomarède | Damage evaluation in woven glass reinforced polyamide 6 . 6/6 composites using ultrasound phaseshift analysis and X-ray tomography[END_REF].

Experiment

We conduct THz imaging of the composite sample with a TPS Spectra 3000 (TeraView, Cambridge, UK) in reflection mode at room temperature (25 • C). The full sample surface was raster scanned with roughly single-cycle THz pulses with a transverse spatial step size of 0.4 mm for both the x and y directions (parallel to the sample's surface), leading to ~ 53 × 193 pixel images. These pulses impinge on the sample nearnormal incidence, and the reflected signal is detected close to the specular direction. Due to refractive-index discontinuities, reflected pulses (echoes) occur. The time delay at their detection allows the construction of the sample's depth profile at the position studied if we know the refractive indices. Repeating this process for all pixels, we can obtain a three-dimensional dataset representative of the sample structure. Finally, we produce a two-dimensional cut of this dataset, showing the reflected intensity at each pixel for a given time delay, which is commonly known as a C-scan.

The spatial resolution of the measurement has two components. The first one is the transverse resolution, which is close to the spatial step size (in our experiments, the step size is slightly larger than the beam spot size -about ~300 μm at 1 THz-on the surface). The second one is the axial resolution, which is linked to half the shortest wavelength present in the high signal-to-noise band within the medium, which in our case is ~46 μm.

Image-processing algorithm

In this subsection, we present our THz imaging workflow and methodology. It relies on the Fourier analysis of raw THz images (Cscans) and the searching of spectral signature (and their associated characteristics) of periodic structures. The results obtained at each step will be presented later -in Section 3.1 -in the same order in which they have been applied. All the image processing has been implemented in Python® and its libraries. Fig. 2 is a schematic representation of the various methods applied, from the raw THz image to the extraction of the desired parameters.

The first step is image averaging. In this process, we produce a new low-pass filtered C-scan where the intensity of each pixel is the average of (2n+1) 2 neighboring pixels (n pixels in each direction) in the raw image. The aim is to smoothen the image by reducing high-spatialfrequency noise. For a number n of neighbors, the averaged intensity I i,j of a pixel with spatial coordinates (i,j) is Fig. 1. Reinforcement architecture of the woven glass fabric-PA 6.6/6 composite showing the warp and weft yarns [START_REF] Miqoi | Detection and evaluation of barely visible impact damage in woven glass fabric reinforced polyamide 6.6/6 composites using ultrasonic imaging , Xray tomography and optical profilometry[END_REF].

I i,j = ∑ j+n y=j-n ∑ i+n x=i-n I x,y (2n + 1) 2 (1)
where I x,y represents the measured intensity at pixel (x,y). We correct for those pixels where i or j < n (i.e., at the edge of the image). The second step is binarization [START_REF] Russ | The image processing handbook[END_REF][START_REF] Gonzalez | Digital image processing[END_REF], which helps distinguish between different materials or objects in the image. We have used a global thresholding method in situations where the distinction between the two materials is evident. As a result, pixels with an intensity larger than a defined threshold T are set to one and zero otherwise. In more complicated situations (such as in fibers-matrix differentiation), the adaptive (also called dynamic) thresholding approach has been used [START_REF]scikit-image[END_REF]. In this case, the local threshold value T i,j is the weighted mean of the neighbors of pixel (i,j). The number of neighbors considered in each direction (i.e., the threshold radius) has been set to 3. Thus, in both binarization methods, each pixel has a final binary value B i,j obtained by applying the criterion

B i,j = { 1 if I i,j > T i,j 0 otherwise (2)
The third step relies on computing the two-dimensional discrete Fourier transform (2D-DFT), a commonly-used processing tool for filtering images and extracting features of interest [START_REF] Russ | The image processing handbook[END_REF][START_REF] Broughton | Discrete fourier analysis and wavelets[END_REF]. Here, we aim at extracting characteristics of periodic patterns (i.e., the weave structure in the composite); the 2D-DFT should exhibit pronounced peaks at spatial frequencies corresponding to these periodicities [START_REF] Marquez | Fourier analysis and automated measurement of cell and fiber angular orientation distributions[END_REF]. The 2D-DFT is given by

X(u, v) = ∑ M-1 m=0 ∑ N-1 n=0 x(m, n)e -j2π ( um M + vn N ) (3)
where x(m, n) is the intensity of pixel (m,n) in the spatial domain in a M x N-pixel image, X(u, v) is the complex-valued magnitude at the (u,v) coordinates in the 2D spectral domain, and

j = ̅̅̅̅̅̅ ̅ -1 √ . The radial (spatial) frequency u r is defined as u r = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ u 2 + v 2 √
, and the angle with respect the horizontal zero frequency axis θ as θ = tan -1 (v /u) [START_REF] Marquez | Fourier analysis and automated measurement of cell and fiber angular orientation distributions[END_REF]. In case that u = 0, θ = 90

• if v > 0 and θ = -90 • if v < 0.
Finally, in the fourth and final step, we implement a high-pass filter [START_REF] Russ | The image processing handbook[END_REF][START_REF] Gonzalez | Digital image processing[END_REF] to remove low-level spatial frequency noise in the obtained 2D spectra. As we apply this filter in just one of the directions (u or v) at each step, we set a wavenumber threshold (wT u/v ) below which all the intensities are set to zero for this specific direction. The criterion is

F u,v = { X u,v if u / v > wT u/v 0otherwise (4)
where F u,v is thus the filtered complex-valued magnitude at (u,v) coordinate in the spectral domain.

To summarize, in steps 1 and 2, the original spatial-domain image is adapted to reveal the matrix-fiber microstructure for a given time slice. Once both phases are differentiated, we proceed in step 3 to apply a 2D-DFT to obtain well-defined peaks in the spectral domain, associated with the mean woven structure. In step 4, filtering of the 2DFTT separates the periodic structure from the noise. The resulting 2D spectrum has two purposes: (i) Locate the maxima and calculate the corresponding distance and angle in the spatial-domain, and (ii) perform the inverse 2D-DFT and obtain a graphical representation of the mean fibers' arrangement in the sample. We have repeated this methodology with multiple C-scans obtained at different time delays to ensure the robustness of our analysis.

It is worthy of note that, in addition to the image processing techniques described here, signal processing of the temporal signal corresponding to each pixel is also applied, as described later, to account for the noise of the experimental THz pulse. Specific parameter values used in the various image and signal processing steps might need to be adapted to the sample under investigation, though, in the following, we guide the reader by discussing how these parameter values are determined.

Results and discussion

The measurement and the subsequent data analysis correspond to a basic THz time-of-flight tomography (TOFT) experiment [START_REF] Zhai | Nondestructive measurement of mill-scale thickness on steel by terahertz time-of-flight tomography[END_REF][START_REF] Zhai | Terahertz time-of-flight tomography beyond the axial resolution limit: autoregressive spectral estimation based on the modified covariance method[END_REF]. Depth information in THz TOFT is extracted from echoes in the reflected signal at a given spatial location, which originates from a discontinuity in the refractive index (such as those induced by the matrix-fiber-bundle interface). The depth information is in the relative time delays between the various echoes (accounting for the refractive indices of the materials). The process is then repeated at all transverse pixels across the sample area, leading to a block of three-dimensional data (two spatial dimensions x and y, one temporal dimension).

One way to represent this information is to take a 2D cut through the data (using one spatial dimension and the time dimension) to produce a B-scan, which provides a cross section through the sample, time serving as a surrogate for the axial coordinate. Fig. 3 shows the B-scan obtained at approximately the center of the sample along the short axis (y = 30 mm). To increase the contrast and facilitate the observation of the The air-sample interface stands out at 4 ps in Fig. 3 as a pronounced red band due to the large refractive-index mismatch between air and the polymer matrix. One observes successive maxima at increasing time delays within the sample, showing a tendency to form horizontal bands in the B-scan. These red bands correspond to abrupt changes in refractive index, as it may happen for interfaces such as air-polymer, polymerglass fibers or inter-ply interfaces. That is, echoes occur due to the stratified nature of the CFRP laminate. Due to the different refractive indices of the matrix and glass-fiber bundles, these bands reveal the woven structure at various depths. However, it is challenging to directly extract clear quantitative information on the material internal structure from this three-dimensional dataset.

To better extract quantitative information about the weave structure, it is more fruitful to consider C-scans, i.e., two-dimensional sections parallel to the surface within a narrow time-delay window. Below, in Sec. 3.1, we first present a complete analysis from the main echo (airsample interface) study, revealing the characteristics of the fabric in the first layer. Section 3.2 shows how the proposed algorithm can examine the fabric in the three deeper layers. Finally, in Sec. 3.3 the accuracy of the method is confirmed by using a larger set of artificially created noisy C-scans.

First echo analysis

The measured THz pulses have been initially processed using typical denoising procedures to overcome experimental artifacts. In our case, we first applied an ideal high-pass filter with a cutoff frequency of 0.2 THz (which removes low-frequency noise) and the resulting pulse was filtered with frequency-wavelet domain deconvolution to reduce noise [START_REF] Zhai | Thickness characterization of multi-layer coated steel by terahertz time-of-flight tomography[END_REF][START_REF] Zhai | Pulsed THz imaging for thickness characterization of plastic sheets[END_REF]. Fig. 4 shows an averaged (see above) C-scan at the sample's surface (time delay ~4 ps) obtained from these pulses. The optimal number of neighbors used in the averaging is determined to be 2 for these samples; this value helps remove image noise while not introducing significant blur that would complicate the distinction of the bundles.

This image reveals the presence of two perpendicular fiber bundles oriented at ±45 • (warp and weft), in good agreement with the manufacturer's description. We observe a magnified amplitude of the reflected signal at the center of each bundle. It results from the enhanced reflection of the THz pulse at this position under normal incidence at the top of their circular cross-section. Subsequently, we binarize this C-scan to identify the fibers better and differentiate them from the matrix. To do so, we have used an adaptive thresholding method, as illustrated in Fig. 5 (A). The definition of the borders between the fibers and the matrix is not well resolved. Assuming a continuous shape of the fibers, we have smoothened the binary C-scan with an image-averaging procedure. The result is shown in Fig. 5 (B).

The improvement in the definition of the fibers is substantial. Indeed, Fig. 5 (B) provides clear identification of the matrix (blue) and fibers (red).

Once the two materials are well-differentiated, we take the 2D-DFT of the image to quantify the periodicities present. The modulus square of the 2D-DFT of Fig. 5 (B) is shown in Fig. 6 (A).

The image exhibits four peaks surrounding the zero-frequency point (itself of no interest) in the two directions of the periodicities seen in the C-scans. Two peaks (corresponding to one of the two bundles' orientations) are stronger than the other two due to the more pronounced appearance of one of the directions in the C-scan (Fig. 4). This difference in magnitude of the two fiber orientations is attributed to birefringence. It is essential to highlight that the quality of the experimental measurements influences the contrast in this analysis. The peaks have a nonzero spectral extent due to (i) the finite sample size (much smaller in the horizontal than vertical directions, given greater horizontal extension of the sample) and (ii) departures from the ideal periodicity of the fiber bundles.

To optimize the extraction of information from the four prominent peaks, we have first applied a combination of a horizontal and a vertical high-pass filter to eliminate all these noisy periodic signals along the vertical and horizontal axes. By this process, we set all pixels' intensity to zero when they are below a defined wavenumber threshold for each axis. In this case, we have used a 4 × 10 -2 mm -1 threshold for both cases, considering the extension in the spatial frequency domain of the central noise in Fig. 6 (A). The result is shown in Fig. 6 (B).

The horizontal and vertical noisy patterns have been eliminated -as well as the central peak. Consequently, the four peaks representing the constant texture in the C-scans can be better differentiated, especially the two weaker ones. The downside of this process is the detection of lower-intensity parasitic signals in the background, which do not represent the main fibers' pattern. Nevertheless, given that the relative intensity of the noise is lower than that of the peaks, it can be eliminated by a simple thresholding process. Due to the difference in relative intensity between the two pairs of peaks, we define a threshold at 80% of the maximum intensity for each of the directions (i.e., at each of the four 2D-DFT quadrants) to ensure correct detection. Smaller threshold values might introduce background noise, while larger ones could lead to the loss of meaningful peaks. The result is shown in Fig. 6 (C), in which the centers of the peaks can be clearly identified. The spectral coordinates are u = -0.19 mm -1 and v = 0.17 mm -1 for the upper-left peak (the weakest, corresponding to the warp direction), while the coordinates are u = 0.18 mm -1 and v = 0.17 mm -1 for the upper-right peak (the most intense, associated to the weft direction). Due to the symmetric nature of the 2D-DFT we only consider the peaks in the upper half of the image. These coordinates are used to calculate the corresponding radial frequency and angle in the frequency domain, which result in u r = 0.26 mm -1 and θ = -41.61 • for the left-side peak, and u r = 0.25 mm -1 and θ = 44.52 • for the right-side peak. The corresponding space-domain distance is 3.86 mm for the warp direction and 4.05 mm for the weft direction. This distance corresponds to the distance between the ridges observed in the C-scans, i.e., the average distance between the top of the fiber bundles. Additionally, the quantitative measurement of the fiber orientations leads to the respective values of 41.61 • and -44.52 • with respect to the vertical axis (see Fig. 7 (C) for a schematic representation). We notice that the weave pattern does not form a perfect right angle. Therefore, with this method, we can also detect assembly irregularities between the different bundles. Overall, the values are in good agreement with the theoretical values provided by the manufacturer, the structural data given in Ref. [START_REF] Pomarède | Damage evaluation in woven glass reinforced polyamide 6 . 6/6 composites using ultrasound phaseshift analysis and X-ray tomography[END_REF], and the C-scan shown in Fig. 4.

To check the accuracy of our results, we have tested this algorithm with the surface C-scans of a set of four different samples made of the same material. We have obtained a mean distance value of 3.90 ± 0.03 mm, an orientation of 42.90 ± 1.82 • for the warp direction, a mean distance value of 4.08 ± 0.08 mm, and an orientation of -45.63 ± 0.92 • for the weft. It reflects the consistency of the extracted data.

Finally, we have used the processed 2D-DFT image, represented in Fig. 6 (C), to perform the inverse 2D-DFT and obtain a clear representation of the mean disposition of the fibers at this specific time delay. In order to allow easy comparison, we provide in Fig. 7 (A) the obtained inverse 2D-DFT together with the averaged C-scan in Fig. 7 (B).

After the image processing, we see that Fig. 7 (A) reveals the fiber structure suggested by the C-scan. As expected from the experimental Cscan, one of the two diagonals (weft) is more intense than the other (warp). In addition, the positions with higher intensity result from the intersection of the two diagonal patterns, leading to a stronger intensity in these points.

In conclusion, the proposed method does not only allow one to extract quantitative information about the woven fabrics, but it also

In-depth analysis

We analyze here multiple C-scans, obtained at various time delays, to extract information about the weave structure in subsurface plies of one laminate. To do so, we select the optical delays corresponding to the three interfaces previously identified (12.0, 19.5, and 27.5 ps), and an inner point of the first layer that has a strong reflection (6.0 ps). In Fig. 8, we present the B-scan with the selected time delays marked and corresponding C-scans.

From Fig. 8, we observe that there is still a preferential fiber-bundle direction, as on the surface, although both may be observed in all cases.

On the other hand, the weave patterns appear to be less clear than in the surface-ply C-scan, and it can even be difficult to identify them visually from the raw C-scans. It is the consequence of the increased signal attenuation at larger depths, which leads to an increased role of noise. Moreover, as one concentrates on longer time delays, shadow effects due to overlying features contribute artifacts to the C-scans. However, the contrast and quality of the C-scan increases again for the back surface of the sample (27.5 ps). It is due to a larger refractive index mismatch between the sample and air. We have applied the same algorithm, with the same parameters as in Section 3.2, to reveal the woven structure and extract quantitative information. The summarized results are presented in Fig. 9.

The background noise is, in all cases, slightly larger than in the 2D-DFT pattern obtained for the C-scan at 4 ps due to the lower signal to noise associated with subsurface plies. Again, the 2D high-pass filter of the 2D-DFT allows us to remove the main peak and noise and thus identify the peaks corresponding to the woven structure. These filtered squared modulus 2D-DFT images have been thresholded at 80% of the maximum intensity in each direction, and the inverse 2D-DFT has been computed to reconstruct the woven fabrics. As the reader may notice, despite the lower definition of the raw C-scans, this process can identify the fiber-bundle directions and reveal the fabric structure. Compared with the surface analysis, these results are of importance due to two aspects: first, we observe the woven fabric deeper than the surface ply thanks to the penetration of THz waves. Second, the proposed image and signal processing enhances the physical measurement and reveals the periodic fabrics from C-scans, where this information is not apparent to the eye in the raw C-scans.

Finally, by computing the coordinates of the modulus squared 2D-DFT, we deduce the fibers' orientation and spacing at different depths (see Table 1). The standard deviation for the orientations of the bundles over all the plies is 1.30 • and 1.48 • , respectively. It is in the expected order of magnitude, when compared with other studies dealing with fiber orientations dispersion through the thickness on similar materials [START_REF] Nguyen | Micro-CT measurement of fibre misalignment: application to carbon/epoxy laminates manufactured in autoclave and by vacuum assisted resin transfer moulding[END_REF][START_REF] Malcom | Three-dimensionally woven glass fiber composite struts : characterization and mechanical response in tension and compression[END_REF]. The consistency of these values-and their agreement with the manufacturer's values-demonstrates that the woven fabrics are the same and have similar orientations in each layer, which is reasonable as the sample has been analyzed as received, without suffering deformation or damage. At the same time, this becomes a powerful tool to measure and control small deviations from the ideal values expected from design at each ply and depth.

Method accuracysimulation analysis

After having proven that the approach works with real samples, we test its accuracy and robustness. To do so, we work with artificially produced C-scans, because we can generate a wider range of sample characteristics and, at the same time, we remove all the experimental error, limitations, or noise.

With this goal, using a Python® code, we have generated a set of binary C-scans with the two fiber-bundle directions (mimicking the sample studied in this work) but varying the two measurable parameters and creating all their possible combinations: (i) orientations from the 10 • to the 80 • from the vertical axis, in steps of 5 • ; (ii) distances from 2 mm to 6 mm, in steps of 0.5 mm. We have therefore tested 135 different binary C-scans. The fibers' width is set to 1 mm, the sample size to 60 mm in the y-direction and 20 mm in the x-direction, and a pixel step size of 0.25 mm in both directions. To provide a general view of the patterns created, Fig. 10 shows the four most extreme cases.

From Fig. 10 it can be noticed how the 2D-DFT peak distribution changes according to the fibers' orientations and distance in the C-scans. Obviously, contrary to the experimental sample, there is no background noise due to the almost perfectly defined texture in the simulated Cscans. These 2D-DFTs have been filtered and thresholded, as we previously did with the experimental sample. Then, the coordinates of the peaks defining the weave fabrics have been determined to calculate their corresponding distance and angle in the space domain. For all the 135 combinations, the obtained results are plotted in Fig. 11 (A) for the measured distances and Fig. 11 (B) for the measured angles. To analyze the accuracy of the measurement, we have compared each measured value with the real one to calculate the relative error (ε r ) as 6, the 2D-DFT patterns have been threshold at 80% of the maximum intensity after a filtering consisting of a 4 × 10 -2 mm -1 cutoff wavenumber on the horizontal and vertical axes. The main peaks are circled in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1

Numerical results of the measurement of the fiber bundles' orientation and distance at each ply. No significant variation of the orientations and spacing of the fibers is noticed over the layers. The "interface" column indicates between which two layers the time delay is placed. The round brackets indicate that the position corresponds to an inner point of a layer.

ε r = 100 ⃒ ⃒ X R -X C ⃒ ⃒ ⃒ ⃒ X R ⃒ ⃒ (5)
where X R and X C refer, respectively, to the real and calculated values of the magnitude X, which can be either the distance or angle. These errors have been estimated, as a function of the simulation distance and simulation angle and are presented in the bottom part of Fig. 11. The top plots in Fig. 11 show that the measured values (y-axis) are almost identical to those used to simulate the C-scan (x-axis), independently of the other magnitude value. This excellent agreement between both values confirms the performance of our method, leading to accurate quantitative data of the former structure. In addition, the accuracy has been quantitatively calculated by computing the relative error. In terms of the distance measurement, the obtained errors range from 1.4% to 2.8%, confirming the precision of the method. In addition, in the bottom plot of Fig. 11 (A), it can be observed how the relative error slightly increases as the distance between fibers increases (i.e., this approach is more accurate for closely packed weave fabrics). This change in the accuracy might be produced because the distance between the 2D-DFT peaks decreases when the distance in the space domain increases, and thus the peaks start to overlap and become more difficult to differentiate.

On the other hand, the relative error on the angle measurements ranges from 2.1% to 7.9%. In this case, the accuracy is better for highangle (i.e., almost horizontal) fibers orientation. It is the result of the error propagation in the angle calculation in the frequency domain [θ = tan -1 (v /u)]. The smaller the fibers' angle (i.e. the more vertical the orientation), the shorter the v/u ratio in the frequency domain. Therefore, a small error in the coordinates' approximation leads to a proportionally larger relative error when computing the ratio v/u. However, the larger the fibers' angular orientation, the larger the v/u ratio and, thus, the smaller the relative error despite having the same absolute error in the determination of the coordinates (u, v). In this regard, one can check that -contrary to Fig. 11 (A) -in Fig. 11 (B) the absolute error (or dispersion of results) is almost constant in all the angular range, and therefore the difference comes from the relative calculation discussed above.

These results confirm the intrinsic good performance of our algorithm when the experimental conditions do not add significant noise to the raw image. However, as it has been proved in Section 3.2, this procedure has the potential to reveal structures that are difficult to identify by eye. With this in mind, we have tested our approach with a batch of artificial C-scans with added noise and blurring effects, in order to replicate real C-scans with different grades of experimental noise. To do so, we have produced binary C-scans with θ SD = 45 • and d SD = 4 mm (the same as the real samples examined before) and added different types of noise and blurring: (i) First, we have added Gaussian noise by replacing random pixels with intensities within a gaussian distribution centered at 0.5 (the average intensity in the binary C-scans). The standard deviation (SD or σ) of the Gaussian distribution has been varied to see how the dispersion of intensities affects the quality of the measurements; (ii) Another strategy has been to add salt and pepper (S&P) noise, where different amounts of binary pixels (white and black, or 0 and 1) are randomly distributed through the image. The ratio between salt and pepper pixels has been kept constant (1:1) while the effect of the number of substituted pixels has been analyzed; (iii) Finally, we have also tested the effect of blurring effects in the image by using a Gaussian kernel with different SD and a constant truncation of 3. Fig. 12 shows the most extreme cases that we have produced by each approach. Each column corresponds to one of the methods (Gaussian noise, S&P noise, and Gaussian blurring, respectively); the top images show the weaker noise cases and in the bottom are represented the cases with the most significant alterations. In total, we have tested 21 noise values for each case to evaluate the accuracy of our approach. The results are depicted in Fig. 13. Furthermore, given that the two first approaches have a random component, we have performed 10 replicates with the aim of obtaining more statistically meaningful data. Concerning the Gaussian-type noise, we may first observe the high degree of perturbation reached in our simulations. However, we must notice that the accuracy is excellent while σ < 2. In this range, the error in the distance measurement is stable around 1.5%, while it oscillates between a 2 and 3% in the case of the angular measurement. Both are very low values and are in excellent agreement with the results shown before in Fig. 11 where the C-scans did not contain any noise. Therefore, we can conclude that the measurement is not affected by the noise until these levels. For σ ≥ 2 the error in the angular measurement remains stable, but it starts to oscillate for the distance. For the case σ = 2.42 (where the error reaches a maximum), the distances measured for the 10 replicates are respectively 3.92, 3.94, 25.84, 3.92, 3.92, 3.99, 3.99, 3.92 and 3.99 mm. 9 out of the 10 cases are close to the simulation value of 4 Fig. 12. Simulated C-scans with added noise and blurring effects. First column (A and B) corresponds to Gaussian noise with σ = 0.1 and σ = 3.0 respectively; Second column (C and D) corresponds to S&P noise with a 1% and 75% of the image pixels substituted; Third column (E and F) corresponds to C-scans blurred by using σ = 0.1 and σ = 3.0, respectively. mm, but there is one result at 25.84 mm which is responsible for the huge increase in the relative error. Given that there is almost no error in the angular measurement, the algorithm is measuring a pattern with a larger period in the same direction, but not having a significant error in the pattern measurement. Thus, for large values of noise there is a higher probability to hide the fiber bundles and make other less frequent periodicities more visible. However, we must highlight that the method is still working in most cases (90%) with high accuracy, while visual inspection is almost impossible.

In the S&P case [Fig. 13 (B)] we have tested how the amount of noise affects to the measurement accuracy. As it may be noticed, despite adding noise up to the 75% of the pixels in the image (which is a very extreme case) the error remains low and stable both for the angle and distance. Again, these error values are in good agreement with those reported in Fig. 11 and we can deduce that the algorithm keeps an excellent performance under all these conditions. Despite the fact that the S&P pixels might be predominant, there are still a number of pixels defining a pattern that can be detected through the 2D-DFT.

Finally, the last case illustrates how the performance of the method is affected by the blurring of the image. It seems clear that it is not affected by it because the relative error remains low, stable, and again in the error range obtained in Fig. 11 for ideal C-scans with the same characteristics. We deduce that blurring effects reduce the resolution of an image but, as long as it does not change the main patterns present, they can still be detected by the 2D-DFT.

With these results, it is seen that the proposed approach is robust against experimental noise and blurring, even in highly affected cases. Only a big dispersion in the intensity of the noise can hinder the detection of the smallest patterns in some cases.

Conclusions

In this work, we have described an image processing methodology that enables the visualization and quantitative measurement of the woven fabrics in GFRP. The method mainly consists of concatenating fast and straightforward averaging, binarization, filtering techniques, and Fourier analysis. It leads to an accurate representation of the composite's reinforcement microstructure from a raw C-scan THz image.

The algorithm has been demonstrated on real composite samples, in this case a 2/2 twill weave reinforcing a polyamide matrix. The analysis has succeeded in characterizing the microstructure-induced process, which resulted in an average bundles' orientation of -44.52 • and 41.61 • and distances of 4.05 and 3.86 mm at the surface ply. These results are in excellent agreement with the manufacturer's information and reveal imperfect alignment induced by the thermo-compression process.

Furthermore, the utility of this approach in extracting information on the inner structure of these laminates has been corroborated. Thanks to the combination of the penetration capability of THz waves and signal/ image processing, the orientation and spacing of the fibers at different plies has been determined, despite the lower SNR, revealing that the misalignment of one of the directions persists in depth. The error induced by the image processing algorithms has been quantified: It is less than 3% for the distance and within 2%-8% for the angle. The proposed algorithm has additionally demonstrated strong capabilities to deal with blurred or noisy raw images without any meaningful increase in the measurement error. Moreover, the robustness of the method and the values optimized for each variable has been demonstrated, as the accuracy is excellent even when a batch of artificially generated C-scans with a wide range of characteristics has been analyzed.

In conclusion, this work reports for the first time the measurement, in the different plies, of the characteristics of the woven fabric of laminated GFRP by a THz image-based method. The described method is of potential interest for the NDE of GFRP composites and opens the door to be used in a more versatile way for extracting quantitative information of the sample for quality control or in situations where the woven fabrics might be distorted. Specifically, the approach presented here gives accurate quantitative information on the weave's tow orientations and spatial period. This information may be helpful to determine whether materials have been fabricated to specification or damaged in service. Fig. 13. Accuracy of the proposed method for different types and levels of noise. (A) Shows the mean and standard deviation (shadowed region) values of the relative error of the measurement when Gaussian noise is added to the C-scan. The inner plot zooms into the region where 0<σ < 2 to facilitate its visualization; (B) Shows the mean and standard deviation (shadowed region) values of the relative error of the measurement when S&P noise is added to the C-scan; (C) Shows the relative error of the measurement when the C-scan is blurred by using a Gaussian kernel.

Fig. 2 .

 2 Fig. 2. Schematic description of the methodology followed by the developed image processing algorithm.

Fig. 3 .

 3 Fig. 3. B-scan obtained for a cross-section along the short axis (y = 30 mm) at the center of the sample. (A) Negative amplitudes have been set to zero in order to increase the contrast. (B) Binary B-scan obtained with a global threshold of 2.5% of the maximum intensity. (C) THz pulse amplitude in time (A-scan) is located at y = 30 mm and x = 10 mm. The THz signal impinges from below (small time delays), and the horizontal dashed lines locate the interfaces found in the sample.

Fig. 4 .

 4 Fig. 4. (A) Raw and (B) averaged with n = 2 neighbors C-scans at the sample surface (4 ps). The two perpendicular fiber bundles (at ±45 • ) are clearly visible.

Fig. 5 .Fig. 6 .

 56 Fig. 5. (A) Binary C-scan obtained with an adaptive thresholding with 3 neighbors, and (B) smoothed binary C-scan obtained after averaging with n = 2 and thresholding at 50% of the intensity to increase the fibers' definition. The matrix part is in blue, and the reinforcement tows are in red in each picture. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .

 7 Fig. 7. Comparison between (A) inverse 2DFFT representing the mean fibers' structure after processing the image, and (B) experimental averaged C-scan, both at 4 ps. (C) graphically shows the space-domain distance (d SD ) and angle (θ SD ) criteria used in the calculation. Similar to the experimental C-scan, one orientation is more intense than the other in the inverse 2DFFT.

Fig. 8 .

 8 Fig. 8. (A) Binary B-scan with different selected time delays marked with dark lines, and its corresponding averaged C-scans at (B) 6.0 ps, (C) 12.0 ps, (D) 19.5 ps, and (E) 27.5 ps. The ±45 • weave pattern is still visible in the C-scans in depth.

Fig. 9 .

 9 Fig. 9. Filtered modulus square 2D-DFT patterns and inverse-2D-DFT sample reconstructions obtained from the averaged C-scans at (A) 6.0 ps, (B) 12.0 ps, (C) 19.5 ps, and (D) 27.5 ps. Similar to what is presented in Fig. 6, the 2D-DFT patterns have been threshold at 80% of the maximum intensity after a filtering consisting of a 4 × 10 -2 mm -1 cutoff wavenumber on the horizontal and vertical axes. The main peaks are circled in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Simulated binary C-scans and their corresponding modulus squared 2DFFT images for the conditions: (A) 10 • and 2 mm, (B) 10 • and 6 mm, (C) 80 • and 2 mm, and (D) 80 • and 6 mm. The matrix is in blue and the fibers in red in the simulated C-scans. The four main peaks are circled in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Acknowledgements

JCR, PA, DR, AL, and DSC gratefully acknowledge the financial support of the Conseil Régional du Grand Est, of the Fonds Européen de Développement Régional (FEDER), of the Institut Carnot ARTS, and of CPER SusChemProc. DR gratefully acknowledges the support of the Chair in Photonics. FM and PP would like to thank Stellantis for their support in the framework of the OpenLab Materials and Processes involving Stellantis group, Arts et Métiers (ENSAM) and Georgia Tech-Europe. They are grateful to DuPont de Nemours for providing the materials of this study.

Data availability

The data that has been used is confidential.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Konstantin Zimenko

Andrey Polyakov

Denis Efimov

Homogeneous Systems Stabilization Based on Convex Embedding ⋆

The paper develops control algorithms for a class of affine nonlinear systems using the so-called canonical homogeneous representation. It is demonstrated that such a representation exists for any homogeneous vector field bounded on the unit sphere. It is shown that canonical homogeneous representation is useful for LMI-based control design and stability analysis of nonlinear systems. Theoretical results are supported by numerical examples.

Introduction

Homogeneous systems are widely used in stability analysis, nonlinear control and observer design (see, for example, [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF]- [START_REF] Zimenko | Robust Feedback Stabilization of Linear MIMO Systems Using Generalized Homogenization[END_REF]). Additional interest to homogeneity (homogenization) based approaches consists in their practically important robustness (see [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF]- [START_REF] Andrieu | Homogeneous approximation, recursive observer design, and output feedback[END_REF]), e.g., Input-to-State Stability, non-Lipschitz disturbance rejection, etc. Stability (stabilizability) analysis for homogeneous systems is the subject of many studies. In [START_REF] Zubov | Methods of A.M. Lyapunov and Their Applications[END_REF], [START_REF] Rosier | Homogeneous Lyapunov function for homogeneous continuous vector field[END_REF] it is shown that for a stable (weighted) homogeneous system there exists a corresponding homogeneous Lyapunov function. The papers [START_REF] Kawski | Homogeneous stabilizing feedback laws[END_REF], [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF] present conditions for the existence of stabilizing homogeneous feedbacks. Note that for stable homogeneous systems rates of convergence (exponential, finite-time, nearly fixed-time) can be determined based on the sign of homogeneity degree (see, e.g., [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF], [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF], [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF]). The paper [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF] develops the result [START_REF] Rosier | Homogeneous Lyapunov function for homogeneous continuous vector field[END_REF] for the class of generalized homogeneous systems. It is shown that any generalized homogeneous system is diffeomorphic to a standard homogeneous one, i.e., there exists a change of coordinates, which transforms an asymptotically stable (generalized) homogeneous system to a quadrati-⋆ The material in this paper was partially presented at the IFAC World Congress 2020, Berlin, Germany [START_REF] Zimenko | Robust Stabilization of Control Affine Systems with Homogeneous Functions[END_REF] and European Control Conference 2021, Rotterdam, the Netherlands [START_REF] Zimenko | On Convex Embedding and Control Design for Nonlinear Homogeneous Systems[END_REF]. This work is supported by RSF under grant 21-71-10032. Corresponding author K. Zimenko.

Email addresses: konstantin.zimenko@itmo.ru (Konstantin Zimenko), andrey.polyakov@inria.fr (Andrey Polyakov), denis.efimov@inria.fr (Denis Efimov).

cally stable one. In some particular cases (mostly for linear and close to linear systems, e.g. [START_REF] Zimenko | Robust Feedback Stabilization of Linear MIMO Systems Using Generalized Homogenization[END_REF]) these results jointly with implicit Lyapunov function approach allow to formulate useful stability (stabilizability) conditions in terms of linear matrix inequalities (LMIs).

In order to obtain simple LMI-based stability (stabilizability) conditions and control tuning algorithms the embedding of a nonlinear system into a linear differential inclusion can be used (e.g., see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], [START_REF] Pyatnitskiy | Criteria of asymptotic stability of differential inclusions and periodic motions of time-varying nonlinear control systems[END_REF], [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H∞ performance analysis[END_REF]). The similar idea was utilized in this paper to introduce a canonical representation for homogeneous systems. Two algorithms for derivation of the canonical homogeneous representation are presented for smooth (differentiable) and non-smooth vector fields, respectively. In the latter case, the boundedness on the unit sphere of the homogeneous vector field is required. Based on this representation a stability analysis method and two control design algorithms with LMI-based stability (stabilizability) conditions are proposed for a class of non-linear homogeneous systems. These results are extended for some non-homogeneous systems as well.

The preliminary results of this paper are discussed at IFAC World Congress [START_REF] Zimenko | Robust Stabilization of Control Affine Systems with Homogeneous Functions[END_REF] and European Control Conference [START_REF] Zimenko | On Convex Embedding and Control Design for Nonlinear Homogeneous Systems[END_REF]. The key differences are as follows:

• the proofs of all claims are given;

• new LMI-based stability conditions for non-linear homogeneous systems are obtained (Lemma 4); • new schemes simplifying control design are presented (Remarks 4, 8, 9, Algorithms 1, 2, 3); • new numerical examples demonstrating efficiency of the proposed control design schemes are considered.

Notation: R + = {x ∈ R : x > 0}, where R is the field of real numbers; the symbol 1, m is used to denote a sequence of integers 1, ..., m; diag{a i } n i=1 is the diagonal matrix with the elements a i on the main diagonal; C n (X, Y ) is the set of continuously differentiable (at least up to the order n) functions X → Y , where X and Y are open subsets of R n ; e s (i) =

  0 • • • 0 ith 1 0 • • • 0 s components   T ∈ R s , s ≥ 1; l n
is the vector of dimension n with all elements equal 1; the inequality P > 0 (P ≥ 0) means that a symmetric matrix

P = P T ∈ R n×n is positive definite (positive semi-definite); ∥ • ∥ denotes a norm; S = {x ∈ R n : ∥x∥ = 1} is the unit sphere in R n .
2 Problem statement

Consider the system ẋ = f 0 (x) + m i=1 u i f i (x), (1)
where x ∈ R n are states, u = (u 1 , ..., u m ) T ∈ R m are controls, f 0 (0) = 0, f 0 and f i are homogeneous (in the sense explained below) or almost homogeneous (presented by product of a homogeneous vector field and a bounded scalar function). The state vector x is assumed to be measured.

The main goal is to derive a LMI-based stability condition for the open-loop system (u i = 0) and to propose LMI-based design of a feedback control algorithm for the system [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF]. In order to obtain constructive stability conditions and control tuning algorithms, a specific canonical representation of a homogeneous system needs to be obtained using the so-called convex embedding (see, e.g., [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]Chapter 4]).

3 Preliminaries

Stability Notions

Consider the system

ẋ(t) = f (x(t)), x(0) = x 0 , t ≥ 0, (2)
where

x(t) ∈ R n is the state vector, f : R n → R n is continuous, f (0) = 0.
Definition 1 [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF], [START_REF] Bhat | Finite-time stability of continuous autonomous systems[END_REF], [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF] A set M ⊂ R n is said to be globally finite-time attractive for (2) if any solution x(t, x 0 ) of (2) reaches M in some finite time moment t = T (x 0 ) and remains there for all t ≥ T (x 0 ), where T : R n → R + ∪{0} is the so-called settling-time function.

It is fixed-time attractive if in addition the settlingtime function T (x 0 ) is globally bounded by some number T max > 0. If M = {0} is globally finite/fixed-time attractive and Lyapunov stable then the origin is said to be globally finite/fixed-time stable. The origin is said to be nearly fixed-time stable if it is globally Lyapunov stable and any neighborhood of the origin is fixed-time attractive.

Theorem 1 [START_REF] Bhat | Finite-time stability of continuous autonomous systems[END_REF], [START_REF] Nakamura | Smooth Lyapunov functions for homogeneous differential inclusions[END_REF] Suppose there exist a positive definite C 1 function V defined on an open neighborhood of the origin D ⊂ R n and real numbers C > 0 and σ ≥ 0, such that for x ∈ D\{0} the condition V (x) ≤ -CV σ (x) is true for the system [START_REF] Ryan | Universal stabilization of a class of nonlinear systems with homogeneous vector fields[END_REF]. Then depending on the value σ the origin is stable with different types of convergence:

• if σ = 1, the origin is asymptotically stable;

• if 0 ≤ σ < 1, the origin is finite-time stable and T (x 0 ) ≤ 1 C(1-σ) V 1-σ (x 0 ); • if σ > 1 the origin is nearly fixed-time stable and for every ε ∈ R + , the set B = {x ∈ D : V (x) ≤ ε} is fixed-time attractive with T max = 1 C(σ-1)ε σ-1 .
If D = R n and V is radially unbounded, then the mentioned stability properties of the system (2) hold globally.

'Universal' Stabilizer for Control Affine Systems

Consider a system in the form [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF], where f 0 , f i ∈ C ∞ , and f 0 (0) = 0. Let this system admit a control-Lyapunov function V that is, a smooth, proper, and positive definite function

R n → R + ∪ {0} so that inf u∈R m ∂V (x) ∂x f0(x)+ ∂V (x) ∂x m i=1 uifi(x) < 0, ∀x ̸ = 0. ( 3
)
Definition 2 [START_REF] Sontag | A 'universal' construction of Artstein's theorem on nonlinear stabilization[END_REF] The control-Lyapunov function V satisfies the small control property if for each ε > 0 there is a δ > 0 such that, if x ̸ = 0 satisfies ∥x∥ < δ, then there is some u ∈ R m with ∥u∥ < ε such that

∂V (x) ∂x f 0 (x) + ∂V (x) ∂x m i=1 u i f i (x) < 0.
Definition 3 [START_REF] Sontag | A 'universal' construction of Artstein's theorem on nonlinear stabilization[END_REF] Let k : R n → R m be a mapping, smooth on R n \ {0} and k(0) = 0. This is a smooth feedback stabilizer provided that the closed-loop system

ẋ = f 0 (x) + m i=1 k i (x)f i (x) is globally asymptotically stable, where k = (k 1 , ..., k m ) T .
Theorem 2 [START_REF] Sontag | A 'universal' construction of Artstein's theorem on nonlinear stabilization[END_REF] Let V be a smooth control-Lyapunov function V . Then (a) there exists a smooth feedback stabilizer k; (b) if V satisfies the small control property, then k can be chosen to be also continuous at the origin, e.g., in the form

k(x) =    a(x)+ √ a 2 (x)+|B(x)| 4 |B(x)| 2 B(x), if |B(x)| ̸ = 0, 0, if |B(x)| = 0, where a(x) = ∂V (x) ∂x f 0 (x), B(x) = (b 1 , ..., b m ) T , b i (x) = ∂V (x)
∂x f i (x) for i = 1, m.

Generalized Homogeneity

Homogeneity is a certain invariance of a mathematical object with respect to a group of transformations called dilations. In this paper we deal with the one-parameter group d(s) : R n → R n of linear dilations given by d(s) =

e G d s = +∞ i=0 s i G i d i! , s ∈ R, where G d ∈ R
n×n is an anti-Hurwitz matrix (i.e., -G d is Hurwitz) called the generator of the dilation.

Definition 4 [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF] The dilation d is said to be strictly monotone if ∃β > 0 : ∥d(s)∥ < e βs for s ≤ 0.

Theorem 3 [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF] If d is a dilation in R n , then it is strictly monotone with respect to the norm ∥x∥ P = √

x T P x for x ∈ R n and P satisfying

P G d + G T d P > 0, P = P T > 0. ( 4
)
Theorem 3 shows that any dilation d is strictly monotone if R n is equipped with the weighted Euclidean norm ∥x∥ P = √ x T P x provided that the matrix P satisfies (4).

Definition 5 [START_REF] Kawski | Homogeneous Feedback Stabilization[END_REF] A vector field f :

R n → R n (a function h: R n → R) is said to be d-homogeneous of degree ν ∈ R if f (d(s)x) = e νs d(s)f (x), ∀x ∈ R n , ∀s ∈ R. (5) (resp. h(d(s)x) = e νs h(x), ∀x ∈ R n , ∀s ∈ R.) (6)
Let F d (R n ) (respectively H d (R n )) be the set of vector fields R n → R n (respectively functions R n → R) satisfying the identity (5) (respectively (6)), which are continuous on R n \ {0}. Let deg

F d (f ) (respectively deg H d (f )) denote the homogeneity degree of f ∈ F d (R n ) (respec- tively h ∈ H d (R n )).
A special case of homogeneous function is a homogeneous norm [START_REF] Kawski | Geometric Homogeneity and Stabilization[END_REF], [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF], [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF]: a continuous positive definite d-homogeneous function of degree 1. For monotone dilations we define the canonical homogeneous norm

∥ • ∥ d : R n → R + ∪ {0} as ∥x∥ d = e sx , where s x ∈ R : ∥d(-s x )x∥ = 1, x ̸ = 0 (7)
and, by continuity, we assign ∥0∥ d = 0 (see, [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF] for more details). Notice that (7) also implies

∥d(-ln ∥x∥ d )x∥ = 1. ( 8
)
Lemma 1 [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF] If d is a strictly monotone dilation, then

• the homogeneous norm ∥ • ∥ d is Lipschitz continuous outside the origin;

• if the norm ∥ • ∥ is smooth outside the origin, then the homogeneous norm ∥ • ∥ d is also smooth outside of the origin and

∂∥x∥ d ∂x = ∥x∥ d ∂∥z∥ ∂z z=d(-s)x ∂∥z∥ ∂z z=d(-s)x G d d(-s)x s=ln ∥x∥ d . (9)
Lemma 2 [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF] The vector field

f ∈ F d (R n ) is Lipschitz continuous (smooth) on R n \ {0}
if and only if it satisfies a Lipschitz condition (it is smooth) on the unit sphere S, provided that d is strictly monotone on R n .

If a function (or a vector field) is smooth, then homogeneity is inherited by its derivatives in a certain way.

The following lemma is an extension of the well-known Euler's homogeneous function theorem to the case of linear homogeneity.

Lemma 3 [1, pages 211, 212], [17, Corollary 2], [18] If f ∈ F d (R n ) ∩ C 1 (R n \ {0}, R n ) and h ∈ H d (R n ) ∩ C 1 (R n \ {0}, R), then e deg H d (h)s ∂h(x) ∂x = ∂h(z) ∂z z=d(s)x d(s), (10)
e deg F d (f )s d(s) ∂f (x) ∂x = ∂f (z) ∂z z=d(s)x d(s), (11) ∂h
(x) ∂x G d x = deg H d (h)h(x), (12)
∂f (x) ∂x G d x = (deg F d (f )I n + G d )f (x) (13)
for x ∈ R n \ {0} and s ∈ R.

Due to topological equivalence of d-homogeneous system to a standard homogeneous one [START_REF] Polyakov | Sliding Mode Control Design Using Canonical Homogeneous Norm[END_REF], all results about stability and robustness of standard and weighted homogeneous systems hold for d-homogeneous systems as well. For example, the convergence rate of homogeneous systems can be assessed by its homogeneity degree:

Theorem 4 [START_REF] Nakamura | Smooth Lyapunov functions for homogeneous differential inclusions[END_REF] An asymptotically stable d-homogeneous

system ẋ = f (x), f : R n → R n is finite-time stable if deg F d (f ) < 0 (nearly fixed-time stable if deg F d (f ) > 0).
The homogeneity provides many other advantages to analysis and design of nonlinear control system. For example, homogeneous control systems are Input-to-State Stable (ISS) with respect to measurement noise and additive disturbances (see, for example, [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF], [START_REF] Ryan | Universal stabilization of a class of nonlinear systems with homogeneous vector fields[END_REF], [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF], [START_REF] Andrieu | Homogeneous approximation, recursive observer design, and output feedback[END_REF]).

Convex embedding

Definition 6 [31, pages 10, 12] A subset S of R n is said to be convex if (1 -λ)x + λy ∈ S whenever x ∈ S, y ∈ S and λ ∈ (0, 1). The convex hull of a given set S may be defined as the intersection of all convex sets containing S and is denoted as co S.

Theorem 5 (Carathéodory's theorem) [32, page 62] For any bounded closed set S ⊂ R n any point x ∈ co S can be represented in the form x = k i=0 α i x i , where

x i ∈ S, α i ∈ R + ∪ {0} : k i=0 α i = 1 and k ≤ n. Let A = a11 ••• a1n . . . . . . . . . an1 ••• ann , where a ij ≤ a ij ≤ a ij , i, j = 1, n.
Let us show how the convex embedding concept can be used in order to represent the matrix A as a convex combination. Let us define the set A = a = (a 11 , ..., a 1n , ..., a nn ) :

a ij ≤ a ij ≤ a ij , i, j = 1, n . The set A ⊂ R n 2 is a bounded convex domain of which the set of vertices is defined by V a = ϑ = (ϑ 11 , ..., ϑ 1n , ..., ϑ nn ) : ϑ jl ∈ {a ij , a ij }, i, j = 1, n .
Then applying the Carathéodory's theorem, for any a ∈ A we obtain a

= k i=0 α i a i , where a i ∈ V a , k i=0 α i = 1, k ≤ n 2 or A = k i=0 α i A i , where A i ∈ V A , V A = ϑ = ϑ11 ••• ϑ1n . . . . . . . . . ϑn1 ••• ϑnn ∈ R n×n : ϑ ij ∈ {a ij , a ij }, i, j = 1, n .
In [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], [START_REF] Zemouche | Observers for a class of Lipschitz systems with extension to H∞ performance analysis[END_REF] the convex embedding concept was used for embedding of a differentiable nonlinear system into a linear differential inclusion.

4 Convex embedding for homogeneous systems

Convex Embedding for Differentiable Homogeneous Systems

Let us consider a system in the form

ẋ(t) = f (x(t)), t ≥ 0 ( 14
)
where

x(t) ∈ R n , f ∈ F d (R n ), deg F d (f ) = ν ∈ R, f is locally bounded. Let f ∈ C 1 (R n \{0}, R n ).
For i, j = 1, n define the bounds

g ij = sup y∈R n :∥y∥=1 ∂f i (z) ∂z j z=y , g ij = inf y∈R n :∥y∥=1 ∂f i (z) ∂z j z=y , that are always exist due to f ∈ C 1 (R n \{0}, R n ).
Denote the set of vertices defined by

Vg = ϑ = ϑ 11 ••• ϑ 1n . . . . . . . . . ϑ n1 ••• ϑnn ∈ R n×n : ϑij ∈ {g ij , g ij }, i, j = 1, n . (15)
Proposition 1 Let the matrix νI n + G d be invertible. There exist

α i : S → R + ∪ {0} : N i=1 α i (x) = 1 and M i ∈ V g , i = 1, ..., N ≤ 2 n 2
such that every trajectory x(t) of ( 14) is also a trajectory of

ẋ = ∥x∥ ν d N i=1 αi(p(x))d(ln ∥x∥ d )Aip(x), t ≥ 0 (16)
as long as x(t) ̸ = 0, where

A i = (νI n + G d ) -1 M i G d and p(x) = d(-ln ∥x∥ d )
x is the so-called homogeneous projection of the vector x to the unit sphere.

Proof. Due to νI n +G d is invertible then by [START_REF] Kawski | Homogeneous Feedback Stabilization[END_REF] we have

f (x) = (νI n + G d ) -1 ∂f (x) ∂x G d x (17)
for all x ∈ R n \ {0}. Consider [START_REF] Kawski | Geometric Homogeneity and Stabilization[END_REF] with s = -ln ∥x∥ d then one can obtain

∂f (x) ∂x = e -νs d(-s) ∂f (z) ∂z z=d(s)x d(s) = e ν ln ∥x∥ d d(ln ∥x∥ d ) ∂f (z) ∂z z=p(x) d(-ln ∥x∥ d ) = ∥x∥ ν d d(ln ∥x∥ d ) ∂f (z) ∂z z=p(x) d(-ln ∥x∥ d )
Substituting the obtained expression for ∂f (x) ∂x to ( 17) and taking into account

(νI n + G d ) -1 d(ln ∥x∥ d ) = d(ln ∥x∥ d ) (νI n + G d ) -1 and d(-ln ∥x∥ d )G d = G d d(-ln ∥x∥ d ) we have f (x) = (νI n + G d ) -1 ∂f (x) ∂x G d x = ∥x∥ ν d (νI n + G d ) -1 d(ln ∥x∥ d ) ∂f (z) ∂z z=p(x) × d(-ln ∥x∥ d )G d x = ∥x∥ ν d d(ln ∥x∥ d ) (νI n +G d ) -1 ∂f (z) ∂z z=p(x) G d p(x) for all x ∈ R n \ {0}. Taking into account that ∂f (z) ∂z z=p(x) = n j=1 n l=1 Q jl ∂f j (z) ∂z l z=p(x)
with

Q jl = e n (j)e T n (l), for θ jl ∈ [g jl , g jl ] we obtain f (x)=∥x∥ ν d d(ln ∥x∥ d ) (νI n + G d ) -1 × n j=1 n l=1 Q jl ∂fj (z) ∂z l z=p(x) G d p(x) = ∥x∥ ν d d(ln ∥x∥ d ) (νI n +G d ) -1 × n j=1 n l=1 Q jl θ jl G d p(x).
Finally, using Carathéodory theorem we can represent the term

n j=1 n l=1 Q jl θ jl (z) ∈ co V g as N i=1 α i (z)M i with M i ∈ V g .
Then all trajectories of ( 14) are also trajectories of the system [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF]. ■

Example 1 Consider the system (14

) with f (x) = -x1x3 -x1 x2
that is d-homogeneous of degree 1, where

G d = 3 0 0 0 2 0 0 0 1
is a generator of the dilation. According to Proposition 1 the system can be represented in the form [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF]

with ∥x∥ = √ x T x, M i ∈ {-1,1} 0 {-1,1} -1 0 0 0 1 0 and A i = {-0.75,0.75} 0 {-0.25,0.25} -1 0 0 0 1 0 for i = 1, 4.
In [START_REF] Zimenko | Robust Stabilization of Control Affine Systems with Homogeneous Functions[END_REF] the representation of the function f in the form like ( 16) was obtained by another constructions. Comparing with [START_REF] Zimenko | Robust Stabilization of Control Affine Systems with Homogeneous Functions[END_REF] the given result allows to obtain vertices A i significantly closer to each other, that may simplify the stability analysis and control design for essentially nonlinear homogeneous systems.

Example 2 Consider the system (14) with f (x) = 16), where

A i = 0 {-0.6,0.6} {-0.4,0.4} 0 0 1 0 0 0 , i = 1, 4.
Another advantage with respect to [START_REF] Zimenko | Robust Stabilization of Control Affine Systems with Homogeneous Functions[END_REF] is that the differentiability at the origin is not required. However, the condition f ∈ C 1 (R n \ {0}, R n ) also restricts sufficiently the class of homogeneous functions for representation in the form [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF]. This condition, however, can be relaxed as shown below.

Convex Embedding for Non-Smooth Homogeneous Systems

Let f ∈ F d (R n ) in ( 14) be locally bounded on R n but possibly non-differentiable on a unit sphere and let P ∈ R n×n satisfy (4). Define the bounds

q ij = sup y∈R n :∥y∥ P =1 e T n (i)f (y)yj y T P G d y , q ij = inf y∈R n :∥y∥ P =1 e T n (i)f (y)yj y T P G d y
and denote the corresponding set of vertices V q as above (see, the formula (15)) using the sets {q ij , q ij } instead of {g ij , g ij }.

Proposition 2 There exist α i : S → R + ∪ {0} : N i=1 α i (x) = 1 and M i ∈ V q , i = 1, ..., N ≤ 2 n 2 such that every trajectory x(t) of ( 14) is also a trajectory of (16) as long as x(t) ̸ = 0, where

A i = M i P G d .
Proof. Since P satisfies (4), then the dilation d is strictly monotone on R n equipped with the norm ∥x∥ P = √

x T P x.

Note that by [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF], for s = -ln ∥x∥ d we have

f (x) = ∥x∥ ν d d(ln ∥x∥ d )f (p(x)) = ∥x∥ ν d d(ln ∥x∥ d )diag{e T n (i)f (p(x))} n i=1 l n .
Then due to ∂∥z∥ d ∂z z=p(x)

G d p(x) = 1 (see [START_REF] Moulay | Finite-time stability and stabilization: State of the art[END_REF] and ( 12))

and ∂∥z∥ d ∂z z=p(x)

= z T P z T P G d z (see ( 9)) we obtain

f (x) = ∥x∥ ν d d(ln ∥x∥ d )diag{e T n (i)f (p(x))} n i=1 l n = ∥x∥ ν d d(ln ∥x∥ d )         e T n (1)f (z)z T z T P G d z z=p(x)
. . .

e T n (n)f (z)z T z T P G d z z=p(x)         P G d p(x).
Finally, using convex embedding procedure on the unit sphere we finish the proof. ■

Remark 1 If one can represent the system as ẋ =

q i=1 f i (x) with f i ∈ F d (R n ), q > 1,
then in order to obtain vertices closer to each other the given results (e.g., based on Proposition 1 and/or Proposition 2) can be used for each term independently. In this case the representation ( 16) stays true with N ≤ q2 n 2 , N i=1 α i (x) = q. In the case q = 2 and f 1 = Cx the representation ( 16 that is d-homogeneous of degree -1/2, where G d = ( 0.5 0.5 0 1 ). Let us rewrite the system as ẋ = f 1 (x) + f 2 (x), where

f 1 (x) = x1-x2 x1-x2 , f 2 (x) = |x2| 1/2 |x2| 1/2 .
Then using the results of Proposition 2 for representing f 2 (x) with ∥x∥ =

√

x T x the following matrices can be obtained [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF], and according to Remark 1 the system under consideration can be represented as

M i = {-1.502,1.502} {-1.235,1.235} {-1.502,1.502} {-1.235,1.235} , i = 1,
ẋ = ∥x∥ ν d 16 i=1 α i (p(x))d(ln ∥x∥ d ) 1 -1 1 -1 + A i p(x), where 16
i=1 α i (x) = 1, A i = M i P G d and ν = -0.5.

Generalization to Some Non-Homogeneous Systems

The results of Propositions 1, 2 can be used for some nonhomogeneous dynamics representation in the form [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF].

Consider the system ẋ = b(x)f (x), (18)
where b :

R n → R is bounded on R n and f ∈ F d (R n ) is bounded on S, f (0) = 0.
Define the bounds

p ij = sup x∈R n b(x) e T n (i)f (z)zj z T P G d z z=p(x) , p ij = inf x∈R n b(x) e T n (i)f (z)zj z T P G d z z=p(x)
and denote the corresponding set of vertices V p (if f (x) is differentiable on the unit sphere define the bounds as

p ij = sup x∈R n b(x) ∂fi(z) ∂zj z=p(x)
,

p ij = inf x∈R n b(x) ∂fi(z) ∂zj z=p(x)
).

Corollary 1 Every trajectory x(t) of ( 18) is also a trajectory of ( 16) as long as x(t) ̸ = 0, where A i is defined as in Proposition 2 (Proposition 1), provided that f (x) is bounded on the unit sphere (resp., f (x) is differentiable on the unit sphere).

Proof. The proof is a straightforward consequence of Proposition 2 (resp., Proposition 1). ■ Remark 2 If one can represent the system under consideration

ẋ = M i=1 b i (x)f i (x), (19)
with bounded b i : R n → R, f i ∈ F d (R n ), M ≤ n, then according to Remark 1 in order to obtain vertices closer to each other the given results can be used for each term independently (see Example 4). Similarly, the systems ẋ = diag{b i (x)} n i=1 f (x) with bounded b i : R n → R,

f i ∈ F d (R n ) can be considered under the commutativity conditions that diag{b i (x)} n i=1 d(ln ∥x∥ d ) = d(ln ∥x∥ d )diag{b i (x)} n i=1 . ( 20
)
Example 4 Consider the system ẋ = , i = 1, 8. Note that the choice of the norm (matrix P in ∥ • ∥ P ) significantly affects the location of vertices, that should be taken into account for control design based on the use of the canonical representation (see Section 5).

cos(x1)(x2x 2 3 +x 2 2 ) x1 x2+x 2 3 . Note that f (x) =

Stability Analysis and Control Design

LMI-based Stability Conditions

The system representation ( 16) can be useful for stability analysis, where stability conditions are presented in the form of LMIs. Consider the system ( 14) represented in the canonical form [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF].

Lemma 4 Let ν < 0 (resp., ν > 0). If H ∈ R n×n , η, δ ∈ R + , θ min , θ max ∈ R satisfy the system of LMIs θmin(P G d +G T d P ) ≤ H -HG d -G T d H ≤ θmax(P G d +G T d P ), (θmaxP + H) Ai + A T i (θmaxP + H) ≤ -δP, (θminP + H) Ai + A T i (θminP + H) ≤ -δP, ηP ≥ H > 0, θmax ≥ θmin > 0 ( 21
)
for i = 1, N , then the system ( 16) is finite-time (resp., nearly fixed-time) stable.

Proof. Let ν ≥ -1. Choose a candidate Lyapunov function in the form

V (x) = ∥x∥ d x T d T (-ln ∥x∥ d )Hd(-ln ∥x∥ d )x. ( 22
)
Then we have:

V = p T (x)Hp(x) 2∥x∥ d p T (x)P d(-ln ∥x∥ d ) p T (x)(P G d +G T d P )p(x) ẋ +2∥x∥ d p T (x)H -∥x∥ -1 d G d p(x) ∂∥x∥ d ∂x +d(-ln ∥x∥ d ) ẋ = 2∥x∥ d p T (x)Hp(x) p T (x)(P G d +G T d P )p(x) p T (x)P d(-ln ∥x∥ d ) ẋ +2∥x∥ d p T (x)H -2G d p(x) p T (x)P p T (x)(P G d +G T d P )p(x) +In × d(-ln ∥x∥ d ) ẋ = 2∥x∥ d s T p T (x)(H-HG d -G T d H)p(x) p T (x)(P G d +G T d P )p(x) P +H d(-ln ∥x∥ d ) ẋ = 2∥x∥ 1+ν d p T (x) p T (x)(H-HG d -G T d H)p(x) p T (x)(P G d +G T d P )p(x) P + H × N i=1 αi(p(x))Aip(x) = 2∥x∥ 1+ν d p T (x) (θ(p(x))P + H) N i=1 αi(p(x))Aip(x),
where

p T (x)(H-HG d -G T d H)p(x) p T (x)(P G d +G T d P )p(x) = θ(p(x)) ∈ [θ min , θ max ].
Using convex embedding one can show that (θ(p(x))P + H)

N i=1 α i (p(x))A i + N i=1 α i (p(x))A T i (θ(p(x))P + H) ≤ -δP, (23)
if the inequalities

(θ max P + H) A i + A T i (θ max P + H) ≤ -δP, (θ min P + H) A i + A T i (θ min P + H) ≤ -δP, (24)
are satisfied for i = 1, N . Finally, if ( 23) is satisfied then we have

V ≤ -δ∥x∥ 1+ν d ≤ - δ η 1+ν V 1+ν . ( 25
)
In the case ν < -1 choose the candidate Lyapunov function in the form Ṽ = V µ , where µ ≥ -ν. Then according to [START_REF] Zimenko | Robust Feedback Stabilization of Linear MIMO Systems Using Generalized Homogenization[END_REF] we can obtain V ≤ -µδ η 1+ν

Ṽ µ+ν µ .

Finally, according to Theorem 1 the system is finite-time (nearly fixed-time) stable if ν < 0 (ν > 0). ■ Remark 3 According to Theorem 1 the following estimates of the settling time function can be obtained:

• for ν < 0 the origin is finite-time stable and the settling time function is bounded as T (x 0 ) ≤ -

η 1+ν V -ν 0 δν
, where V 0 = V (x 0 ); • for ν > 0 the system is nearly fixed-time stable and for any ε ∈ R + , the set

B = {x ∈ R n : V (x) ≤ ε} is fixed-time attractive with T max = η 1+ν δν ε -ν .
Remark 4 Note that the matrix P cannot be used as a solution of the inequalities ( 21) while it is used for representation of f (x) in the form ( 16), i.e., the matrices A i depend on P satisfying (4). In this regard, feasibility of ( 21) may depend on the choice of the matrix P , then the following heuristic algorithm can be aplied:

Algorithm 1 INITIALIZATION: take P = I n ; compute matrices A i ; STEP: While ( 21) is not feasible Solve the LMI

X > 0, A i X + XA T i + δX < 0, i = 1, N (26)
for some δ ∈ R + with respect to X ∈ R n×n . take P = X -1 ; recompute matrices A i ; end Example 5 Consider the double integrator system ẋ1 ẋ2 = ( x2 u ) with the finite-time stabilizing control [START_REF] Zimenko | On simple scheme of finite/fixed-time control design[END_REF] u

(x) = c l-µ (x)K c -l-µ (x) 0 0 c -l (x) x, c(x) = |x 1 | 2l+2µ+0.1 l+µ + |x 2 | 2l+2µ+0.1 l 1 2l+2µ+0.1 , K = ( -12 -6 ) , µ ∈ (0, 1], l ∈ R + .
The system is homogeneous of degree -µ. In the stability analysis of such systems, LMI-based conditions are often satisfied for a sufficiently small µ only (i.e. for the case the control is close to the linear one) or for a sufficiently big l (which can lead to a conservative estimate of T (x 0 )). Applying Proposition 1 for the system representation in the form ( 16) and using Lemma 4 with Algorithm 1 it can be shown that the closed-loop system is finite-time stable, in particular, for µ = 1 and l = 4.5 with the estimation of the settling time T (x 0 ) ≤ 4.925 for x 0 = (0 1) T (according to simulation T (x 0 ) ≈ 1.65).

Applying the result of [START_REF] Zimenko | On simple scheme of finite/fixed-time control design[END_REF] for the same parameters the obtained estimation of the settling time is too conservative: T (x 0 ) ≤ 1128.5.

LMI-based Design of Sontag's Universal Controller

In this section, we demonstrate how the canonical form can be utilized for control design of the system (1), where f 0 , f i are homogeneous functions. Using the representation ( 16) the control for (1) can be obtained based on the scheme of a 'universal' stabilizing control design given in the work [START_REF] Sontag | A 'universal' construction of Artstein's theorem on nonlinear stabilization[END_REF]. Let us define the matrices A i,j , i = 0, m, j = 1, N , and α i,j (p(x)), such that f i in (1) can be represented as

f i (x) = Ni j=1 α i,j (p(x))∥x∥ deg F d (fi) d d(ln∥x∥ d )A i,j p(x), (27)
where

N i ≤ 2 n 2 . Theorem 6 Let deg F d (f 0 ) ≥ -1, and deg F d (f 0 ) ≥ deg F d (f i ) for i = 1, m. If the system of LMIs P A0,j +A T 0,j P < m i=1 τi(P A i,l i +A T i,l i P ), j = 1, N0, li = 1, Ni (28)
is feasible for some τ i ∈ R, then the control in the form

ui(x) =    -∥x∥ ν i d bi(x) a(x)+ √ a 2 (x)+|B(x)| 4 |B(x)| 2 , if |B(x)| ̸ = 0, 0, if |B(x)| = 0, (29)
where

i = 1, m, ν i = deg F d (f 0 ) -deg F d (f i ), B(•) = (b 1 (•), ..., b m (•)) T , and
a(x) = 2p T (x)P f0(p(x)) p T (x)(P G d +G T d P )p(x) , b i (x) = 2p T (x)P fi(p(x)) p T (x)(P G d +G T d P )p(x) ,
globally finite-time (nearly fixed-time) stabilizes the system (1) at the origin if deg

F d (f 0 ) < 0 (deg F d (f 0 ) > 0). Proof. Assume that V = ∥x∥ d is a candidate control- Lyapunov function. Since ∥x∥ d = e s : ∥d(-s)x∥ = 1, ∂∥d(-s)x∥ ∂s = -∂∥z∥ ∂z z=d(-s)x G d d(-s)x = -1 2 ∥d(-s)x∥ -1 x T d T (-s) P G d + G T d P d(-s)x, ∂∥d(-s)x∥ ∂x = ∥d(-s)x∥ -1 x T d T (-s)P d(-s), then implying ∂s ∂x = -∂∥d(-s)x∥ ∂s -1 ∂∥d(-s)x∥ ∂x
(by means of Implicit Function Theorem [START_REF] Courant | Introduction to Calculus and Analysis[END_REF] for the function s : R n → R implicitly defined by ∥d(-s)x∥ = 1) we obtain

∂V (x) ∂x f 0 (x) = ∥x∥ 1+deg F d (f0) d a(x), ∂V (x) ∂x f i (x) = ∥x∥ 1+deg F d (fi) d b i (x), i = 1, m,
and for x ̸ = 0 we have

∂V (x) ∂x f 0 (x) + ∂V (x) ∂x m i=1 u i f i (x) = -∥x∥ 1+deg F d (f0) d a 2 (x) + |B(x)| 4
and the right-hand side of this inequality is negative provided that the condition

|B(x)| = 0 ⇒ a(x) < 0 ( 30
)
is satisfied, which according to (3) follows from the assumption that this V is a control-Lyapunov function.

Let us verify this hypothesis.

Since d is strictly monotone according to Theorem 3 we have P G d +G T d P > 0. Then due to the function f i can be represented as [START_REF] Pyatnitskiy | Criteria of asymptotic stability of differential inclusions and periodic motions of time-varying nonlinear control systems[END_REF], when we substitute f i from ( 27) into the expressions of a(x) and b i (x), then the condition [START_REF] Courant | Introduction to Calculus and Analysis[END_REF] takes the form

S ⊇ ∩ i Q i , i = 1, m, Q i = {x ∈ S : x T N
j=1 α i,j (x) P A i,j + A T i,j P x = 0}, S = {x ∈ S : x T N j=1 α 0,j (x) P A 0,j + A T 0,j P x < 0}. (31) Then applying S-procedure, the condition (31) is satisfied if the system of LMIs ( 28) is feasible for some τ i ∈ R, i = 1, m. Finally, since the controller (29) homogenizes the closed-loop system with degree deg F d (f 0 ), then according to Theorem 4 the closed-loop system is nearly fixed-time (finite-time) stable, if deg

F d (f 0 ) > 0 (deg F d (f 0 ) < 0). ■ Remark 5 If in Theorem 6 the strict inequality deg F d (f 0 ) > deg F d (f i
) is satisfied, then the control-Lyapunov function satisfies the small control property and the control ( 29) is continuous at 0.

Remark 6

The system of LMI ( 28) is just a sufficient condition to satisfy [START_REF] Courant | Introduction to Calculus and Analysis[END_REF]. In general, since

a(x) = 2p T (x)P f0(p(x)) p T (x)(P G d +G T d P )p(x) , b i (x) = 2p T (x)P fi(p(x)) p T (x)(P G d +G T d P )p(x) , i = 1, m,
the condition (30) can be numerically checked using a numerical grid on a sphere. Explicitly calculating a(x) and B(x) we can obtain that

∂V (x) ∂x f 0 (x) + ∂V (x) ∂x m i=1 u i f i (x) ≤ -CV 1+deg F d (f0) ,
where C = min y∈S a 2 (y)+|B(y)| 4 . Then, utilizing Theorem 1 we can get estimates for settling-time functions.

Example 6 Consider the system [START_REF] Polyakov | Generalized Homogeneity in Systems and Control[END_REF]

with x ∈ R 2 , m = 2, where f 0 = -|x1| 1/2 x 2 2 0 , f 1 = x 2 2 +x1 x2 , f 2 = x 3 2 x1+x 2 2
. In this system f 0 , f 1 and f 2 are d-homogeneous of degrees 1, 0 and 1, respectively, where d is strictly monotone dilation equipped with the norm ∥x∥ = √

x T x, and G d = ( 2 0 0 1 ) is the generator of the dilation. For these functions, using Proposition 1 we have representations in the form [START_REF] Pyatnitskiy | Criteria of asymptotic stability of differential inclusions and periodic motions of time-varying nonlinear control systems[END_REF] with the matrices A 0i = {-0.438,0.438} {-0.437,0.437} 0 0

, A 1j = 1 {-1,1} 0 1
,

A 2l = {0,2} {-0.333,0.333} 1 {-1,1}
, where i = 1, 4, j = 1, 2, l = 1, 8. For these matrices the system ( 28) is feasible with P = I 2 , τ 1 = 13.8175, τ 2 = 0.9915, i.e., the closedloop system with the control ( 29) is nearly fixed-time stable. Numerical simulations (Fig. 1) have been done using explicit Euler method with the step h = 10 -3 and

x 0 = 5 -1 T
. The results of simulation for five different initial conditions with using the logarithmic scale are shown in Fig. 2 in order to demonstrate nearly fixedtime convergence rate (e.g., in accordance with Theorem 1 and Remark 6 the set B = {x ∈ R 2 : V (x) ≤ 1} is fixed-time attractive with T max = 4).

where x ∈ R n , u ∈ R m , B ∈ R n×m , and f can be represented in the form

f (x) = ∥x∥ ν d N i=1 α i (p(x))d(ln ∥x∥ d )A i p(x), (33)
with corresponding matrices

A i ∈ R n×n , i = 1, N , a generator G d ∈ R n×n , ν ∈ R,
where the dilation d is strictly monotone on R n equipped with the norm ∥x∥ P = √ x T P x, for P > 0 satisfying [START_REF] Andrieu | Homogeneous approximation, recursive observer design, and output feedback[END_REF].

Theorem 7 Let G d B = γB, where γ ∈ R : γ + ν ≥ 0. Let K ∈ R m×n , H ∈ R n×n , η, δ ∈ R + be a solution of the bilinear matrix inequalities (BMIs) ηP ≥ H > 0, θmin(P G d +G T d P ) ≤ H -HG d -G T d H ≤ θmax(P G d +G T d P ), (θmaxP +H) (Ai +BK)+(Ai +BK) T (θmaxP +H) ≤ -δP, (θminP +H) (Ai +BK)+(Ai +BK) T (θminP +H) ≤ -δP, (34)
for i = 1, N and some θ min , θ max ∈ R, then the system [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] with

u(x) = ∥x∥ γ+ν d Kp(x) (35)
is finite-time (nearly fixed-time) stable if ν < 0 (ν > 0) .

Proof. Let us show that the equality

d(-ln ∥x∥ d )∥x∥ ν+γ d BKp(x) = ∥x∥ ν d BKp(x) (36)
holds. To do this, it is enough to show that (e γs d(-s) -I n )BK = 0.

Since (e γs d(-s) -

I n )BK = +∞ i=1 s i (γIn-G d ) i i!
BK, the equality (37) holds if (G d -γI n )BK = 0 which is implied by G d B = γB. Then, using (37) it is easy to show that (36) holds.

Let ν ≥ -1. Choose a candidate Lyapunov function in the form [START_REF] Bhat | Finite-time stability of continuous autonomous systems[END_REF]. Due to (36) we have that

d(-ln ∥x∥ d ) ẋ = ∥x∥ ν d N i=1 αi(p(x))Ai + BK p(x). (38)
Then following the proof of Lemma 4 and applying (38) we have

V = 2∥x∥ 1+ν d p T (x) p T (x)(H-HG d -G T d H)p(x) p T (x)(P G d +G T d P )p(x) P + H × N i=1 α i (p(x))A i + BK p(x) = 2∥x∥ 1+ν d p T (x) (θ(p(x))P + H) × N i=1 α i (p(x))A i + BK p(x), where p T (x)(H-HG d -G T d H)p(x) p T (x)(P G d +G T d P )p(x) = θ(p(x)) ∈ [θ min , θ max ].
Using convex embedding one can show that Ṽ µ+ν µ . Finally, applying Theorem 1 one can show that the feedback law (35) stabilizes the system (32) with different types of convergence. ■ Remark 7 For ν < 0 the closed-loop system (32), ( 35) is finite-time stable with the following settling time estimate T (x 0 ) ≤ -

η 1+ν V -ν 0 δν
, where V 0 = V (x 0 ). For ν > 0 the closed-loop system is nearly fixed-time stable and for any ε

∈ R + , the set B = {x ∈ R n : V (x) ≤ ε} is fixed-time attractive with T max = η 1+ν δν ε -ν .
Remark 8 Analogously to Remark 4 and Algorithm 1, feasibility of BMI (34) may depend on the choice of the matrix P , and the following LMI-based heuristic algorithm can be applied:

Algorithm 2 INITIALIZATION: take P = I n , K = (0 0 0); compute matrices A i ; STEP: While LMI (34) is not feasible with respect to H ∈ R n×n , θ min , θ max ∈ R, η, δ ∈ R + Solve the LMI system A i X +XA T i +BY +Y T B T +δX < 0, i = 1, N , G d X + XG T d > 0, X > 0 (41) for some δ ∈ R + with respect to X ∈ R n×n , Y ∈ R m×n .
take P = X -1 , K = Y P ; recompute matrices A i ; end Note that if Algorithm 2 (Algorithm 1) does not find a solution, it is advisable to reduce δ ∈ R + . The logic of these algorithms can be augmented by the maximum admissible number of iterations with stop afterwards avoiding an infinite loop.

Example 7 Consider the system (32) with f (x) as in Example 4 and B = 1 0 0 T . The condition G d B = γB is satisfied for γ = 3. Applying Algorithm 2 for η = 0.2, δ = 0.1 we obtain K = -5.5055 -15.8387 -16.3807 , and the closed-loop system is nearly fixed-time stable due to ν = 1. Numerical simulations (Fig. 345) have been done using explicit Euler method with the step h = 10 -3 . Due to homogeneity the closed-loop system is ISS with respect to measurement noise and additive disturbances. Fig. 4 shows the results of simulation with measurement band limited noise of power 10 -6 and additive disturbances d = 0.5 sin(x 1 ) cos(x 2 ) 0 0 T . The results of simulation for three different initial conditions using the logarithmic scale of V are shown in Fig. 5 in order to demonstrate nearly fixed-time convergence rate (e.g., in accordance with Remark 7 the set B = {x ∈ R 3 : V (x) ≤ 1} is fixed-time attractive with T max = 0.4). Note, that the given control is a generalized version of [START_REF] Zimenko | Robust Feedback Stabilization of Linear MIMO Systems Using Generalized Homogenization[END_REF]. Indeed, if f (x) in ( 32) is linear, then the proposed result coincides with [START_REF] Zimenko | Robust Feedback Stabilization of Linear MIMO Systems Using Generalized Homogenization[END_REF].

Remark 9 Note that the function ∥x∥ d is defined implicitly by [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF]. In order to realize the control [START_REF] Sontag | A 'universal' construction of Artstein's theorem on nonlinear stabilization[END_REF] or (35) the bisection method may be utilized: Here ϵ ∈ R + is a minimal admissible value of ∥x i ∥ d . If STEP is applied recurrently many times to the same vector x i then it allows to localize the unique positive root of the equation ∥d(-ln ∥x i ∥ d )x∥ = 1 as follows a < ∥x i ∥ d < b, and next, to improve the latter estimate using the bisection method.

Conclusions

The paper presents methods for representation of homogeneous systems in the canonical form [START_REF] Grüne | Homogeneous state feedback stabilization of homogeneous systems[END_REF]. These methods were extended for some non-homogeneous systems as well. It is demonstrated that the representation in the canonical form allows to present stability conditions in the form of LMIs even for sufficiently nonlinear systems. In the paper, the stabilizing control algorithms for affine in control nonlinear systems are presented utilizing the representation in the canonical form. Moreover, due to homogeneity property the proposed controls are robust (in ISS sense) to a large class of disturbances and measurement noises. Numerical examples demonstrate effectiveness of the proposed control algorithms.

x2x3 x3 0

 0 that is d-homogeneous with negative degreedeg F d = -1 and G d = ∥x∥ = √x T x one can obtain the representation (

Example 3

 3 ) stays true with N ≤ 2 n 2 , N i=1 α i (x) = 1 and A i = C + Ãi , where Ãi = (νI n + G d ) -1 Mi G d ( Ãi = Mi P G d ) and Mi are vertices obtained from representation of f 2 based on Proposition 1 (Proposition 2). Consider the system ẋ = x1-x2+|x2| 1/2 x1-x2+|x2| 1/2

,

  and for b 1 (x) = cos(x 1 ), b 2 (x) = b 3 (x) = 1 the condition (20) is satisfied. Then according to Proposition 1 for ∥x∥ = √ x T P x, P =

2 Fig. 1 .Fig. 2 . 5 . 3

 21253 Fig. 1. System states x versus time

(

  θ(p(x))P + H) N i=1 αi(p(x))Ai + BK + N i=1 αi(p(x))Ai +BK T (θ(p(x))P +H) ≤ -δP, (39)if the inequalities(θmaxP +H) (Ai +BK)+(Ai +BK) T (θmaxP +H) ≤ -δP, (θminP +H) (Ai +BK)+(Ai +BK) T (θminP +H) ≤ -δP, are satisfied for i = 1, N . If (39) is satisfied then we have V ≤ -δ∥x∥ 1+ν d ≤ -δ η 1+ν V 1+ν . (40)Similarly to the proof of Lemma 4, in the case ν < -1 choose the candidate Lyapunov function in the form Ṽ = V µ , where µ ≥ -ν. Then according to (40) we can obtain V ≤ -µδ η 1+ν

3 Fig. 3 . 3 Fig. 4 .Fig. 5 .

 33345 Fig. 3. System states x versus time

Algorithm 3

 3 INITIALIZATION: ∥x∥ d = 1; a = ϵ min ; b = 1; STEP : If x T i d(ln b) T P d(ln b)x i > 1 then a = b; b = 2b; elseif x T i d(-ln a) T P d(ln a)x i < 1 then b = a; a = max a 2 , ϵ min ; else c = a+b 2 ; If x T i d(ln c) T P d(ln c)x i < 1 then b = c; else a = max{ϵ min , c}; endif; endif; ∥x i ∥ d = b;