Peter Markowich
email: peter.markowich@kaust.edu.sa

Benoît Perthame
email: benoit.perthame@sorbonne-universite.fr

Di Fang

Jin ‡ Shi
email: shijin-m@sjtu.edu.cn

Implicit and Semi-implicit Numerical Schemes for the Gradient Flow of the Formation of Biological Transport Networks *

Implicit and Semi-implicit Numerical Schemes for the Gradient Flow of the Formation of Biological Transport

Networks

Introduction

Biological transport networks have drawn extensive research interests due to their ubiquitous existence in living organisms and rich phenomenon observed in nature, two common examples of which are the leaf venations and blood flow. It has been hypothesized in the literature that the branching structures of the network is governed by an optimization of the energy consumption of the living system, as a consequence of natural selection (see for example [START_REF] Nelson | Leaf vascular pattern formation[END_REF][START_REF] Laguna | The role of elastic stresses on leaf venation morphogenesis[END_REF][START_REF] Katifori | Damage and fluctuations induce loops in optimal transport networks[END_REF][START_REF] Corson | Fluctuations and redundancy in optimal transport networks[END_REF][START_REF] Bohn | Structure, scaling, and phase transition in the optimal transport network[END_REF][START_REF] Dimitrov | A constant production hypothesis guides leaf venation patterning[END_REF]). However, the internal mechanism may appear rather counter-intuitive from a biological point of view, simply due to the global nature of the optimization approach in general. One naturally questions how one blood vessel, for example, with only the "knowledge" of local information, is able to participate in the global energy optimization task of the entire network. The link in between was made transparent by Hu and Cai in 2013 [START_REF] Hu | Adaptation and optimization of biological transport networks[END_REF]. In their work, a global energy functional approach is suggested considering both material and metabolism costs, but nevertheless the corresponding gradient flow ends up driven by the wall shear stress on the tube walls, which is a local information and has been observed in experiments that can be sensed by the tissue [START_REF] Pohl | Crucial role of endothelium in the vasodilator response to increased flow in vivo[END_REF][START_REF] Koller | Endothelial regulation of wall shear stress and blood flow in skeletal muscle microcirculation[END_REF][START_REF] Lu | Role of shear stress and stretch in vascular mechanobiology[END_REF]. This model unifies the global picture with the local one, and is referred to later as the Cai-Hu model. The original Cai-Hu model is in ODE form, which unfortunately cannot describe the growth or initiation of a network. A generalized PDE version is henceforth proposed and studied in [START_REF] Haskovec | Mathematical analysis of a pde system for biological network formation[END_REF][START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF][START_REF] Albi | Biological transportation networks: Modeling and simulation[END_REF][START_REF] Haskovec | ODE and PDE based modeling of biological transportation networks[END_REF][START_REF] Burger | A mesoscopic model of biological transportation networks[END_REF][START_REF] Haskovec | Rigorous Continuum Limit for the Discrete Network Formation Problem[END_REF]. To aim at the understanding of the growth and formation of these biological networks, we consider the PDE version in this work, which is in fact a lot more challenging in terms of numerical simulations.

In the PDE model, one has a coupled system in terms of the conductance vector m(t, x) ∈ R d (with space dimension d ≤ 2) whose direction coincides with the allowable flow direction and whose modulus indicates how strongly the tissue passes the flow at position x, and the pressure p(t, x) ∈ R as follows:

-∇ • [(r(x)I + m ⊗ m)∇p] = S, (1.1)

∂m ∂t = D 2 ∆m + c 2 (∇p ⊗ ∇p)m -α|m| 2(γ-1) m, (1.2)
where r(x) denotes the background permeability, I is the identity matrix, S(x) is the distribution of fluid sources, D > 0 the diffusive constant (usually very small), c > 0 the activation constant driving the network adaptation, α > 0 is the metabolic constant, and γ ∈ [1/2, 1] is the metabolic rate of the bio-organism according to Murray's Law [START_REF] Murray | The Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume[END_REF]. To be specific, it has been shown that for blood vessels, γ = 1/2 while γ ∈ (1/2, 1] for leaf venation [START_REF] Hu | Adaptation and optimization of biological transport networks[END_REF]. Notice that choosing the time scale allows to fix a parameter (for instance α = 1 as later). Also note that the relationship between the pressure p(t, x) and the conductivity vector m(t, x) follows Darcy's Law in porous media written as the Poisson equation (1.1), where r(x)I + m ⊗ m is the permeability tensor. The second equation (1.2) is a reaction diffusion equation derived as a gradient flow by minimizing the energy cost functional consumed by the system. Now let us make explicit the total energy functional

(1.3) E(m) := 1 2 Ω (D 2 |∇m| 2 + α γ |m| 2γ + c 2 |m • ∇p[m]| 2 + c 2 r(x)|∇p[m]| 2 ) dx,
where p[m] is the unique solution of the Poisson equation (1.1) subject to appropriate boundary conditions, the first term represents the diffusive energy, the second the metabolic cost consumed by the network (according to the celebrated Murray's Law in mathematical biology) to keep the tissue alive, and the last two terms correspond to the energy consumption incurred by the fluid itself -the former for the network while latter for the background. Again the pressure p(t, x) here is given according to Kirchoff's law (1.1), where the flux follows Darcy's Law. We point out that this is clearly a non-convex optimization problem in general, where the first two terms are the convex part in the physical cases γ ≥ 1/2 while the coupling with the Poisson equation gives rise to the non-convexity. We remark that the reaction term in (1.2) illustrates a competition between the contribution by activation forces (by flow inside the network) c 2 (∇p ⊗ ∇p)m and the metabolic part -α|m| 2(γ-1) m. And it is exactly this competition that gives rise to interesting branching phenomena and nontrivial equilibria.

Although model (1.1)-(1.2) may appear "standard" at a first glance in the form of a Poisson equation coupled with an diffusion-reaction equation, it is by no means easy to compute in practice. The major numerical difficulty of this system lies in the stiffness of basically all terms in (1.2). Unlike the rather standard reaction-diffusion equation (such as the Allen-Cahn equation) where the reaction term is smooth, we point out that with the physical parameter γ ∈ (1/2, 1), the activation term α|m|2(γ-1) m in (1.2) is not smooth and in fact stiff since m can be close to zero (to be made more transparent in later discussions with the 1D case). This feature makes our problem more difficult than well-studied reaction diffusion equations where the smoothness of the reaction term is guaranteed. Hence a direct application of the currently-developed tricks for reaction diffusion equations, such as linear penalization [START_REF] Du | Stability analysis and application of the exponential time differencing schemes[END_REF][START_REF] Du | Analysis and applications of the exponential time differencing schemes and their contour integration modifications[END_REF][START_REF] Feng | Nonlinear stability of the implicit-explicit methods for the allen-cahn equation[END_REF] and its resulting Integration Factor (IF) method or Exponential Time Differencing (ETD) method [START_REF] Du | Stability analysis and application of the exponential time differencing schemes[END_REF][START_REF] Du | Analysis and applications of the exponential time differencing schemes and their contour integration modifications[END_REF] remain not clear. Before proceeding any further we summarize the main numerical difficulties of this system as follows 1. The diffusion term D 2 ∆m is stiff as in any reaction-diffusion system. This is actually not a real difficulty at all, since the stiffness of the diffusion can be removed by treating the term exactly or implicitly. Different approximations on the temporal integral involving reaction term in this approach result in either Integration Factor (IF) methods or Exponential Time Differencing (ETD) methods [START_REF] Du | Stability analysis and application of the exponential time differencing schemes[END_REF][START_REF] Du | Analysis and applications of the exponential time differencing schemes and their contour integration modifications[END_REF].

the physical energy in 1D for the semi-implicit solver, numerically one observes that the physical energy decays in time for both 1D and 2D cases. The rest of the paper is organized as follows: In Section 2, we briefly revisit the 1D Cai-Hu model, and propose an implicit algorithm whose unconditional stability is proved. Motivated by the practical calculation for higher dimensional cases, we henceforth propose a semi-implicit scheme, and prove that it decays the physical energy provided the stability condition is satisfied in 1D. Section 3 presents both the fully implicit and semi-implicit schemes for 2D in analogy to the 1D case. In Section 4, numerical tests are conducted in 1D with both the implicit and semi-implicit solvers, and 2D with the semi-implicit one. We observe numerically in both cases the energy decreases in time. The 2D examples present branching phenomena as occurred in leaf venations in nature. Section 5 contains conclusions.

2 The 1D Cai-Hu Model and Numerical Approximations

The 1D Cai-Hu Model

As has been shown in [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF], the 1D Cai-Hu model has the form

(2.1) ∂ t m -D 2 ∂ 2 xx m = c 2 B(x) 2 (1 + m 2 ) 2 -|m| 2(γ-1) m
with proper boundary conditions (see [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF] for the detailed derivation), where the parameters are already chosen as r(x) = α = 1 for simplicity and B(x) is the given function

B(x) := x 0 S(y)dy.
The energy in the 1D setting is then simplified as

(2.2) E(m) = 1 2 Ω D 2 |∂ x m| 2 + 1 γ |m| 2γ + c 2 B(x) 2 1 + m 2 dx.
As has been pointed out in the introduction, the nonlinear reaction term is also stiff. To see it, simply defines

(2.3) H(m) = c 2 B(x) 2 (1 + m 2 ) 2 -|m| 2(γ-1) m,
where -1 ≤ 2(γ -1) ≤ 0, since 1 2 ≤ γ ≤ 1, and

(2.4) H (m) =    c 2 B(x) 2 1-3m 2 (1+m 2 ) 3 -2γ-1 |m| 2(1-γ) , if m = 0, -∞, if m = 0.
Note that |H (m)| is very large as m approaches zero.

An Implicit Treatment and Energy Decay

Consider the 1d setting in [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF], namely, the interval (0, 1) and homogenous Neumann Boundary Conditions for m. We propose the following scheme:

m n+1 j -m n j ∆t = D 2 2∆x 2 m n+1 j+1 -2m n+1 j + m n+1 j-1 + m n j+1 -2m n j + m n j-1 (2.5) + c 2 B 2 j (m n+1 j + m n j ) 2(1 + (m n+1 j ) 2 )(1 + (m n j ) 2 ) - (m n+1 j ) 2γ -(m n j ) 2γ 2γ(m n+1 j -m n j ) . (2.6) m -1 = m 1 , m N -1 = m N +1 .
Note that this scheme is second order in time, and the diffusion term is treated by the Crank-Nicolson method. We define the discrete energy as (2.7)

E n = 1 2 M j=1 D 2 ∆x 2 m n j -m n j-1 2 + 1 γ |m n j | 2γ + c 2 B 2 j 1 + (m n j ) 2 ,
which is a discrete analog of the energy functional (2.2). We shall prove in the following theorem that our scheme decays the discrete energy and is hence unconditionally stable.

Theorem 2.1. The scheme (2.6) satisfies the discrete energy dissipation equation

(2.8) E n+1 -E n = - 1 ∆t M j=1 (m n+1 j -m n j ) 2 ,
and hence decays the energy (2.9)

E n+1 ≤ E n .
Remark 2.2. This result implies that the scheme is unconditionally stable, i.e. there is no restriction on ∆t. It shows the scheme preserves in the discrete fashion the energy equation given in Lemma 1 of [START_REF] Haskovec | Mathematical analysis of a pde system for biological network formation[END_REF] (2.10)

d dt E = - Ω ∂m ∂t 2 dx.
In fact, one should simply view the scheme as a discretization of the energy equation.

Proof. The lemma follows a straightforward energy estimate, namely, multiply m n+1 j -m n j on both sides of the scheme (2.6) and sum over j,

M j=1 (m n+1 j -m n j ) 2 ∆t = (I) + (II) + (III),
where

(I) = M j=1 D 2 2∆x 2 m n+1 j+1 -2m n+1 j + m n+1 j-1 + m n j+1 -2m n j + m n j-1 (m n+1 j -m n j ), (II) = M j=1 c 2 B 2 j 2 (m n+1 j ) 2 -(m n j ) 2 (1 + (m n+1 j ) 2 )(1 + (m n j ) 2 ) , (III) = - M j=1 1 2γ (m n+1 j ) 2γ -(m n j ) 2γ .
It follows from summation by parts that

j m n+1 j+1 -2m n+1 j + m n+1 j-1 m n+1 j = j m n+1 j+1 -m n+1 j m n+1 j - j m n+1 j -m n+1 j-1 m n+1 j (2.11) = - j m n+1 j -m n+1 j-1 2 . (2.12)
Similarly applying summation by parts to other terms, one finds

(I) = D 2 2∆x 2 M j=1 m n+1 j+1 -2m n+1 j + m n+1 j-1 m n+1 j - D 2 2∆x 2 M j=1 m n j+1 -2m n j + m n j-1 m n j + D 2 2∆x 2 M j=1 m n j+1 -2m n j + m n j-1 m n+1 j - D 2 2∆x 2 M j=1 m n+1 j+1 -2m n+1 j + m n+1 j-1 m n j = - D 2 2∆x 2 j m n+1 j -m n+1 j-1 2 + D 2 2∆x 2 j m n j -m n j-1 2 .
A simple calculation shows

(II) = M j=1 c 2 B 2 j 2 (m n+1 j ) 2 -(m n j ) 2 (1 + (m n+1 j ) 2 )(1 + (m n j ) 2 ) = - 1 2 M j=1 c 2 B 2 j 1 + (m n+1 j ) 2 + 1 2 M j=1 c 2 B 2 j 1 + (m n j ) 2 .
Hence,

(I) + (II) + (III) = -E n+1 + E n ,
which completes the proof.

In numerical implementations, we use the Newton method to solve the implicit algebraic equation (2.6) in every time step. Note that supposing D = 0, then the algebraic equation has only one root, and sufficiently away from 0 the function in consideration is actually monotone. Numerically, the Newton iteration converges just in a few steps taking the initial guess sufficiently away from 0.

A Semi-implicit Treatment and Energy Decay

To avoid the nonlinear Newton solver, which may become impractical for higher dimension, we propose the following semi-implicit treatment in preparation for the 2D practical simulations.

(2.13)

m n+1 j -m n j ∆t - D 2 ∆x 2 m n+1 j+1 -2m n+1 j + m n+1 j-1 = c 2 B(x) 2 1 + (m n j ) 2 2 m n j -|m n j | 2(γ-1) m n+1 j , m -1 = m 1 , m N -1 = m N +1 .
Note that in (2.13), the activation part on the right hand side is treated fully explicitly, while the metabolic term is treated in a semi-implicit fashion. Since |m n j | 2(γ-1) is strictly positive, this term plays the role of stabilizing the scheme. Moreover, one only needs to invert a tridiagonal matrix where standard fast algorithms can be used. In the following theorem, we establish the conditional stability of the semi-implicit solver in the semi-discrete set-up which guarantees the physical energy does not increase. Note that the semi-discretization is used for a clear presentation, the fully-discrete stability should be similar, which is omitted here. Theorem 2.3. Assume B ∈ L ∞ , consider the semi-implicit scheme (2.13) in the semi-discrete set-up, i.e.,

(2.14) m n+1 -m n ∆t -D 2 ∂ xx m n+1 = c 2 B(x) 2 (1 + (m n ) 2 ) 2 m n -|m n | 2(γ-1) m n+1 ,
then the physcial energy as defined in (2.2) decays

E(m n+1 ) ≤ E(m n ),
provided the condition

(2.15) c 2 B 2 ∞ ∆t ≤ 2, is satisfied. Proof. Denote f (m) = c 2 B(x) 2 (1 + m 2 ) 2 m, F (m) = c 2 B(x) 2 2(1 + m 2 ) , then clearly one has F (m) = -f (m) and (2.16) |f (m)| = c 2 B(x) 2 1 -3m 2 (1 + m 2 ) 3 ≤ c 2 B 2 ∞ .
In the following, we use •, • to denote the inner product and • for the L 2 norm in Ω.

On one hand, multiplying the scheme (2.14) by m n+1 -m n and integrating in x, one obtains (2.17)

1 ∆t m n+1 -m n 2 + D 2 2 ∂ x m n+1 2 - D 2 2 ∂ x m n 2 + D 2 2 ∂ x m n+1 -∂ x m n 2 = f (m n ), m n+1 -m n - Ω |m n | 2γ-2 (m n+1 ) 2 -m n+1 m n dx,
where integration by parts and the equality a -b, a = 1 2 a 2 -1 2 b 2 + 1 2 a -b 2 are used. On the other hand, it holds that

(2.18) E(m n+1 ) -E(m n ) = D 2 2 ∂ x m n+1 2 - D 2 2 ∂ x m n 2 + Ω 1 2γ |m n+1 | 2γ - 1 2γ |m n | 2γ dx + Ω F (m n+1 ) -F (m n )dx,
by the definition of the physical energy (2.2). Notice that

F (m n+1 ) -F (m n ) = -f (m n )(m n+1 -m n ) - m n+1 m n f (τ )(m n+1 -τ )dτ,
and hence combining with (2.17), we have

E(m n+1 ) -E(m n ) + 1 ∆t m n+1 -m n 2 + D 2 2 ∂ x m n+1 -∂ x m n 2 = - Ω |m n | 2γ-2 (m n+1 ) 2 -m n+1 m n - 1 2γ |m n+1 | 2γ + 1 2γ |m n | 2γ dx- Ω m n+1 m n f (τ )(m n+1 -τ )dτ dx := I 1 + I 2
By (2.16), clearly

I 2 ≤ 1 2 c 2 B 2 ∞ m n+1 -m n 2 .
Next, we shall show that I 1 ≤ 0. Case 1: If m n+1 m n ≥ 0, it suffices to show that the function

(2.19) G(x, y) = x 2γ-2 y 2 -x 2γ-1 y - 1 2γ y 2γ + 1 2γ
x 2γ ≥ 0 for all x, y ≥ 0. (Note that I 1 corresponds to the case when x = |m n | and y = |m n+1 |.) Though one could prove the result by proving the minimum of G(x, y) is no less than 0, an easier way is to use the homogeneity of the inequality and convert the problem into a single variable function.

To be specific, note that

(2.20) G(x, y) = x 2γ ( y x ) 2 - y x - 1 2γ ( y x ) 2γ + 1 2γ .
Denote z = y/x ≥ 0, it suffices to show that

A(z) = z 2 -z - 1 2γ z 2γ + 1 2γ ≥ 0
for all z ≥ 0. It clearly holds for z = 0 and for γ = 1/2 respectively, hence we only need to consider z > 0 and γ ∈ (1/2, 1). Straightforward computation gives

A (z) = 2z -1 -z 2γ-1 , A (z) = 2 - 2γ -1 z 2-2γ ,
where A (z) is monotone increasing with a unique zero at z = ( 2 2γ-1 ) 2γ-2 := z 0 , which means A (z) decreases in (0, z 0 ) and increases in (z 0 , ∞). On one hand, one has A (0) = -1 < 0, which implies that A (z) has at most one zero. On the other hand, z = 1 is clearly a critical point, and therefore is the unique critical point of A(z). One could easily check the second derivatives and the limits on the boundary to see that the critical point z = 1 is the global minimum. It follows that

A(z) ≥ A(1) = 0, G(x, y) ≥ 0. Case 2: When m n+1 m n < 0, it suffices to show that G(x, y) = x 2γ-2 y 2 + x 2γ-1 y - 1 2γ y 2γ + 1 2γ x 2γ ≥ 0 for all x, y ≥ 0, which clearly holds since G(x, y) = G(x, y) + 2x 2γ-1 y ≥ 0.
In summary, we have

E(m n+1 ) -E(m n ) + D 2 2 ∂ x m n+1 -∂ x m n 2 ≤ ( c 2 B 2 ∞ 2 - 1 ∆t ) m n+1 -m n 2 ≤ 0,
where the stability condition (2.15) is used. This completes the proof.

The 2D Cai-Hu Model and Numerical Schemes

The ideas developed in the 1D case can be extended to treat the spatially discretized 2D Cai-Hu model when set on a rectangular grid. For completeness, we present both the fully implicit and implicit-explicit schemes.

An Implicit Approach and Energy Decay

Similar as the 1D case, an implicit scheme can be constructed as

-∇ h • r 2 ∇ h p n+1 + ∇ h p n + 1 4 m n+1 • ∇ h p n+1 + m n • ∇ h p n m n+1 + m n = S, (3.1) m n+1 l -m n l ∆t = D 2 2 ∇ 2 h m n+1 l + m n l + c 2 4 m n+1 • ∇ h p n+1 + m n • ∇ h p n ∇ h p n+1 l + ∇ h p n l (3.2) - α 2γ |m n+1 | 2(γ-1) (m n+1 l ) 2 -|m n | 2(γ-1) (m n l ) 2 m n+1 l -m n l ,
where the subscript l = 1, 2 labels the l-th component of the vector m or ∇ h p; the spatial discretizations are all performed at the grid point with index-pair (k, j) and hence spatial indices are omitted for simplicity; and ∇ 2 h denotes the discretized Laplacian operator, and hence the diffusion term is treated in a Crank-Nicolson matter (similarly as 1D case). This scheme is clearly second order in time.

However, this scheme gives rise to an nonlinear algebra problem, which in practice is very difficult to solve. Note that a direct application of the Newton iteration may fail to converge depending on initial guesses. Hence in the next section, we propose a more practical numerical method.

A Semi-implicit Scheme

The key goal here is to avoid a nonlinear algebraic problem which requires the Newton iterations, in which case hopefully one only deals with a linear algebra problem that provides enough efficiency in practical calculation. Following the semi-implicit idea for the 1D Cai-Hu model, we propose the following scheme for multi-dimensional cases:

-D 1 • [(rI + m n ⊗ m n )D 1 p n+1 ] = S, (3.3) m n+1 1 -m n 1 ∆t = D 2 D 2 m n+1 1 + c 2 (D x p n+1 ) 2 m n 1 + D x p n+1 D y p n+1 m n 2 -α|m n | 2(γ-1) m n+1 1 , (3.4)
m n+1 2 -m n 2 ∆t = D 2 D 2 m n+1 2 + c 2 (D y p n+1 ) 2 m n 2 + D x p n+1 D y p n+1 m n 1 -α|m n | 2(γ-1) m n+1 2 ,
where D x and D y denote the spatial discretizations in x and y directions, respectively, which could be chosen as either the forward difference or the central difference, D 1 = (D x , D y ), and

D 2 = D 2 x + D 2 y .
Note that similarly to the 1D scheme (2.13), the diffusion term is treated implicitly, the metabolism term semi-implicitly and the activation term explicitly in the sense that the pressure p n+1 used here is computed from the Poisson equation (3.3) using the information of m n . This semi-implict scheme is a 2D version of scheme (2.13) that has been proven to decay the physical energy in 1D, and is efficient in actual computations, since the semi-implicit treatment of the last term in (3.4) helps to stabilize the whole procedure and create positive definite matrices to be inverted. For numerical tests, we use the staggered grids, as is widely used in computational fluid dynamics to avoid the famous "checkerboard" problem.

Numerical Results

We now illustrate the theoretical results with numerical simulations. We compare the semiimplicit scheme and the fully implicit which is more complex due to Newton iterations but avoids the restriction on time steps. We pay a special attention to the decay of energy as it holds for the continuous model, and on the ability of the scheme to generate patterns.

The 1D Cai-Hu Model

In this section, we present our numerical results of the 1D Cai-Hu model (2.1) using the implicit solver (2.6) and the semi-implicit scheme (2.13), respectively. Denote the nonlinear system of equations resulting from the implicit scheme as

(4.1) F (m n+1 ) := m n+1 j -m n j ∆t - D 2 2∆x 2 m n+1 j+1 -2m n+1 j + m n+1 j-1 + m n j+1 -2m n j + m n j-1 - c 2 B 2 j (m n+1 j + m n j ) 2(1 + (m n+1 j ) 2 )(1 + (m n j ) 2 ) + (m n+1 j ) 2γ -(m n j ) 2γ 2γ(m n+1 j -m n j ) = 0,
where F : R N → R N . The Newton method reads

J F m n+1 (k) m n+1 (k + 1) -m n+1 (k) = -F m n+1 (k) ,
with J F (m n+1 ) as the N × N Jacobian matrix of F and k as the k-th step of iteration. In order to solve this linear system, the Thomas algorithm [START_REF] Conte | Elementary Numerical Analysis: An Algorithmic Approach[END_REF] is implemented since the Jacobian is a tri-diagonal matrix.

Example 4.1 (1D Extinction of Solutions). It has been shown in Lemma 4 of [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF] that for 1/2 ≤ γ ≤ 1, the initial data m 0 ∈ L ∞ (0, 1) and the parameters satisfying

c B(x) L ∞ ≤ Z γ := 2 γ + 1 1 -γ 1 + γ γ-1 2
, one has the extinction of solutions, namely, the solution converges to zero as t → ∞. In other words, the solution converges to zero in infinite time for 1/2 ≤ γ ≤ 1 (in contrast to the case -1 ≤ γ < 1/2 where there exists a finite break down time). In this example, we try to observe the vanishing effect numerically. Consider γ = 0.75 corresponding to the leaf venation case, D = 0.01, cB(x) = Z γ /2, and initial datum as (4.2) m(0, x) = 10(sin(2πx) + 1) cos(2πx) -3, which is plotted in Fig. 1a. We choose this initial datum so that it is a bit more complicated than simple algebraic functions, and here one has different combinations of monotonicity and convexity in the sub-intervals. Both the implict and semi-implict solvers are implemented with ∆t = 0.01 and ∆x = 0.005 until the final time T = 10. We first show the results of the implicit solver: Fig 1 plots the solutions at various time t = 0, 2, 5, 10 and the discrete physical energy defined in (2.7). It can be seen that the discrete physical energy decays throughout the simulation in time as is proved in Theorem 2.1. Next, the solutions computed by the semi-implicit scheme at different time and its discrete physical energy are plotted in Fig. 2. In the numerics one indeed can see the decay of the discrete physical energy as is proved in Theorem 2.3. In both methods, the solution converges finally to zero and we numerically recover the theoretical results proved in Lemma 4 of [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF].

The 2D Cai-Hu Model

In this section, we present a number of numerical examples for the 2D Cai-Hu model, using the semi-implicit solver. The linear algebraic problems resulting from the scheme are solved using a preconditioned conjugate gradient method. The discrete physical energy considered here is defined as E(m n ) at each time t = t n , with E defined in (1.3). We choose γ = 0.75, which corresponds to leaf venations according to [START_REF] Hu | Adaptation and optimization of biological transport networks[END_REF]. The initial data and source term are given in Fig. 5.

In this example we use ∆t = 0.005, and ∆x = ∆y = 0.002. We plot the long-time behavior of the system for various time T = 20, 40, 60, 80. To be specific, Fig. 6 shows |m| 2 in the log-scale, which corresponds to the structure of the leaf. Note that for m = 0, the log function is not defined, and hence hereafter log(|m| 2 + 0 ) is plotted in numerics, where 0 = 4e-16 is of size of machine precision. In this example, we also plot the volumetric flux u (proportional to the flow velocity in leaf venations), defined as

(4.3) u = (rI + m ⊗ m)∇p
in Fig. 7. It can be seen that the leaf venation grows as the fluid spreads, and finally achieve a tree-like structure. Numerically, we see that the discrete physical energy decays with respect to time through out the simulation as is shown 8. In fact, in the last several steps, the energy keeps the decay behavior even if the first five decimal places in the energy value no longer change. To be specific, the last seven energy values are 2.40638937, 2.40638876, 2.40638813, 2.40638744, 2.40638683, 2.40638619, 2.40638551.

Example 4.4 (Varying D). In this example, we investigate how the diffusion coefficients may affect the formation of the biological network. Consider the same set-up as in Example 4.3, but with different diffusion coefficients D = 0.5, 0.1, 0.02, 0.005, 0.002. The solutions are computed till time T = 40. It can be seen from Fig. 9 that the smaller the diffusion coefficient is, the more detailed structure the network tends to exhibit and also the faster it reaches the equilibrium. Example 4.5 (Varying spatial mesh size). In this example, we test the semi-implicit scheme with different spatial mesh size. In this example, we first take asymmetric initial datum as in Fig. 5, and denote the non-zero region are Ω 0 . We take with a small perturbation ∼ 0.1Uniform(Ω 0 ) on the initial datum in the region Ω 0 . We take D = 0.01 and the number of grid points in each spatial direction as 500, 800 and 1000. It can been see in Fig 10, the differences resulting from the varying of the spatial mesh size are of small magnitude as can be seen in both the log scale and the regular scale. the second row 800 by 800; the third row: 500 by 500. Though the pattern seem a bit different in the log scale, it can be seen from the plots that the differences between the three results are in fact minor in the normal scale. Now, we take symmetric initial data and symmetric source S. For γ = 0.75, D = 0.01, the solutions are computed until T = 20 using the number of grid points in each spatial direction as 500, 800 and 1000, respectively. It can been see in Fig 12, the pattern seems symmetric, and varying different mesh sizes does not seem to break the symmetry. It can be seen in the middle column of Fig. 12 that the differences resulting from the change of the grid size are minor, which are more clearly illustrated in the log-scale plots (first and last columns) that the difference is of order 1e-4 marked as dark yellow in the color bar. Example 4.6 (Asymmetric Grids). An interesting phenomenon been pointed out in [START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF] that for asymmetric grids, the solution could form an asymmetric pattern even with symmetric initial data and source. In this last example, we change to an asymmetric grids, namely, dx = dy. We consider the symmetric initial data as in Example 4.5 with γ = 0.75 and D = 0.01, but use 500 grid points in x direction while 800 grid points in y direction. As can be see in Fig. 12, the pattern becomes asymmetric in the log-scale plot, but in fact in the regular scale the solution does not seem to differ much from the results of symmetric grids as plotted in Fig. 11 middle column. We point out that just as the checkerboard problem in computational fluid dynamics, to avoid some potential instability caused by the "natural grids", the choice of grids should certainly be carefully designed. Here we used the staggered grids, and a better gridding strategy for both finite difference and finite element schemes of this model is a very interesting problem, which is left for future study.

Conclusion

In this work, we propose two possible approaches to the Cai-Hu model for biological transport networks -the fully implicit treatment and a semi-implicit approach. While the fully implicit one is proven to be unconditionally stable in 1D, it results in a nonlinear algebraic system that could be difficult and costly to implement in higher dimension, and convergence of the Newton solver remains unclear. The semi-implicit scheme avoids the Newton iteration, but has a stability constraint on time step which depends on some coefficients of the equation (still independent of the mesh size though). For the 2D model, the semi-implicit scheme is of more practical potential, since it only deals with linear algebraic problems by inverting positive definite matrices that can be solved efficiently by CG or PCG even in higher dimensions. Our numerical examples show the stabilities and energy decay of both schemes.

Figure 1 :Figure 2 :

 12 Figure 1: Example 4.1 computed by the implicit solver: the solution m at different times t = 0, 2, 5, 10. It can be seen that the solution converges to 0 eventually, and the discrete physical energy decreases in time.

Figure 3 :

 3 Figure 3: Example 4.2 with D = 0.01. Left panel: computed by the implicit solver; right panel: computed by semi-implicit solver. The equilibrium does not vanish and the discrete physical energy is observed to decay in time.

Example 4 . 3 .

 43 Consider the equation with the Neumann boundary condition for p and Dirichlet boundary condition for m. The background permeability r(x) is chosen as r(x) = 0.1, the activation constant c = 50, the diffusivity D = 0.01, and α = 1.

  1D, T = 20 (b) Semi-implicit scheme: m at t = 20 Semi-implicit scheme: zoom in of discrete energy after t = 10

Figure 4 :

 4 Figure 4: Example 4.2 with D = 0.001. Left panel: computed by the implicit solver; right panel: computed by semi-implicit solver. The equilibrium does not vanish and the discrete physical energy is observed to decay in time for both solvers.

Figure 5 :(

 5 Figure 5: Plot of initial data in the (x,y) plane and the source term.

Figure 6 :

 6 Figure 6: Example 4.3 computed by the semi-implicit solver: the leaf structure |m| 2 in the log scale at different times t = 20, 40, 60, 80. It can be seen that the flow spreads and finally achieves a tree-like structure.

Figure 7 :

 7 Figure 7: Example 4.3 computed by the semi-implicit solver: the flux u defined in(4.3) at different times t = 20, 40, 60, 80. It can be seen that the leaf grows and finally achieves a tree-like structure.

  Zoom in to the tail for 60 ≤ t ≤ 80

Figure 8 :

 8 Figure 8: Example 4.3 computed by the semi-implicit solver: The discrete physical energy decays all the time, and it can be seen in the zoom-in figure (on the right) that it keeps the decay behavior even if the first two decimal places no longer change.

  Discrete Physical energy versus time (d) leaf structure |m| 2 in log-scale (f) Discrete Physical energy versus time (g) leaf structure |m| 2 in log-scale Discrete Physical energy versus time (j) leaf structure |m| 2 in log-scale Discrete Physical energy versus time (m) leaf structure |m| 2 in log-scale (n) flux u in log-scale (o) Discrete Physical energy versus time

Figure 9 :

 9 Figure 9: Example 4.4 computed by the semi-implicit solver for various D till T = 40. From first row to the last row, the values of D are 0.5, 0.1, 0.02, 0.005 and 0.002, respectively. It can be seen that the smaller the diffusion coefficient is, the more detailed the network structure becomes.

|m| 2 (

 2 (a) leaf structure |m| 2 in log-scale 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2|m| 2 (c) leaf structure |m| 2 in log-scale 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2e) leaf structure |m| 2 in log-scale 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2

Figure 10 :

 10 Figure 10: Example 4.5 (asymmetric initial datum) computed by the semi-implicit solver for D = 0.01 till T = 40 for various grid points in each spatial direction. The first row: 1000 by 1000;the second row 800 by 800; the third row: 500 by 500. Though the pattern seem a bit different in the log scale, it can be seen from the plots that the differences between the three results are in fact minor in the normal scale.

  |m| 2 (c) flux u in log-scale (d) leaf structure |m| 2 in log-scale flux u in log-scale (g) leaf structure |m| 2 in log-scale flux u in log-scale

Figure 11 :

 11 Figure 11: Example 4.5 (symmetric initial datum) computed by the semi-implicit solver for D = 0.01 till T = 20 for various grid points in each spatial direction. The first row: 1000 by 1000; the second row 800 by 800; the third row: 500 by 500. The symmetry of the patterns is not broken by varying of the mesh size.

  (a) |m| 2 in log-scale

Figure 12 :

 12 Figure 12: Example 4.6 symmetric initial data with asymmetric grids dx = dy. The pattern becomes asymmetric as well. But as can be seen in the regular scale (on the right), the solution does not differ much from the results of symmetric grids (see Fig. 11 middle column).

The metabolism term |m| 2(γ-1) m is nonlinear and non-differentiable in general for physical most relevant cases[START_REF] Albi | Biological transportation networks: Modeling and simulation[END_REF] 2 ≤ γ ≤ 1, and is also stiff since m can be close to zero. This term is in fact much harder to treat than the diffusion.

The activation matrix (in 2D) c 2 ∇p ⊗ ∇p has eigenvalues 0 and c 2 |∇p| 2 , which makes the activation term stiff when |∇p| is large. This is precisely the case 'inside' the network.In the literature, a number of efforts have been devoted to construct effective numerical schemes for this system.[START_REF] Haskovec | Notes on a pde system for biological network formation[END_REF] and[START_REF] Albi | Biological transportation networks: Modeling and simulation[END_REF] develop numerical schemes with only the diffusion term implicit formulated in either a forward Euler or a Crank Nicolson fashion and other terms explicit, and[START_REF] Albi | Continuum Modeling of Biological Network Formation[END_REF] presents a time-splitting strategy, where the stability of the former remain unclear and the latter has severe time-step constraints due to the stiffness. The goal of this paper is to develop efficient implicit and semi-implicit schemes with good stability properties. First, to effectively take care of the stiffness, we propose a fully implicit scheme for both 1D and 2D Cai-Hu models. For 1D, we prove that the energy dissipation equation is preserved which gives the unconditional stability of the implicit solver. However, as any implicit solver to nonlinear equations, the scheme gives rise to a nonlinear algebraic problem and hence requires the Newton iterations, which makes the implementation for multi-dimensional cases utterly costly and complicated. Even worse, the convergence of the Newton solver becomes unclear for the 2D case. To achieve a practical algorithm with only linear algebraic solvers needed, semi-implicit treatment is desired. However, as is well known, unlike the fully implicit treatment, it is usually difficult to construct a semi-implicit solver that still decays the energy functional (physical energy) of the gradient flow. As the second and main part, we propose a semi-implicit scheme for the Cai-Hu model in both 1D and 2D that does not result in nonlinear algebraic problems. Moreover, we prove for 1D that this scheme indeed decays the physical energy as long as a condition independent of the mesh size is satisfied. With this new solver, one only needs to deal with linear algebraic problems that can be easily solved either by conjugate gradient (CG) or preconditioned conjugate gradient (PCG). This makes the computation in 2D efficient and feasible. Though we only prove the decay of

* Research supported by NSF grants no. DMS-1522184, DMS-1107291: RNMS KI-Net, and NSFC grant No. 31571071.


André F Caldeira
email: caldeiraaf@gmail.com

Christophe Prieur
email: christophe.prieur@gipsa-lab.fr

Valter Daniel Coutinho
email: coutinho@das.ufsc.br

J S Leite

Regional stability and stabilization of a class of linear hyperbolic systems with nonlinear quadratic dynamic boundary conditions

Keywords: dynamic boundary conditions, Poiseuille flow, linear hyperbolic systems, robust control, LMIs

  L'archive ouverte pluridisciplinaire

Introduction

In industrial processes, a large variety of physical systems is governed by hyperbolic partial differential equations (PDEs) such as hydraulic networks [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2× 2 hyperbolic systems over a bounded interval[END_REF][START_REF] Diagne | Lyapunov exponential stability of linear hyperbolic systems of balance laws[END_REF] and gas devices [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF][START_REF] Castillo | Fresh air fraction control in engines using dynamic boundary stabilization of LPV hyperbolic systems[END_REF]. In particular, fluid transport is often modeled by balance laws (or conservation laws when additive and dissipative terms are disregarded) which are hyperbolic PDEs normally used to express the fundamental dynamics of open conservative systems. However, infinite dimensional systems introduce variable time-delays making the closed-loop control much more challenging and, moreover, distributed measurements and actuators are not usually available. As a consequence, it is more common that actuators and measurements are located at the boundaries which is desired in practical applications as, for instance, in the references previously cited. In addition, a large number of numerical techniques based on finite-dimensional tools, which are often used to the stability analysis of PDE systems, provides only approximate solutions.

On the other hand, Lyapunov theory has been largely applied for several decades to deal with the stability analysis and control design of finite dimensional systems described by ordinary differential equations (ODEs) [START_REF] Vidyasagar | Nonlinear Systems Analysis, 2nd Edition[END_REF]. In the particular case of linear dynamical systems, a large number of stability and stabilization results are cast in terms of linear matrix inequality (LMI) constraints [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], which are numerically solved using dedicated software [START_REF] Lofberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF]. The LMI framework is a powerful tool for linear and nonlinear finitedimensional systems, since it can deal with a large diversity of control and systems theory problems such as robust stability, domain of attraction estimation, input-to-output performance, state or dynamic output-feedback control, and state estimation (see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF][START_REF] Boyd | Convex optimization[END_REF][START_REF] Papachristodoulou | A tutorial on sum of squares techniques for systems analysis[END_REF][START_REF] Coutinho | L 2 -gain analysis and control of uncertain nonlinear systems with bounded disturbance inputs[END_REF][START_REF] Chesi | LMI techniques for optimization over polynomials in control: A survey[END_REF][START_REF] Duan | LMIs in Control Systems[END_REF] among other references). As a result, LMIs have been successfully applied in a wide diversity of control oriented applications as, for instance, wind turbine operation, satellite attitude regulation, turbo-charged combustion engine control and bioprocess control and estimation [START_REF] Mohammadpour | Control of Linear Parameter Varying Systems with Applications[END_REF][START_REF] Coutinho | A robust non-linear feedback control strategy for a class of bioprocesses[END_REF].

In the context of infinite dimensional systems, the stability problem of boundary control for first-order hyperbolic systems using quadratic strict Lyapunov functions was more recently stated; see, for instance, [START_REF] Coron | A Strict Lyapunov Function for Boundary Control of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF][START_REF] Prieur | Stability of switched linear hyperbolic systems by lyapunov techniques[END_REF]. A common assumption in most of available results is that the boundary action is faster than the wave travel making possible to establish a static relationship between the control input and the boundary condition. However, there are applications where the dynamics associated to the boundary control action cannot be neglected as, for instance, in the temperature control of an airflow in a heating column. To deal with dynamic boundary actions, the infinitedimensional system discretization together with the use of finite-dimensional tools are often employed. On the other hand, several approaches considering infinite-dimensional based techniques have been recently proposed. To cite a few, some sufficient conditions for the exponential stability of hyperbolic systems with linear time invariant dynamic boundary conditions were proposed in [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF], the backstepping technique for the boundary control of hyperbolic and time-delayed systems was considered in [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic pdes and application to systems with actuator and sensor delays[END_REF], the dynamic boundary stabilization of linear parameter-varying (LPV) hyperbolic systems using the LMI framework was studied in [START_REF] Castillo | Fresh air fraction control in engines using dynamic boundary stabilization of LPV hyperbolic systems[END_REF], the problem of stability analysis and control synthesis for first-order hyperbolic linear PDEs over a bounded interval with spatially varying coefficients was the focus of [START_REF] Lamare | An optimisation approach for stability analysis and controller synthesis of linear hyperbolic systems[END_REF], and a multimodel approach with a bilinear matrix inequality to deal with the stability of a hyperbolic system representing the flow in a open channel was proposed in [START_REF] Santos Martins | A multimodel approach to Saint-Venant equations: A stability study by LMIs[END_REF]. It turns out that the majority of the latter references consider only dynamic boundary actions described by linear differential equations.

This paper addresses the boundary stabilization of uncertain first-order hyperbolic systems with boundary actions governed by an uncertain nonlinear quadratic differential equation. It should be emphasized that the class of quadratic systems can describe state-space models containing quadratic nonlinearities in the state variables and bilinear terms involving the state and control signal. Moreover, this class of systems can represent a large number of processes such as distillation and heating columns [START_REF] Espana | Reduced order bilinear models for distillation columns[END_REF][START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF], induction motors [START_REF] Arguello-Serrano | Nonlinear control of a Heating, Ventilating, and Air conditioning system with thermal load estimation[END_REF] and DC-DC converters [START_REF] Thounthong | A new control law based on the differential flatness principle for multiphase interleaved DC-DC converter[END_REF], which has attracted recurring interest of the control practitioners; see, for instance, [START_REF] Amato | On the region of attraction of nonlinear quadratic systems[END_REF][START_REF] Amato | Stabilization of bilinear systems via linear state-feedback control[END_REF][START_REF] Valmórbida | State feedback design for input-saturating quadratic systems[END_REF][START_REF] Coutinho | Nonlinear state feedback design with a guaranteed stability domain for locally stabilizable unstable quadratic systems[END_REF] and references therein. Notice, in the context of infinite dimensional systems, that boundary actions are often nonlinear as, for instance, in the geometric average particle volume regulation of plug-low reactors, flow regulation in open channels delimited by overflow spillways and the vibration rejection in flexible marine risers [START_REF] Christofides | Nonlinear and Robust Control of PDE Systems -Methods and Applications to Transport-Reaction Processes[END_REF][START_REF] Coron | A Strict Lyapunov Function for Boundary Control of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Ge | Boundary Control of a Coupled Nonlinear Flexible Marine Riser[END_REF]. In particular, local boundary stabilization conditions for the coupled linear first-order hyperbolic system and nonlinear quadratic dynamic boundary actuation are in this paper proposed in terms of a finite set of LMI constraints based on the Lyapunov stability theory. The local stability is characterized in a regional context, that is, bounded initial states will imply that the state trajectories remain bounded and they converge to the equilibrium point as the time goes to infinity. Hence, an optimization problem is proposed for the boundary feedback control design in order to de-termine a quadratic boundary state feedback control law that maximizes the set of admissible boundary initial conditions while guaranteeing the stability of the coupled ODE-PDE system.

The remaining of this paper is organized as follows. The problem of interest is stated in Section 2. Then, Section 3 presents the main stability results for uncertain first-order linear hyperbolic system subject to nonlinear boundary conditions. Next, the stability results are adapted in Section 4 to control synthesis by means of an appropriate similarity transformation where a stabilizing nonlinear boundary state feedback control law is designed in order to maximize the region of admissible initial conditions. Then, Section 5 provides two application-based examples in order to demonstrated the potentials of the proposed approach. Finally, Section 6 ends the paper giving some concluding remarks.

Notation: R + = [0, +∞), R n is the n-dimensional Euclidean space, R m×n is the set of m × n real matrices with real entries, I n in the n × n identity matrix, 0 n is the n × n zeros matrix and diag{• • •} denotes a block-diagonal matrix. For a real matrix S, S denotes its transpose, He{S} stands for S + S , and S > 0 means that S is symmetric and positive-definite. For two polytopes X and Y, X × Y denotes the meta-polytope obtained by cartesian product and V(X ) is the set of all vertices of X . For a vector ξ ∈ R n , the k-component of ξ is denoted by ξ (k) . The usual Euclidian norm in R n is denoted by • , whereas the set of all functions φ : (0, L) → R n such that L 0 φ(x) φ(x) dx < ∞ is denoted by L 2 ((0, L); R n ) that is equipped with the norm • L 2 ((0,L);R n ) . Given a topological set J , and an interval I in R + , the set C 0 (I, J ) is the set of continuous functions φ : I → J .

Problem Statement

Let n be a positive integer, Ω an open non-empty set of R n and ∆ a non empty convex set of R n δ . Consider the general class of first-order uncertain hyperbolic systems of order n defined as follows:

∂ t ξ(t, x) + Λ(δ)∂ x ξ(t, x) = 0, t ∈ R + , x ∈ [0, L], (1)
where ξ : R + × [0, L] → Ω, δ ∈ ∆ is a bounded sufficiently smooth vector function of time; ∂ t and ∂ x represent (first order) partial derivatives with respect to time and space, respectively; Λ(δ) : ∆ → R n×n is a diagonal and continuous matrix function (called the characteristic matrix), i.e., Λ(δ) = diag{λ 1 (δ), λ 2 (δ), . . . , λ n (δ)} with λ i (δ) being positive continuous functions, for all δ ∈ ∆ and i = 1, . . . , n.

Assumption 1. The diagonal elements of Λ(δ) ∈ R n×n satisfy the following inequalities for all δ ∈ ∆:

0 < λ 1 (δ) < • • • < λ n (δ).
Remark 1. Assumption 1 implies that system (1) has only positive convecting speeds, however it does not add any loss of generality. Note that, for first order hyperbolic systems with both negative and positive convecting speeds, there always exists a variable transformation to obtain a system as in (1) satisfying Assumption 1 (see the beginning of proof of Theorem 1 in [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF], or Page 7 of [START_REF] Tang | Tikhonov theorem for linear hyperbolic systems[END_REF] or Remark 1 of [START_REF] Castillo | Fresh air fraction control in engines using dynamic boundary stabilization of LPV hyperbolic systems[END_REF] or Remark 3 of [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF]). In particular, consider system (1) with the diagonal elements of Λ(δ) ordered as follows

λ 1 (δ) < • • • < λ m (δ) < 0 < λ m+1 (δ) < • • • < λ n (δ), ∀ δ ∈ ∆,
and let ξ : R + × [0, L] → Ω ∈ R n be partitioned as

ξ(t, x) = ξ -(t, x) ξ + (t, x) , ξ -∈ R m , ξ + ∈ R n-m .
Then, applying the transformation

ξ(t, x) = ξ -(t, L -x) ξ + (t, x)
yields a system as in (1) with only positive convection speeds.

Associated with (1), consider the following nonlinear quadratic dynamic boundary action:

ξin (t) = A(ξ in (t), δ)ξ in (t) + B(ξ in (t), δ)u(t), u(t) = G(ξ in (t))ξ in (t) + Kξ out (t), (2)
where

ξ in (t) = ξ(t, 0) , ξ out (t) = ξ(t, L) (3)
are the boundary conditions of (1) interconnecting the boundary dynamics of (1) with (2), ξ in ∈ Ξ ⊂ R n and ξ out ∈ R n are the state vectors of the dynamic boundary condition at x = 0 and x = L, respectively, and assumed to be measurable, u(t) ∈ R nu is the boundary control input, Ξ is compact region of the boundary state-space containing ξ in = 0 and to be specified later in this paper. Notice to ease the presentation that the system and boundary dynamics uncertainties are merged into the vector δ, i.e., the elements of δ comprise both system and boundary dynamics modeling errors. Further, K ∈ R nu×n is a real constant matrix and

A(•) ∈ R n×n , B(•) ∈ R n×nu , G(•) ∈ R nu×n
are affine matrix functions on their arguments, that is:

[A(ξ in , δ) B(ξ in , δ)] = [A 0 B 0 ] + n i=1 ξ in (i) [A i B i ] + n δ j=1 δ (j) Ȃj Bj , (4) G(ξ in ) = G 0 + n i=1 ξ in (i) G i , (5)
with A i , B i , and G i , i = 0, 1, . . . , n, and Ȃj , Bj , j = 1, . . . , n δ , being given constant real matrices with appropriate dimensions. It is assumed that the pair (A 0 , B 0 ) is stabilizable and the unforced system of ( 2) is allowed to be unstable.

In addition, the initial conditions of the coupled PDE-ODE system of (1) and ( 2), namely

ξ(0, x) = ξ(x) and ξ in (0) = ξ in , (6)
are assumed to satisfy:

ξ ∈ D 1 and ξ in ∈ D 2 , (7)
where the sets D 1 and D 2 are defined as follows

D 1 := ξ ∈ L 2 ((0, L); R n ) : L 0 ξ(x) ξ(x) dx ≤ σ , (8)
D 2 := ξ in ∈ R n : ξ in P 1 ξ in ≤ 1 ,
with the scalar σ > 0 and the matrix P 1 > 0 defining the sizes of D 1 and D 2 , respectively. In the particular case of linear dynamic boundary conditions, that is, when A i , B i , Ǎi and Bi are null matrices, for i = 1, . . . , n, the existence and uniqueness of solutions to (2), ( 3) and ( 6) for initial conditions satisfying ( 7) is ensured, see e.g. [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF]. As these solutions may not be differentiable everywhere, the concept of weak solutions to partial differential equations has to be used; see [START_REF] Diagne | Lyapunov exponential stability of linear hyperbolic systems of balance laws[END_REF] and references therein for further details.

In view of the above scenario, this paper is concerned in obtaining sufficient conditions to the regional stability and stabilization problems for coupled PDE-ODE system of ( 1) and ( 2) as stated below: P 1 For given state feedback gain matrices G(ξ in ) and K(ξ out ) as defined in ( 5), derive stability analysis conditions ensuring the robust local stability of the coupled PDE-ODE system while determining a maximized set of initial boundary states D 2 for a given set D 1 . P 2 Design the state feedback gain matrices G(ξ in ) and K(ξ out ) as defined in ( 5), ensuring the robust local stability of the coupled PDE-ODE system while maximizing the set of initial boundary states D 2 for a given set D 1 .

Before ending this section, the notion of exponential stability to be considered in this paper is introduced for the coupled PDE-ODE system of ( 1) and ( 2). Definition 1. The coupled PDE-ODE system of (1) and (2), with initial conditions ξ and ξ in satisfying [START_REF] Boyd | Convex optimization[END_REF], is said to be locally robustly exponentially stable if there exist positive scalars α and β such that the following holds:

ξ in (t) + ξ(t, •) L 2 ((0,L);R n ) ≤ βe -αt ξ in + ξ L 2 ((0,L);R n ) , ∀ t ∈ R + , δ ∈ ∆. (9)

Local Stability Analysis

In this section, it is developed an LMI approach to derive a numerical and tractable solution to the robust regional stability analysis problem P 1 as defined in Section 2. To this end, consider the following Lyapunov function candidate defined for all continuously differentiable functions ξ : R + ×[0, L] → Ω:

V (ξ) = ξ in P 1 ξ in + L 0 e -µx ξ(x) P 2 ξ(x) dx, (10)
where P 1 , P 2 ∈ R n×n are positive definite diagonal matrices and µ is a positive scalar.

Then, evaluating the time derivative of V (•) along the solutions to (1) and ( 2) with (7) leads to (see, for instance, [START_REF] Coron | A Strict Lyapunov Function for Boundary Control of Hyperbolic Systems of Conservation Laws[END_REF] and [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF] for further details):

V (ξ(t, •)) = 2ξ in (t) P 1 A(ξ in , δ) + B(ξ in , δ)G(ξ in ) ξ in (t) + 2ξ in (t) P 1 B(ξ in , δ)Kξ out (t) -e -µx ξ(t, x) Λ(δ)P 2 ξ(t, x) x=L x=0 -µ L 0 e -µx ξ(t, x) Λ(δ)P 2 ξ(t, x) dx
In light of the above and taking (3) into account, V (ξ(t, •)) can be cast as follows:

V (ξ(t, •)) = -m 1 (ξ, δ) + m 2 (ξ in , ξ out , δ) (11)
where

m 1 (ξ, δ) = µ ξ in Λ(δ)P 1 ξ in + L 0 e -µx ξ(t, x) Λ(δ)P 2 ξ(t, x) dx , m 2 (ξ in , ξ out , δ) = ξ a Φ(ξ in , δ)ξ a , ξ a = ξ in ξ out , Φ(ξ in , δ) =       A(ξ in , δ) P 1 +P 1 A(ξ in , δ) +G(ξ in ) B(ξ in , δ) P 1 +P 1 B(ξ in , δ)G(ξ in ) +Λ(δ)P 2 + µΛ(δ)P 1 P 1 B(ξ in , δ)K K B(ξ in , δ) P 1 -e -µ Λ(δ)P 2       .
Due to Assumption 1 and since Λ(•) is a continuous function on the compact set ∆, there exists a sufficiently small positive scalar such that ≤ λ 1 (δ), for all δ ∈ ∆, which implies that I n ≤ Λ(δ) for all admissible δ.

Then, provided that

m 2 (ξ in , ξ out , δ) ≤ 0 , ∀ (ξ in , δ) ∈ Ξ × ∆, ξ out ∈ R n (12)
for a set Ξ (that will be defined latter), the following holds

V (ξ(t, •)) ≤ -µ V (ξ(t)), ∀ t ∈ R + , (13)
from the fact that m 1 (ξ(t, •), δ) ≥ µ V (ξ(t, •)).

Notice that the inequality [START_REF] Chesi | LMI techniques for optimization over polynomials in control: A survey[END_REF] implies that the coupled PDE-ODE system is locally exponentially stable if the condition ( 12) is satisfied along the solutions to (1) and ( 2) for all ξ ∈ D 1 and ξin ∈ D 2 . Hence, we have to guarantee that the solution to the coupled PDE-ODE system is confined to a region R(γ

) := {ξ : V (ξ) ≤ γ} (14)
for a certain γ > 0.

Firstly, from ( 13), the following holds:

V (ξ(t, •)) ≤ e -µ t V ( ξ) , ∀ δ ∈ ∆.
Next, provided that ξ ∈ D 1 , ξ in ∈ D 2 and µ > 0, there exists a sufficiently large positive scalar such that the following holds:

L 0 e -µx ξ(x) P 2 ξ(x) dx ≤ L 0 ξ(x) P 2 ξ(x) dx ≤ L 0 ξ(x) ξ(x) dx ≤ σ.
In the above inequality, observe that is any scalar equal or larger than the largest eigenvalue of P 2 .

Then, the following upper bound on V ( ξ) is obtained

V ( ξ) ≤ 1 + σ,
as soon as ξ ∈ D 1 , ξ in ∈ D 2 which, in view of the above, implies that:

V (ξ(t, •)) ≤ γ , γ = 1 + σ, ∀ t ∈ R + , δ ∈ ∆. (15)
Now, let V 1 (ξ in ) = ξ in P 1 ξ in , then it follows from (10) that:

V 1 (ξ in (t)) ≤ V (ξ(t, •)) ≤ γ, ∀ t ∈ R + , δ ∈ ∆. ( 16
)
Let the following set:

R 1 := {ξ in : V 1 (ξ in ) = ξ in P 1 ξ in ≤ γ}. (17)
Thus, the condition in [START_REF] Coron | A Strict Lyapunov Function for Boundary Control of Hyperbolic Systems of Conservation Laws[END_REF] guarantees that

ξ in ∈ R 1 for all t ∈ R + and δ ∈ ∆, since V 1 (ξ in (t)) ≤ V (ξ(t, •)) ≤ V ( ξ) ≤ γ for all t ∈ R + and δ ∈ ∆.
Then, provided that R 1 ⊂ Ξ, the condition in [START_REF] Castillo | Fresh air fraction control in engines using dynamic boundary stabilization of LPV hyperbolic systems[END_REF] holds along the solutions to (1) and ( 2) for all ξ ∈ D 1 and ξ in ∈ D 2 .

In the sequel, it is introduced the main result of this section which proposes a numerical and tractable solution to problem P 1 in terms of a finite set of LMI constraints. To this end, it is assumed that Ξ is a given symmetric hyper-rectangle (i.e., a polytopic region) containing ξ in = 0 with known vertices. Moreover, Ξ can be also equivalently defined in terms of its faces as below

Ξ = {ξ in ∈ R n : |c i ξ in | ≤ 1, i = 1, . . . , n f }, (18)
with c i ∈ R n , i = 1, . . . , n f , defining the faces of Ξ.

Theorem 1 (Stability Analysis). Consider the PDE-ODE system (1) and ( 2), with the initial conditions defined by ( 7) and ( 8), and given control gains G(ξ in ) and K. Let Ξ and ∆ be given polytopes, and σ, µ be given positive scalars. Suppose there exist diagonal matrices P 1 and P 2 , a matrix W , with appropriate dimensions, and a positive scalar satisfying the following:

P 1 > 0, P 2 > 0, I n -P 2 ≥ 0, ( 19
)
γ γc j γc j P 1 ≥ 0, j = 1, . . . , n f , (20)
Θ(ξ in , δ, µ) + He{W M(ξ in )} < 0, ∀ (ξ in , δ) ∈ V(Ξ × ∆), (21)
where γ = 1 + σ and

Θ(ξ in , δ, µ) =    N A(ξ in , δ) P 1 +P 1 A(ξ in , δ) + Λ(δ)P 2 + µΛ(δ)P 1 N +G B(ξ in , δ) P 1 N +N P 1 B(ξ in , δ)G N P 1 B(ξ in , δ)K K B(ξ in , δ) P 1 N N (-e -µ Λ(δ)P 2 )N    , N = I n 0 n×n 2 , G = G 0 G 1 • • • G n , M(ξ in ) = Π(ξ in ) -I n 2 0 0 N (ξ in ) 0 , Π(ξ) =    ξ (1) ⊗ I n . . . ξ (n) ⊗ I n    , N (ξ) =         ξ (2) I n -ξ (1) I n 0 • • • • • • 0 0 ξ (3) I n -ξ (2) I n 0 • • • 0 . . . . . . . . . . . . • • • . . . . . . • • • . . . . . . . . . . . . 0 • • • • • • 0 ξ (n) I n -ξ (n-1) I n         . ( 22
)
Then, the origin of the coupled PDE-ODE system of ( 1) and ( 2) is locally robustly exponentially stable in the sense of Definition 1. Moreover, for any ξ ∈ D 1 and ξ in ∈ D 2 , the system trajectories remain bounded to R 1 as defined in [START_REF] Coutinho | Nonlinear state feedback design with a guaranteed stability domain for locally stabilizable unstable quadratic systems[END_REF], for all t ≥ 0, and vanish to zero as the time goes to infinity.

Proof 1. Suppose the conditions of Theorem 1 are satisfied, then from convexity arguments the condition in ( 21) is also satisfied for all

(ξ in , δ) ∈ Ξ × ∆, ξ out ∈ R n .
Firstly, let us show that (21) implies that (12) holds. To do that, let

Ψ(ξ in ) = Π a (ξ in ) 0 0 I n , Π a (ξ) = I n Π(ξ)
.

Hence, pre-and post-multiplying (21) by Ψ(ξ in ) and Ψ(ξ in ), respectively, yields

Φ(ξ in , δ) < 0, ∀ (ξ in , δ) ∈ Ξ × ∆,
since M(ξ in )Ψ(ξ in ) = 0 by construction, and thus

m 2 (ξ in , ξ out , δ) ≤ 0, ∀ (ξ in , δ) ∈ Ξ × ∆, ξ out ∈ R n .
Next, in light of the definition of V (ξ) in [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF] and taking (11)-( 12) into account, the following holds:

V (ξ(t, •)) ≤ -µ V (ξ(t)), ∀ t ∈ R + .
Now, notice that (20) implies R 1 ⊂ Ξ; see, e.g., [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. Hence, for all ξ ∈ D 1 and

ξ in ∈ D 2 , it follows that V 1 (ξ in (t)) ≤ V (ξ(t, •)) ≤ V ( ξ) ≤ γ for t ∈ R + and δ ∈ ∆, which completes the proof.
Theorem 1 can be applied to obtain the largest estimate of the set D 2 of admissible initial boundary states assuming that D 1 and Ξ are given a priori. For instance, the volume of D 2 can be approximately maximized by minimizing the trace of P 1 (since the trace of P -1 1 is the sum of the squared semi-axis lengths of D 2 ) leading to the following optimization problem min µ, ,P 1 ,P 2 ,W trace{P 1 } subject to ( 19)- [START_REF] Diagne | Lyapunov exponential stability of linear hyperbolic systems of balance laws[END_REF]. [START_REF] Santos Martins | A multimodel approach to Saint-Venant equations: A stability study by LMIs[END_REF] Notice that the matrix inequalities in ( 19)-( 21) become LMIs when µ is given. Hence, it can be applied a line search over µ in order to obtain a solution to optimization problem in [START_REF] Santos Martins | A multimodel approach to Saint-Venant equations: A stability study by LMIs[END_REF] via semi-definite programming.

Remark 2. To be precise, the Lyapunov function in [START_REF] Castillo | Dynamic boundary stabilization of linear and quasi-linear hyperbolic systems[END_REF] is not differentiable along all L 2 solutions. However it is continuously differentiable along the smooth solutions, and thus, by density, the analysis developed in this section is also correct for L 2 solutions . See, for instance, [5, page 67] for further details.

Regional Stabilization

This section addresses the problem of designing the state feedback gain matrices G(•) and K(•) such that the coupled PDE-ODE system of ( 1) and ( 2) with ( 7) is regionally exponentially stable for all δ ∈ ∆. The proposed design will be based on Theorem 1 together with a similarity transformation and matrix parameterizations involving the Lyapunov matrix P 1 .

Thus, consider the following definitions:

η in = P 1 ξ in , η out = P 1 ξ out , η = P 1 ξ, Q 1 = P 1 -1 , Q 2 = Q 1 P 2 Q 1 , G p (ξ in ) = G(ξ in )Q 1 , K p = KQ 1 , (24)
with Q 1 > 0. Now, observe that m 1 (ξ, δ) and m 2 (ξ in , ξ out , δ) in ( 11) can be respectively recast as follows:

m 1 (ξ, δ) = µ η in Q 1 Λ(δ)η in + L 0 e -µx η(t, x) Q 2 Λ(δ)η(t, x) dx , (25) m 2 (ξ in , ξ out , δ) = η a Φ p (ξ in , δ)η a , (26)
where

η a = η in η out , Φ p (ξ in , δ) =     Q 1 A(ξ in , δ) +A(ξ in , δ)Q 1 + G p (ξ in ) B(ξ in ) + B(ξ in )G p (ξ in ) +Q 2 Λ(δ) + µQ 1 Λ(δ) B(ξ in , δ)K p K p B(ξ in , δ) -e -µ Λ(δ)Q 2     .
Hence, using similar arguments to Section 3, the condition V (ξ(t, •)) ≤ -µ V (ξ(t, •)) will hold for t ∈ R + if m 2 (ξ in , ξ out , δ) as defined in (26) satisfies

m 2 (ξ in , ξ out , δ) ≤ 0, ∀ (ξ in , δ) ∈ Ξ × ∆, ξ out ∈ R n . (27)
In addition, notice that the condition I n -P 2 ≥ 0 was considered in Section 3 to obtain a bound γ on V (ξ(t, •)) such that ( 16) holds for all t ∈ R + and δ ∈ ∆. To this end, the following result will be instrumental to obtain a numerical and tractable solution.

Lemma 1 ([20]

). Let P and R be real square matrices with P > 0 and R nonsingular. Then, the following holds:

R P -1 R ≥ R + R -P (28)
Because of the parametrization of P 2 in (24), it follows that I n -P 2 ≥ 0 can be equivalently written as follows:

Q 1 Q -1 2 Q 1 -ϕI n ≥ 0, (29)
with ϕ = -1 , and noting that

P 2 = Q -1 1 Q 2 Q -1 1 .
Hence, the following inequality is a sufficient condition for (29) to hold:

2Q 1 -Q 2 -ϕI n ≥ 0, (30)
by applying Lemma 1 with R = Q 1 and P = Q 2 .

In view of the above developments, it is proposed in the sequel an LMIbased result for designing the state feedback gain matrices G(ξ in ) and K such that the origin of the coupled PDE-ODE of (1) and ( 2) is robustly regionally stable in closed-loop.

Theorem 2 (Regional Stabilization). Consider the PDE-ODE system (1) and (2), with the initial conditions defined by [START_REF] Boyd | Convex optimization[END_REF] and [START_REF] Caldeira | Modeling and control of flow with dynamical boundary actions[END_REF]. Let Ξ and ∆ be given polytopes, and ϕ, σ and µ be given positive scalars. Suppose there exist diagonal matrices Q 1 and Q 2 , and matrices W , S, F i , i = 0, 1, . . . , n, with appropriate dimensions, satisfying the following:

Q 1 > 0, Q 2 > 0, 2Q 1 -Q 2 -ϕI n ≥ 0, ( 31
) 1 -γc j Q 1 c j ≥ 0, j = 1, . . . , n f , ( 32
) Θ p (ξ in , δ, µ) + He{W M(ξ in )} < 0, ∀ (ξ in , δ) ∈ V(Ξ × ∆), ( 33
)
where γ = 1 + σϕ -1 and

Θ p (ξ in , δ, µ) =    N Q 1 A(ξ in , δ) +A(ξ in , δ)Q 1 +F B(ξ in , δ) N +N B(ξ in , δ)F + N Λ(δ)Q 2 + µΛ(δ)Q 1 N N B(ξ in , δ)S S B(ξ in , δ) N -e -µ Λ(δ)Q 2    , F = F 0 F 1 • • • F n . ( 34
)
Then, the origin of the coupled PDE-ODE system (1) and ( 2), with

G(ξ in ) = F 0 Q -1 1 + n i=1 ξ in (i) F i Q -1 1 , K = SQ -1 1 , (35)
is locally robustly exponentially stable in the sense of Definition 1. Moreover, for any ξ ∈ D 1 and ξ in ∈ D 2 , the solution to the coupled PDE-ODE system remains bounded to R 1 as defined in [START_REF] Coutinho | Nonlinear state feedback design with a guaranteed stability domain for locally stabilizable unstable quadratic systems[END_REF], with

P 1 = Q -1 1 and P 2 = Q -1 1 Q 2 Q -1
1 , for all t ≥ 0, and vanish to zero as the time goes to infinity.

Proof 2. Suppose the conditions of Theorem 2 are satisfied, then from convexity arguments the condition in [START_REF] Prieur | Stability of switched linear hyperbolic systems by lyapunov techniques[END_REF] is also satisfied for all (ξ in , δ)

∈ Ξ × ∆, ξ out ∈ R n .
Firstly, let us show that [START_REF] Prieur | Stability of switched linear hyperbolic systems by lyapunov techniques[END_REF] implies that (27) holds. To do that, let

Ψ(ξ in ) = Π a (ξ in ) 0 0 I n , Π a (ξ) = I n Π(ξ)
.

Hence, pre-and post-multiplying (33) by Ψ(ξ in ) and Ψ(ξ in ), respectively, yields

Φ p (ξ in , δ, µ) < 0, ∀ (ξ in , δ) ∈ Ξ × ∆,
since M(ξ in )Ψ(ξ in ) = 0 by construction. Now, pre-and post multiplying the above by η a and η a , respectively, leads to [START_REF] Castillo | Fresh air fraction control in engines using dynamic boundary stabilization of LPV hyperbolic systems[END_REF] taking the definitions of η in and η out in (24) into account and the fact that:

G(ξ in ) = FΠ a (ξ in )P 1 , K = SP 1 .
Next, notice that (32) implies D ∈ Ξ; see, e.g., [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. Then, consider the condition 2Q 1 -Q 2 -ϕI n ≥ 0 on the right-hand side of [START_REF] Mohammadpour | Control of Linear Parameter Varying Systems with Applications[END_REF]. From Lemma 1, it follows that:

P -1 2 -ϕI n ≥ 0, since P 2 = Q -1 1 Q 2 Q -1 1 .
Hence, multiplying the above by P 2 with = ϕ -1 , the following is obtained

I n -P 2 ≥ 0.
Thus, the rest of this proof follows straightforwardly from Theorem 1.

Similarly to the stability analysis counterpart, Theorem 2 can be applied for designing the state feedback gain matrices G(ξ in ) and K, as defined in [START_REF] Bastin | Stability and Boundary Stabilization of 1-D Hyperbolic Systems[END_REF], in order to robustly regionally stabilize the coupled PDE-ODE system of ( 1) and ( 2) in closed-loop while maximizing the volume of D 2 for a given set D 1 . Thus, the following optimization problem solves problem P 2: 31)- [START_REF] Prieur | Stability of switched linear hyperbolic systems by lyapunov techniques[END_REF]. [START_REF] Valmórbida | State feedback design for input-saturating quadratic systems[END_REF] Notice that the matrix inequalities in ( 31)-( 33) are LMIs when ϕ and µ are given a priori. Hence, a solution to the above optimization problem is obtained via semi-definite programming by applying a gridding technique over ϕ and µ.

max µ,ϕ,Q 1 ,Q 2 ,F 0 ,...,S,W logdet(Q 1 ) subject to (

Illustrative Examples

In this section, two application-based examples are presented to illustrate the results derived in the paper. The first one consists on designing a boundary controller to regulate the outlet temperature in a Poiseuille flow system. Then, in the second one, it is performed a local stability analysis of a lossless transmission line connected to a constant power load.

Poiseuille Flow System

The fluid transport system is normally used for industrial applications such as the ventilation system of mining industry [START_REF] Witrant | Air flow modelling in deep wells: aplication to mining ventilation[END_REF] and hydraulic networks [START_REF] Santos | Boundary control of open channels with numerical and experimental validations[END_REF]. To investigate the phenomenon of fluid transport in a Poiseuille flow subject to dynamic boundary conditions, an experimental setup has been designed to test and validate the proposed results.

The experimental setup consists on a heating column encasing a resistor, a tube, two ventilators, a gas speed meter and three distributed temperature sensors [START_REF] Caldeira | Modeling and control of flow with dynamical boundary actions[END_REF]. The control objective is tracking the outlet temperature by driving the power dissipated on the heating resistor at different air flow speeds through the tube. Only the outlet temperature and the flow speed will be considered as measurements for a closed-loop boundary control strategy and thus G(ξ in ) is set to be zero. In this work, it is assumed that the output ventilator is off.

The modeling of the experimental setup is done by considering two subsystems: the heating column and the tube. One-dimensional transport model is used to describe the gas density, speed and pressure variations in the tube. For the dynamic boundary conditions, it is considered a zero-dimensional model of control volume approach with heat exchanges coming from the heating resistor in the column. The perfect gases law is used to convert density on temperature.

System description and modeling are presented in [START_REF] Castillo | Dynamic boundary stabilization of linear parameter varying hyperbolic systems: Application to a Poiseuille flow[END_REF][START_REF] Caldeira | Modeling and control of flow with dynamical boundary actions[END_REF]. Applying some simplifications (according [9, Section 4]) the system of PDEs reduces to the following convection equation

∂ρ(t, x) ∂t + v(t) ∂ρ(t, x) ∂x = 0, (37)
where v(t) is the time-varying convection parameter of [START_REF] Vidyasagar | Nonlinear Systems Analysis, 2nd Edition[END_REF]. The gas speed v(t) in the tube is measured.

For the dynamic boundary condition, considering that the input mass flow rate to be expressed as ṁin = v(t)ρ(t, 0)A t . Thus, the dynamic boundary condition can be expressed as:

ρ(t, 0) = - RγT in v(t)A t p in V 0 ρ(t, 0) 2 + γv(t)A t V 0 ρ(t, 0) - R p in V 0 C v ρ(t, 0)dQ, ( 38
)
where A t is the tube cross section area.

Temperature Boundary Control

The output temperature boundary control is designed for [START_REF] Vidyasagar | Nonlinear Systems Analysis, 2nd Edition[END_REF] with boundary conditions [START_REF] Witrant | Air flow modelling in deep wells: aplication to mining ventilation[END_REF]. Let define the density error as:

ξ = ρ -ρ ref , (39)
where ρ ref is the desired output density. It is easy to show that system [START_REF] Vidyasagar | Nonlinear Systems Analysis, 2nd Edition[END_REF] with boundary conditions [START_REF] Witrant | Air flow modelling in deep wells: aplication to mining ventilation[END_REF] can be expressed as follows:

∂ξ(t, x) ∂t + Λ(δ) ∂ξ(t, x) ∂x = 0, (40)
with the boundary conditions:

ξin = -A(δ)ξ in 2 + a(δ)ξ in -Bξ in dQ, ξ 0 = ξ in , (41)
where

A(δ) = RγT in v(δ)A t p in V 0 , a(δ) = γv(δ)A t V 0 and B = R p in V 0 C v
The speed v(δ) is measured and is considered an uncertain parameter. The temperature inside the heating column is not measurable. Let define the control input using a feedback (partial) and feedforward linearization as:

dQ = F (δ, ρ ref ) + Kξ out , ξ out = ξ(t, L), (42)
with

F (δ, ρ ref ) = C v γv(δ)A t p in R -T in ρ ref and K to be designed.
This yields to the system (40) with boundary condition (41), which corresponds to the system considered in Theorem 1 and Theorem 2. Let define the convex set ∆ = [START_REF] Amato | Stabilization of bilinear systems via linear state-feedback control[END_REF][START_REF] Arguello-Serrano | Nonlinear control of a Heating, Ventilating, and Air conditioning system with thermal load estimation[END_REF] as the flow speed 1 m/s and 3 m/s in order to operate the experimental setup between the temperatures of 300 K and 325 K.

In order to design the boundary control for system (40) with boundary condition (41), consider the control architecture given in Figure 1. In addition, let the following system parameters: adiabatic constant γ = 1.4; molar mass of dry air M = 28.97 g/mol.K; ideal gas constant R = 8.3143 J/(mol.K); specific heat constant for constant pressure C p = 1.005 KJ/Kg.K for constant volume C v = 0.718 KJ/Kg.K; initial pressure p in = 1 Bar (or, equivalently, 1 × 10 5 Pa in SI units); initial temperature T in = 300 K; initial density ρ = 1.1768 Kg/m 3 column volume V 0 = 4 × 10 -3 m 3 ; tube cross section area A t = 6.4 × 10 -3 m 2 ; and the tube length L = 1 m.

In this scenario, we have applied Theorem 1 with [START_REF] Santos Martins | A multimodel approach to Saint-Venant equations: A stability study by LMIs[END_REF] to show that the closed-loop system illustrated in Figure 1 In a scenario similar to that shown in Figure 3, now considering v(t) = 2 + sin(12.56 t). Note that in this case, the couple PDE-ODE system is time variant. Without lost of generality, and using gain scheduling, the Theorem 2 is applied with [START_REF] Valmórbida | State feedback design for input-saturating quadratic systems[END_REF] and considering K(ξ out (t)) = K 0 + δ(t)ξ out K 1 , µ = 0.5

K 0 = -0.6., K 1 = 7.8. ( 44
)
Figure 5 presents the closed-loop simulations for the latter control gains. Notice that the outlet temperature's trajectory robustly follows the reference with a small degradation of transient behavior.

Transmission line connected to a Constant Power Load (CPL)

Consider the lossless electrical line connecting a power supply to a constant power load as shown in Figure 6. The electrical line dynamics (repre-

∂ t I(t, x) + 1 L ∂ x V (t, x) = 0, ∂ t V (t, x) + 1 C ∂ x I(t, x) = 0, (45)
where I(t, x) is the current intensity along the line, V (t, x) is the voltage along the line, L is the line self-inductance per unit length, C is the line capacitance per unit length.

For the circuit represented in Figure 6, the following dynamic boundary conditions are considered:

L 0 dI(t, 0) dt + R 0 I(t, 0) + V (t, 0) = U (t), (46)
C L dV (t, L) dt + V (t, L) R L -ϕ 1 (V (t, L) -V * ) + ϕ 2 (V (t, L) -V * ) 2 = I(t, L),
where R 0 is the internal resistance of the power supply, and the parallel connection of R L , C L and the CPL represent the load. Notice in the above boundary dynamics that the CPL is approximated by considering the first two terms of a Taylor expansion, with V * being the operation point and

ϕ 1 = P (V * ) 2 , ϕ 2 = P (V * ) 3 .
The operating point is defined by means of a constant input voltage U * leading to the following steady state conditions:

V * 2 - R L U * R 0 + R L V * + R 0 R L P R 0 + R L = 0, U * 2 - 4R 0 P (R 0 + R L ) R L ≥ 0, I * = U * -V * R 0 . ( 47
)
Following the modeling technique provided in [START_REF] Bastin | Stability and Boundary Stabilization of 1-D Hyperbolic Systems[END_REF], the Riemann coordinates are defined as follows

ξ 1 (t, x) (V (t, x) -V * ) + (I(t, x) -I * ) L C , ξ 2 (t, x) (V (t, x) -V * ) -(I(t, x) -I * ) L C , (48)
and thus the inverse coordinates are given by

I(t, x) =I * + ξ 1 (t, x) -ξ 2 (t, x) 2 C L , V (t, x) =V * + ξ 1 (t, x) + ξ 2 (t, x) 2 . (49)
Next, expressing the dynamics of (45) and considering Assumption 1 (see Remark 1) and the boundary conditions (46) in Riemann coordinates yields the following system of first-order Hyperbolic equations: ∂ t ξ(t, x) + Λ∂ x ξ(t, x) = 0, ξ(t, x) = ξ 1 (t, x) ξ 2 (t, x) ,

ξ 2 (t, L -x) = ξ 2 (t, x), (50)
Λ = λ 1 + δ 0 0 λ 2 + δ , λ 1 = λ 2 1 √ L C , Ẋ1 (t) Ẋ2 (t) = -α 1 -δ 0 0 -α 2 -ζ 2 X 2 X 1 (t) X 2 (t) + 0 -β 1 β 2 0 ξ 1 (t, L) ξ 2 (t, L) + 0 β 1 + δ 0 0 Û (t). ξ 1 (t, 0) ξ 2 (t, 0) = 1 0 0 1 X 1 (t) X 2 (t) + 0 1 -1 0 ξ 1 (t, L) ξ 2 (t, L) , (51)
with

α 1 = 1 L 0 L C + R 0 L 0 , α 2 = 1 C L C L + 1 R L C L - ϕ 1 C L , β 1 = 2 L 0 L C , β 2 = 2 C L C L , ζ 2 = ϕ 2 C L , Û (t) = U (t) -U * ,
where the uncertainty δ has been added to the model to cope with possible modeling errors.

For controlling the output voltage of the circuit represented in Fig. 6, consider the following control law:

Û (t) = K(ξ out )ξ out , K(ξ out ) = K 0 + 2 i=1 ξ out (i) K i , ξ out = ξ 1 (t, L) ξ 2 (t, L) , K(ξ out ) = 0 k 1 0 0 K 0 + ξ 1 (t, L) 0 0 0 0 K 1 + ξ 2 (t, L) 0 k 2 0 0 , K 2 .
with k 1 and k 2 being tuning parameters.

To evaluate the local stability of the closed-loop system, it is assumed normalized values of the system parameters, i.e., R 0 = 1, L 0 = 1, R L = 1, C L = 1, L = 1 and C ) = 1. In addition, it is considered that U * = 1 and the constant output power is set to P = 0.1 to satisfy the equilibrium conditions given in (47). Also, for this example, the domains of ξ(t, x) = ξ, X(t) and δ are supposed to be known a priori.

In order to apply the results given in Section 3, the following domains were considered D1 := ξ ∈ L 2 ((0, L); R n ) :

is defined in terms of the auxiliary variable X(t). The projection of the above set in Riemann coordinates can be easily obtained by noticing that X(t) = M ξ(t), where

M = 1 0 0 -1 0 1 1 0 , ξ(t) =     ξ 1 (t,

Concluding Remarks

This paper has focused on the boundary control of uncertain first-order hyperbolic system subject to nonlinear boundary actuation dynamics. An LMI based result has been proposed for assessing the regional robust exponential stability of the closed-loop system. The stability result is then extended to cope with nonlinear boundary control design based on strict Lyapunov functions. Two application-based examples were considered to illustrate the proposed results, where simulation results have shown the effectiveness of the proposed technique.

Figure 1 :

 1 Figure 1: Flow tube control architecture.

  is locally stable considering K = -0.8, µ = 0.3 and Ξ =[-40, 40] . It is displayed in Figure2the time response of the outlet temperature's error e L (t) considering a step change of 25 K on the outlet temperature's reference. It should be noticed the nice performance achieved by the proposed controller. Nevertheless, the result obtained by the linear controller u L (t) = -0.8ξ out can be further improved if a polytopic set Ξ out bounding the outlet temperature's trajectory ξ out (t), for all t ≥ 0, is known a priori which is the case in this example. Then, Theorem 1 with (23), K(ξ out ) = -0.8 + 8ξ out , µ = 0.3 and ξ out ∈ Ξ = [-40, 40] is applied to demonstrate that the closed loop system is locally stable with the nonlinear control law u N L = K(ξ out )ξ out . For illustrative purposes, the closed-loop response of the outlet temperature's error e N L (t) considering the latter control law is also given in Figure2. It is observed that the convergence of the temperature error with the non-linear gain outperformed the linear counterpart.

Figure 2 :Figure 3

 23 Figure 2: Temperature error trajectories e L (t) and e N L (t) for flow speed 1 m/s. Now, in order to illustrate the effectiveness of the proposed nonlinear

Figure 3 :

 3 Figure 3: Output temperature boundary control results using the proposed control (43) to the coupled PDE-ODE experimental setup.

Figure 4 :

 4 Figure 4: Control inputs for v = 1 m/s, v = 2 m/s and v = 3 m/s.

Figure 5 :

 5 Figure 5: Output temperature boundary control results for v(t) = 2 + sin(12.56 t).

Figure 6 :

 6 Figure 6: Transmission line connecting an inductive power supply to a Constant Power Load (CPL)

L 0 ξ

 0 (x) ξ(x) dx ≤ 0.2 , Ξ = X ∈ R 2 : |X i | ≤ 0.5, i = 1, 2 , ∆ = [-0.1 0.1],(52)leading to the domain of attraction estimate given in Figure7obtained from optimization problem[START_REF] Santos Martins | A multimodel approach to Saint-Venant equations: A stability study by LMIs[END_REF] with µ = 0.2. Notice that the derived estimate D2 = X ∈ R 2 : X P1 X ≤ 1 , P 1 = 4.23 0.00 0.00 4.25 ,

Figure 7 :

 7 Figure 7: Stability domain .