domain: spi phys math info

text:
C Sarotte

J Marzat

H Piet Lahanier

A Iannetti

M Galeotta

G Ordonneau

Actuator Fault Tolerant System For Cryogenic Combustion Bench Cooling Circuit

Keywords: Change detection, actuator fault accommodation, cryogenic system, fault tolerant control, fault detection and isolation

In this paper a method is proposed to design a fault tolerant control strategy based on a pressure and mass flow rate model, as well as a fault detection and isolation scheme to improve the reliability of a cryogenic bench operation. The fault is detected by the fault detection and isolation scheme composed of an extended unknown input observer, a cumulative sum algorithm and an exponentially moving average chart. Then the goal is to provide a fault tolerant system reconfiguration mechanism with a control law which compensates for the estimated actuator additive faults to maintain the overall system stability. For that purpose we use a linear quadratic regulator on an equivalent system where the unknown input is expressed as a function of the known state and known input vectors. The model and the estimation part were validated on the real Mascotte test bench (ONERA/CNES) data, and the reconfiguration control law was validated in realistic simulations.

INTRODUCTION

The need for increased launch safety and launcher's engines reliability leads to the development of health monitoring systems. The experience acquired during the years of Ariane 5 system's exploitation has pointed out the complexity of the implementation of cryogenic propulsive systems as well as the necessity to get a specialized expertise on physical phenomenon. In terms of traditional engines the objective is the improvement and the reliability of the implementation (see [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for mascotte cryogenic test bench[END_REF]). It is then necessary to be able to handle emergency situations arising from actuator failures that can affect the engine performance. The failure should be detected quickly, then isolated and its causes should be identified. For that purpose, Fault detection and isolation (FDI) methods have been developed to evaluate failures and take a decision using all available information with the help of explicit or implicit models. The most common model-based approach makes use of observers to generate residuals. Faults are detected by setting a fixed or variable threshold on each residual signal as in [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. Those FDI methods assume that the mathematical model used is representative of the system dynamics. This is challenging in practice because of the presence of modeling uncertainties and unknown disturbances. To tackle the problem of those unknown disturbances, a simple class of full order observers for linear systems with unknown inputs. It consists in a coordinate system transformation that decouples the disturbance effect on the system outputs. The observer resulting from such an approach is called Unknown input observer (UIO) (see for example [START_REF] Darouach | Fullorder observers for linear systems with unknown inputs[END_REF]).

An UIO is used to estimate the unknown state of the system independently of the unknown input. In the case of non-linear systems one of the developed techniques is to linearize and design an Extended unknown input observer (EUIO) as described in [START_REF] Witczak | Modelling and estimation strategies for fault diagnosis of non-linear systems: from analytical to soft computing approaches[END_REF]. Once the fault is detected by the FDI method based on observers, one solution to accommodate the presence of fault is to adopt a Fault tolerant control (FTC) strategy, for example see [START_REF] Cieslak | Development of an active fault-tolerant flight control strategy[END_REF]; [START_REF] Hamayun | A fault tolerant control allocation scheme with output integral sliding modes[END_REF]. The main objective of a Fault tolerant control system (FTCS) is to maintain, with a control reconfiguration mechanism, current performances close to the desirable ones and preserve stability conditions in the presence of component and/or instrument faults. An active FTCS is characterized by an on-line FDI process [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF]) which detects and estimates the fault, the second step is to achieve a steady-state tracking of the reference input by compensating the fault [START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]). In this paper a FTC method for the cooling circuit of a cryogenic combustion bench, Mascotte (CNES/ONERA), is studied. This bench has been developed to study heat transfers in the combustion chamber and jet separation in nozzles in the same conditions as for Vulcain 2 motor. The method proposed here is based on a physical model giving the output pressure and mass flow rate of the cooling circuit (Section 2). This model presenting non-linearity and unknown terms is more complete than the previous one established in [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for mascotte cryogenic test bench[END_REF]. To generate residuals we estimate the state with the help of an EUIO (Section 3) and reconstruct the unknown input with the help of a high order sliding modes observer proposed by [START_REF] Kalsi | Slidingmode observers for systems with unknown inputs: A high-gain approach[END_REF]; [START_REF] Zhu | State estimation and unknown input reconstruction via both reduced-order and high-order sliding mode observers[END_REF]. Then the residual is analyzed with a Cumulative sum (CUSUM) algorithm using an Exponentially weighted moving average (EWMA-C) chart to detect a mean shift [START_REF] Jiang | Adaptive cusum procedures with ewma-based shift estimators[END_REF]; [START_REF] Ryu | Optimal design of a cusum chart for a mean shift of unknown size[END_REF]; [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]). If an actuator failure is detected we switch to a FTCS (Section 4) to design a controller based on an unknown input observer by considering the fault to be the unknown input similar to [START_REF] Hamayun | A fault tolerant control allocation scheme with output integral sliding modes[END_REF]. This complete FDI/FTC loop permits to localize an actuator fault thanks to an estimation of the full state in spite of unmeasured information and then to compensate on-line for the fault. The results are obtained with off-line tests based on real experimental data and the reconfiguration control law was validated on realistic simulations based on the established model.

SYSTEM DESCRIPTION

The cryogenic combustion bench Mascotte (Figure 1) performs an oxygen / hydrogen operation with pressures and mass flow rates comparable to an injection element of the Vulcain 2 engine. The system studied here is its water cooling circuit. This part of the circuit feeds the combustion chamber. The detection of a leak or an obstruction is a critical safety task for the bench operation. The water cooling circuit consists in different pipes sections with multiple pressure release valves and a tank at the inlet. The available measurements are pressure, mass flow and temperature. Sections are separated by sliding valves with additional pressure measurements. The circuit between two ferrules can be modeled by two cavities defined in pressure and temperature linked by an orifice where friction forces and heat flux exchanges are taken into account, see [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for mascotte cryogenic test bench[END_REF]. A pressure regulator (actuator) permits to regulate the first cavity pressure of the water cooling circuit.

The method used here (Figure 2) consists of a loop composed of a first EUIO (unknown output mass flow rate) for fault detection purpose and another loop with a second EUIO (actuator failure) for FTC purpose. Once and only if a change has been detected by the FDI function, the FTC is included in the closed-loop system. The fault estimate is calculated from the estimated unknown mean shift provided by the ACUSUM. Then the second EUIO ensures the estimation error convergence (and stability) and the feedback controller ensures the reconfiguration error convergence (and stability).

Fig. 2. Diagram of the complete loop

The performances of the algorithms are evaluated by two means. The first one is based on off-line tests realized with real measurements of the project Conforth (CNES/ONERA). During those trials, pressures (cavities 1 and 2), temperatures (cavities 1, 2 and wall) and output mass flow rate have been recorded. The second one is based on a simulator, Carins (CNES/ONERA). This simulator permits to simulate failures for the evaluation of performances with various faults.

Model of the cooling circuit

In this section we denote ṁ the mass flow rate (kg/s), ρ the density (kg/m 3 ), S the surface (m), c the velocity of sound (m/s), u the fluid velocity (m/s), P the pressure (Pa), D the orifice diameter (m), D h the hydraulic diameter (m), L the length (m), r the radius (m), µ the dynamic viscosity (Pa.s) and V the volume (m 3 ), d t the time step (s).

The flow is assumed to stay monophasic, is ideal (no force due to viscosity acts) and incompressible following Euler equations. The cavity section is assumed constant and the velocity of sound is defined as for an isentropic reaction in the orifice. We assume that the fluid flow velocity is small in comparison to the velocity of sound. The flow crossing cavities respects the mass balance equation, after integrating this equation along the cavity length, we obtain:

∂P ∂t = c 2 V ( ṁe -ṁs ) (1)
The flow crossing the orifice between the two cavities respects the momentum balance equation with friction forces, expressed with the Darcy-Weisbach and Blasius equations for moderate turbulent flows in a smooth pipe [START_REF] Nakayama | Introduction to fluid mechanics[END_REF]). After integrating this equation along the orifice length we obtain:

1 S 2 ∂ ṁ ∂t + ∆P V = -0.316( 4 ṁ πDµ ) -1 4 L D h ṁ2 2ρV S 2
(2)

The model of this part of the cooling circuit is then:

     ∂ ṁ2e ∂t = θ 1 ṁ 7 4 2e -θ 2 ∆P ∂P 2 ∂t = -θ 3 ∆ ṁ (3) θ 1 := -0.316 4 πDµ -1 4 L D h 1 2ρV θ 2 := S 2 V , θ 3 := c 2 V
The parameter θ 1 must be identified since the distance L is unknown. We can assume here that the density and the viscosity remain constant for the considered pressures and temperature ranges. A first model of the cooling circuit has been proposed in [START_REF] Iannetti | Development of model-based fault diagnosis algorithms for mascotte cryogenic test bench[END_REF] with a constant mass flow rate.

This model (Model 1) presented approximations in the transient of the test bench inducing the presence of a pressure difference square root (see Figure 3). Moreover the mass flow rate dynamics was not modeled. The model (Model 2) presented here permits to determine the pressure but also the mass flow rate and it is now possible to model their evolution during the motor speed transients.

The model was tested off-line with real measurements of the project Conforth as inputs. Those trials last 60 seconds. The final evolution of the pressure dynamics is well reconstituted (Table 1).

L = - ρS 32µ ∆P ṁav D 2 (4)
Finally,

θ 1 := 0.316 4 πDµ -1 4 SD∆P 64V µ ṁav
For a circular pipe the hydraulic diameter D h is equal to the pipe diameter.

The set of parameters is chosen in order to fit with the measurements and in accordance with the known geometric properties of the test bench (see Table 2). Since Carins simulations and the Project Conforth do not have the same pressure reference, the parameter θ 1 varies.

FAULT DETECTION AND ISOLATION

Unknown Input Observer design

The first step is to design an observer to estimate the state in the presence of unknown inputs. The system (3), can be rewritten as a linear time-varying system with an unknown input by linearizing around a steady-state equilibrium trajectory corresponding to the mass flow trajectory to estimate the engine speed transients. Then, the system can be transformed into an equivalent discrete-time state space system with an Euler explicit scheme:

X k+1 = A k ( X)X k + BU k + ED k Y k+1 = CX k+1 (5)
where X k is the state vector, Y k is the measured output vector, U k is the known measured input vector, D k is the unknown input vector and X is the equilibrium state.

X k := ṁ2e,k P 2,k , Y k := P 2,k , U k := P 1,k , D k := ṁ2s,k
Where ṁ2e is the cavity 2 input mass flow rate, P 2 is the cavity 2 pressure, P 1 is the cavity 1 pressure and ṁ2s is the cavity 2 unknown output mass flow rate.

With A k the state matrix, B the known input distribution matrix, E the unknown input distribution matrix and C the output distribution matrix.

A k :=   1 + d t 7θ 1 m 3 4 2e,k 4 -d t θ 2 d t θ 3 1   B := [ d t θ 2 0 ] T E := [ 0 -d t θ 3 ] T C := [ 0 1 ]
The first objective is to design an observer depending only on known input and output measurements. We propose to use an EUIO with the following structure [START_REF] Witczak | Modelling and estimation strategies for fault diagnosis of non-linear systems: from analytical to soft computing approaches[END_REF]):

Z k+1 = N k+1 Z k + K k+1 Y k + GU k Xk+1 = Z k+1 + HY k+1 (6)
The above matrices are designed in such a way as to ensure unknown input decoupling as well as the minimization of the state estimate error.

e k = Xk -X k = Z k -X k + HY k (7)
To reduce its expression to a homogeneous equation we impose:

G = T B (8) T A k -N k+1 T -K k+1 C = 0 (9) T E = 0 ( 10
)
With T = I n -HC and n the dimension of the state.

A necessary condition for the existence of a solution is rank(CE) = rank(E). A particular solution is then:

H = E((CE) T (CE)) -1 (CE) T (11)
The matrix N k+1 should be Hurwitz to make the observer converge asymptotically. This is the case if we choose:

N k+1 = T A k -K k+1 C (12)
The gain matrix K k+1 is chosen to minimize the variance of the state estimation error. For a linear estimator under gaussian hypotheses (the standard deviation is denoted σ), this translates into:

K k+1 = T A k+1 P k C T (C k P k C T -R k ) -1 (13)
The covariance matrix is then given by:

P k+1 = T A k+1 P k T A T k+1 -K k+1 CP k T A T k+1 + HR k+1 H T + T Q k T T (14)
The estimation cadence used on real measurements of the project Conforth is fixed to 0.03 second, the acquisition machine acquires one point each 0.01 second and delivers information to the surveillance machine at a rate of one point each 0.03 second. The state estimation error ( 7) is taken as a residual. Figure 4 and Table 3 report the estimation results of the UIO, which are very satisfactory. Moreover in the case where it is not possible to measure the mass flow rates we can obtain an accurate estimate, in the permanent regime of the engine. The first peak in Figure 4 corresponds to the transient.

FDI with an adaptive CUSUM and EWMA-C shift estimator

The FDI mechanism is supposed to detect and diagnose any relevant failure that could lead to engine performance degradations. This shall be done sufficiently early to set up timely safe recovery. One way to proceed to detect faults is to evaluate the residual corresponding to our state estimator error. The objective is then to be able to detect a residual mean shift for a nominal behavior, see [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF].

The two hypotheses considered are then: H0: The mean value of the residual is nominal µ = µ 0 H1: The mean value of the residual has a shift µ = µ 1

In most common practical cases µ 1 is unknown. This can be overcome by using Adaptive CUSUM algorithm (ACUSUM) which estimates this value as in [START_REF] Jiang | Adaptive cusum procedures with ewma-based shift estimators[END_REF].

As mean shift amplitudes can vary drastically for a class of failure, the estimator designed for δ is a generalization of the Exponentially Weighted Moving Average (EWMA) which presents enhanced efficiency for estimation of large mean shifts: δk = δk-1 + Φ γ (e p,k )

With e p,k = r k -δk-1 the prediction error, Φ γ is defined as a Huber score function.

Φ γ :=

e p + (1 -λ)γ , e p < -γ λe p , |e p | ≤ γ e p -(1 -λ)γ , e p > γ
With γ ≥ 0, usually fixed constant. This leads to the following ACUSUM Statistic:

s k = ± δ± σ 2   r k -µ 0 ± δ± 2   ( 16
)
where for a mean shift increase or decrease: δ+ := max (δ +,min , δk ), δ-:= min (δ -,min , δk ) δ +,min and δ -,min are here the minimum mean shifts amplitudes to detect. The threshold is chosen to be a security coefficient times δ+ (Table 4). γ is defined here at each step by γ :=| r k -δk-1 | /2.

With this choice of γ, the algorithm is more efficient for the detection of small shifts. This generalization is referred to as an EWMA-C statistic. To select the coefficients values and test the algorithm performance, three faults have been simulated with Carins for different profiles of the cooling circuit inflow valves closures and openings. The objective of this FDI system is to be able to detect abrupt changes and to differentiate state perturbations and speed transients characterized by slower variations from a failure.

The first fault simulated is abrupt with a large mean shift (Figure 5), the second one has a slow variation with also a large mean shift (Figure 6) and the third one contains two faults one with a small mean shift, another one with a large mean shift (Figure 7). The first one has a slow shift then an abrupt recovery; the second one has an abrupt shift and a slow recovery. The total time of the simulation is 60 seconds with a time step of 1 millisecond (Table 5).

The cadence of the estimation and the detection is 1 time step per 30 milliseconds. The EUIO from the previous subsection permits to estimate outputs and generate the residual as the state estimator error defined by r k := Y k -C Xk . Those rates have been calculated from ten runs for each simulation and the settings have been chosen to optimize the good detection rate and minimize the false detection rate of abrupt mean shifts. Results on Fault 2 are satisfactory since it is not mandatory to detect slow variations that can be confused with transients (see Table 6). Good results are obtained for Faults 1 and 3. The last case permits to evaluate the algorithm performance for successive faults of different sizes and this is better to avoid false alarms. In some rare cases the system behavior between two faults can be considered to be faulty but in most cases the two faults in Fault 3 are well detected.

ACTIVE FTCS FOR THE INPUT PRESSURE

We consider in this part an additive actuator failure on the system. Once the fault has been detected by an online and real-time FDI mechanism the goal is to maintain the overall system stability and an acceptable performance in spite of the occurrence of faults by reconfiguring the nominal control law.

Since the system is controllable in pressure and a pressure failure induces a mass flow rate failure, it is necessary to be able to decouple the fault effect on the system dynamics from the unknown input one. The first step is then to decouple the reconstructed cavity 2 output mass flow rate from the system dynamics in order to evaluate the actuator fault and compensate it.

This method can also be useful in the case of Vulcain 2 engine during an Ariane flight, where it is difficult or expensive to measure the mass flow rate.

Unknown Input Reconstruction via high order observer

In [START_REF] Zhu | State estimation and unknown input reconstruction via both reduced-order and high-order sliding mode observers[END_REF] and [START_REF] Kalsi | Slidingmode observers for systems with unknown inputs: A high-gain approach[END_REF], an auxiliary output vector is introduced so that the observer matching condition is satisfied and is used as the new system output to asymptotically estimate the system state without suffering the influence of the unknown inputs. From this result, it is possible to build an unknown input reconstruction method based on both the state and the auxiliary output derivative estimates.

The auxiliary output is defined as:

Y i a,k := C i a,k X k
with i = 1, ..., p and p is the number of rows of Y k . The auxiliary output vector contains the output information of the original system.

If we denote:

C a,k := C 1 ... C 1 A γ1-1 k ... C p ... C p A γp-1 k T with 1 ≤ γ i ≤ n i i = 1, ..., p
where n i is defined as the smallest integer such that:

c i A γi k E = 0 γ i = 0, 1, ..., n i -2 c i A ni-1 k E = 0 (17)
and

C i the i-th row of C. C i a,k := C i ... C i A γi-1 k T
Since the auxiliary output vector depends on unmeasured variables, we can design a high-order sliding mode observer to get the estimates of both the auxiliary output vector and its derivative.

After discretization we have:

Y i a,k+1 = C i a,k+1 (A k X k + BU k + ED k ) (18)
If we denote:

Λ i := 0 I γi-1 0 0 , r i := 0 (γi-1)×1 1 , Ψ i k := C i a,k B
Then (18) can be written as:

Y i a,k+1 = Λ i Y i a,k + r i f i k (X k , D k ) + Ψ i k U k (19)
where

f i k (X k , D k ) := C i A γi-1 k (A k X k + ED k )
The last equation of this n i size system is:

C i A γi-1 k ED k = Y γi a,k+1 -C i A γi-1 k (A k X k + BU k ) (20)
The above p equations can be unified into a single matrix:

M k D k = ξ k+1 -C k (A k X k + BU k ) (21)
if we denote

M k := C k E C k := (C 1 A γ1-1 k ) T (C 2 A γ2-1 k ) T . . . (C p A γp-1 k ) T T ξ k+1 := (Y γ1 a,k+1 ) T (Y γ2 a,k+1 ) T . . . (Y γp a,k+1 ) T T Since rank(M k ) = rank(C a,k D k ) = rank(D k ) = q, M T k M k is invertible because M k
has full column rank. So the input vector satisfies:

D k = (M T k M k ) -1 M T k (ξ k+1 -C k (A k X k + BU k )) (22)
An estimation of it is then:

Dk = (M T k M k ) -1 M T k ( ξk+1 -C k (A k Xk + BU k )) (23) with ξk+1 :=     C 1 A γ1+1 k Xk + C 1 A γ1-1 k BU k C 2 A γ2+1 k Xk + C 2 A γ2-1 k BU k . . . C p A γp+1 k Xk + C p A γp-1 k BU k    
To validate the result, the unknown input reconstruction is compared to the cavity 2 output mass flow rate measurements available for this trial. Results are reported in Figure 8 and Table 7 and show a correct convergence after the transient phase. Now that the unknown input expression is available, we can rewrite the first system without the mass flow rate as an unknown input. Then, in order to annihilate the actuator fault effect on the system, another EUIO is used to estimate the fault magnitude. A control law has then to compensate the fault and be computed such that the faulty system is as close as possible to the nominal one.

We can use the previous result ( 22) from the unknown input reconstruction part to rewrite the system under a second form only depending on known inputs for control purposes:

X k+1 = AX k + BU k + E((M T M ) -1 M T (ξ k+1 -C(AX k + BU k ))) (24)
X k+1 = (I n -E(M T M ) -1 M T CA) -1 [(A -E(M T M ) -1 M T CA)X k + (B -E(M T M ) -1 M T CB)U k ] ( 25
)
The system is linearized around a steady state equilibrium, the nominal state to reach, the matrix A is then constant in time.

This method requires matrix inversions, which may be numerically unstable due to possible ill-conditioning. In the problem considered, the matrices were invertible.

We obtain a new system under the form:

X k+1 = A c X k + B c U k Y k+1 = CX k+1 (26)
with A c the new state matrix and B c the new known input distribution matrix:

A c := (I n -E(M T M ) -1 M T CA) -1 (A -E((M T M ) -1 M T CA)) B c := (I n -E(M T M ) -1 M T CA) -1 (B -E((M T M ) -1 M T CB))
The estimate of the state is given by: Xc,k+1 = η k+1 + e c,k+1 + X k+1

An additive actuator failure with a control law can be modeled as:

X k+1 = A c X k + B c U n,k + B c (f k + U c,k ) Y k+1 = CX k+1 (28)
where we assume U n,k , the nominal input, to be known, U c,k is the control law and f k is the faulty part of the input. U n,k and f k are components of U k .

It is then possible to design a second unknown input observer for the reconfiguration part, where f k + U c,k is considered to be the unknown input, with the following structure [START_REF] Hamayun | A fault tolerant control allocation scheme with output integral sliding modes[END_REF]):

Z c,k+1 = N c,k+1 Z c,k + K c,k+1 Y k Xc,k+1 = Z c,k+1 + H c Y k+1 (29)
The above matrices are designed in such a way as to ensure unknown input decoupling from the estimation error dynamic as well as the minimization of the state estimator error variance as previously.

e c,k = Xc,k -X k = Z c,k -X k + H c Y k (30)
To reduce this expression to a homogeneous equation we impose:

H c = B c ((CB c ) T (CB c )) -1 (CB c ) T (31) N c,k = T c A c -K c,k C (32)
To give the state estimator error the minimum variance, the gain matrix should be determined to minimize the covariance matrix:

K c,k = (T A) c P k C T (CP k C T -R k ) -1 (33)
The EUIO stability is adressed in [START_REF] Witczak | Modelling and estimation strategies for fault diagnosis of non-linear systems: from analytical to soft computing approaches[END_REF].

We also have to ensure the convergence of the regulation error η k .

η k+1 = X k+1 -X k+1 (34)
η k+1 = A c (X k -X k ) + B c f k + B c U c,k (35) η k+1 = A c η k + B c f k + B c U c,k (36)
Then we can use a control law of the form:

U c,k := -B + c B c fk -W c ( Xc,k -X k )
where -B + c B c fk is the fault compensation part and W c ( Xc,k -X k ) is the reconfiguration part. The fault magnitude estimation fk is obtained from the fault diagnosis part (Section 3.2 Residual) and unknown input reconstruction (Section 4.1 Unknown input estimate). Which can be alternatively written as

U c,k := -B + c B c fk -W c e c,k -W c η k
Where we denote B + c as the pseudo-inverse of B c [START_REF] Theilliol | Actuator fault tolerant control design based on a reconfigurable reference input[END_REF]). Then we have:

ζ k+1 = A c -B c W c -B c W c 0 N c ζ k where ζ k := [ η k e c,k ] T .
For the nominal system, the gain W c must stabilize (A c -B c W c ). Since the pair (A c , B c ) is assumed to be controllable, a Linear-Quadratic Regulator (LQR) formulation can be adopted where W c is selected to minimize

J k := k X T k Q c,k X k + U T c,k R c,k U c,k
Where Q c,k and R c,k are symmetric positive definite design matrices. It is also possible to proceed to a pole placement for the continous time system (small time constant), we can choose to fix a damping ratio and a natural frequency. The desired transient behavior depends on the gain choice (Table 8), in our case we have to limit the overshoots to maintain the cooling circuit performances (Figure 9). The fault was implemented as in the previous section. The aim of this simulation is to see if the controller is able to stabilize the closed-loop system after the detection, see Table 9. If a fault is detected then the system switches to the closed-loop one. When the fault is detected the system switches to the FTCS. The fault is compensated and it can be seen that the control law for the rewritten system permits to stabilize the system around the reference steady-state equilibrium with sufficient precision.

CONCLUSION

In this paper a new model was proposed for the evolution of pressure and mass flow rates in the cooling circuit of a cryogenic test bench. Once the fault in the actuator has been detected by the FDI method composed of a first EUIO, the designed FTCS based on a fault estimator and a second EUIO permits to compensate the failure and to converge if necessary to a chosen steady state. This FTCS consist in a linear quadratic regulator on an equivalent system where the unknown input is expressed as a function of the known state and known input vectors in order to decouple only the fault effect on the system. Future work will address the design of a method to calculate another steady point which may be reachable in the case where the previous nominal steady point cannot be reached because of the actuator failure and the effect of the saturation. Being able to choose the nominal behavior of the system may then be useful to take into account actuator saturation.

Fig. 1 .

 1 Fig. 1. Mascotte test bench -Ferrules (ONERA -CNES)

Fig. 3 .

 3 Fig. 3. Pressure model -Project Conforth -Real data

Fig. 4 .

 4 Fig. 4. Pressure residual -Project Conforth

Fig. 5 .

 5 Fig. 5. Fault 1 residual -Carins simulation

Fig. 8 .

 8 Fig. 8. Unknown input reconstruction -Project Conforth

Fig. 9 .

 9 Fig. 9. Pressure and mass flow rate control -Carins simulator -ATAC configuration

Table 1 .

 1 Relative errors of the pressure model 1 and 2

	Model	Total	Transient Permanent
		(%)	(%)	(%)
	Pressure (1)	10.58	13.35	5.04
	Pressure (2)	5.44	8.01	0.31
	Input mass flow rate (2) 3.31e-5	4.97e-5	6.17e-8

Table 2 .

 2 Parameters for estimation and control

	Parameters	Carins	Mascotte -Conforth
	θ 1	-0.1858	-0.2864
	θ 2	6.07	6.07
	θ 3	2.25	2.25

Table 3 .

 3 Relative errors of the pressure and input mass flow rate estimations

	Model	Total	Transient Permanent
		(%)	(%)	(%)
	Pressure (Pa)	9.92e-2	7.00e-2	1.16e-2
	Mass flow rate (kg/s)	6.27	31.4	1.18e-2

Table 4 .

 4 Parameters for fault detection

	Parameters	Values
	δ +,min	4e-2
	δ -,min	-4e-2
	λ	1.1055
	Threshold coefficient	4.5e4

Table 5 .

 5 Simulated faults

					end
	Fault 1	41010	46200	1367	1540
	Fault 2	30960	37560	1032	1252
	Fault 3(1)	39300	41040	1310	1368
	Fault 3(2)	45960	60000	1532	2000

Fault

t begin (s) t end (s) N begin N

Table 6 .

 6 Detection rates

	Fault	GDR	FDR
	Fault 1 98.8%	0.0%
	Fault 2 27.4%	0.0%
	Fault 3 98.5%	14.5%

Table 7 .

 7 Relative error of the output mass flow rate reconstruction

	Model	Total Transient Permanent
		(%)	(%)	(%)
	Output mass flow rate (kg/s)	17.6	35.2	1.19e-2
	4.2 Fault tolerant control system design

Table 8 .

 8 Gain matrix choice

	Pressure Mass flow Damping	Natural
	part	rate part	ratio	frequency
	-0.3668	-0.9956	2	1e-1

Table 9 .

 9 Relative errors of the simulated pressure and input mass flow rate from references

	Model (2) with control	Permanent (from detection time)
		(%)
	Pressure (Pa)	1.6e-1
	Input mass flow rate (kg/s)	7e-2


HAL is

Deepening analysis of uncertain categorical inputs using

Gaussian processes -application to marine flooding J. Rohmer 1 , D. Idier 1 , S. Lecacheux 1 , R. Pedreros 1

[1]{BRGM, 3 av. C. Guillemin -BP 36009 -45060 Orléans Cédex 2 -France} Correspondence to: J. Rohmer (j.rohmer@brgm.fr)

Marine flooding assessment relies on the use of high-resolution hydrodynamic numerical models to understand the relationships between inputs (e.g. offshore sea conditions) and the output variables (e.g. flood spatial extent), and to characterize the uncertainty. The inputoutput mapping is, however, hindered by two aspects: 1. the high computation time cost of the numerical model limits the number of model runs (usually a few 100s); 2. some of the input variables are often categorical in nature, because different modelling assumptions are equally appropriate. In this communication, we propose to explore how recent advances in Gaussian Processes with mixed continuous/categorical inputs can bring valuable insights into the influence of such categorical inputs. Depending on how they are modelled (nominal, ordinal) and integrated in the kernel (level-related heteroscedasticity, compound symmetry, group structure and interactions, etc), different viewpoints can be brought. Three application cases with gradual complexity are investigated: i. nominal input to analyse the scenarios of rupture geometry for tsunamis; ii. ordinal variable to reveal the latent influence of the angle of approach for cyclone tracks; iii. group kernels to reveal the interactions between artificial and natural coastal defences' failure scenarios for storm tides. Finally, the problem of selecting the most appropriate kernel structure using exploratory analysis and cross-validation procedures is discussed.


text:
P Apkarian
email: pierre.apkarian@onera.fr

D Noll

Boundary control of partial differential equations using frequency domain optimization techniques

Keywords: math Boundary control of PDEs, frequency-domain design, convection-diffusion, wave equation, H ∞, structured feedback, infinite-dimensional systems

We present a frequency domain based H ∞ -control strategy to solve boundary control problems for systems governed by parabolic or hyperbolic partial differential equation, where controllers are constrained to be physically implementable and of simple structure suited for practical applications. The efficiency of our technique is demonstrated by controlling a reaction-diffusion equation with input delay, and a wave equation with boundary anti-damping.

Introduction

A recurrent issue in system control is whether, or to what extent, frequency-domain based H ∞ -control strategies originally developed for real-rational systems expand to infinite-dimensional processes. Success in rendering H ∞optimization fit to provide practically implementable controllers for infinite-dimensional systems should substantially foster the acceptance of PDE-modeling as a tool for control.

In response to this quest, we present a frequency domain based method to control infinite-dimensional LTIsystems, which is in particular suited for H ∞ -boundary control of parabolic and hyperbolic partial differential equations. Our method leads to practically implementable structured output feedback controllers for PDEs in such a way that the typical work-flow in control design is respected.

After briefly outlining our method, we will apply it in more detail to two infinite-dimensional H ∞ -control problems: boundary control of a reaction-diffusion system with input delay, as discussed in [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], and boundary control of an anti-stable wave equation to control noise and disturbance effects on duct combustion dynamics in a drilling pipe system [2, p.6], [START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF]. While the first study leads to a parabolic equation of retarded type, the second study leads to a system of neutral type, which poses new challenges to our frequency approach.

The structure of the paper is as follows. In section 2 we give the principal steps of our method. Stability is discussed in section 3, the role of the Nyquist test in optimization in section 3.1, its implementation in section 3.2. Sampling for performance is addressed in section 4. In section 6, we discuss the application of our method to a reaction-diffusion equation, and in section 7 to a wave equation.

Outline of the method

We start out with an infinite-dimensional LTI-system represented by a transfer function G(s) with p inputs and m outputs, assumed well-posed in the sense of Salamon-Weiss [START_REF] Salamon | Infinite dimensional linear systems with unbounded control and observation: a functional analytic approach[END_REF][START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF][START_REF] Chang | Well-posedness of boundary control systems[END_REF][START_REF] Curtain | Well posedness of triples of operators (in the sense of linear system theory)[END_REF]. As principal application, we consider the case of a linearized parabolic or hyperbolic boundary control problem in state-space form

Γ :          ẋ = Ax Px = u y = Cx (1)
with operators A ∈ L(X, H), P ∈ L(X, R p ), C ∈ L(X, R m ) on Hilbert spaces H, X and finite-dimensional input and output spaces, where X is densely embedded in H. Then under natural assumptions specified in [START_REF] Chang | Well-posedness of boundary control systems[END_REF]Sect. 2] the transfer function G(s) of ( 1) is well-posed and obtained by applying the Laplace transform to [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], where every function evaluation G(s) requires solving an elliptic boundary control problem

Γ s :          sx(s) = Ax(s) Px(s) = u(s) y(s) = Cx(s) (2)
Well-posedness means that G(s) is holomorphic on a half-plane Re(s) > σ, but it may be convenient to require a little more, namely, that G(s) extends meromorphically over a domain containing C + . This is satisfied in all cases of practical interest, and guaranteed theoretically e.g. when G is exponentially input/output stabilizable, see [8, Lemma 8.2.9 (i)(b), (ii)]. The meromorphic form of the transfer function is a necessary requirement for applicability of the Nyquist stability test.

After embedding G(s) in a plant P(s) with one or several closed-loop performance and robustness channels T wz (P, K), we set up the infinite-dimensional H ∞ -optimization problem minimize T wz (P, K) ∞ subject to K stabilizes G in closed loop

K ∈ K (3)
where K represents a suitably chosen class of structured controllers with m inputs and p outputs. In this work we understand the term structured in the sense that controllers K(x, s) depend differentiably on a vector x ∈ R n of tunable parameters, and have well-posed transfer functions K i j (x, s), typically with quasi-polynomial numerators and denominators. Such control laws combine real-rational elements with input and output delays, and can therefore be physically implemented. In the optimization procedure it will also be necessary to know the finite number of rhp poles of K(x) for every x. Since (3) as a rule cannot be solved exactly, we use an inexact bundle trust-region method as in [START_REF] Noll | Bundle method for non-convex minimization with inexact subgradients and function values, Computational and Analytical Mathematics[END_REF][START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF][START_REF] Apkarian | Non-smooth optimization for robust control of infinite-dimensional systems[END_REF], which guarantees stability of the closed loop, and approximates H ∞ -performance up to a user specified precision. The following scheme presents our method in a more formal way.

Algorithm 1. H ∞ -control of infinite-dimensional systems

Parameters: ϑ > 0.

Step 1 (Prepare). Linearize system about steady-state and pre-compute open-loop transfer function G(s).

Step 2 (Initialize). Choose controller structure and find initial closed-loop stabilizing controller K(x 0 ) of that structure. Let G 0 = feedback(G, K(x 0 )).

Step 3 (Plant). Embed G 0 into plant P representing desired closed-loop performance specifications.

Step 4 (Non-smooth optimization). Run inexact bundle trust-region method [START_REF] Noll | Bundle method for non-convex minimization with inexact subgradients and function values, Computational and Analytical Mathematics[END_REF] with starting point x 0 , discretizing (3) at each iterate x j so that Nyquist test guarantees stability of the loop, and H ∞ -performance up to tolerance ϑ.

Several steps of this scheme require further explanations, which we provide in the following sections.

Stability test

Let us recall that with the definitions

F(s) = I G(s) -K(s) I , f (s) = det F(s)
the inverse T (s) = F(s) -1 is given as

T = (I + KG) -1 -K (I + GK) -1 (I + GK) -1 G (I + GK) -1 , (4)
and we call the closed-loop system (G, K) stable in the H ∞ -sense, or simply stable, if the transfer function T belong to the Hardy space H ∞ (C + , C (m+p)×(m+p) ). As is well-known, H ∞ -stability is equivalent to the absence of unstable poles in tandem with boundedness of T (s) on jR. We are interested in situations, where absence of unstable poles of T (s) can be verified by the Nyquist stability test. Systems arising from parabolic equations are of retarded type and typically satisfy the spectrum decomposition assumption, which means that they have only a finite number of unstable poles. The Nyquist stability test may therefore be applied directly, to the effect that in order to guarantee absence of unstable poles in the loop f ( jω) = det(I + G( jω)K(x 0 , jω)) has to wind n p times around the origin in the clockwise sense, where n p the number of rhp poles of G and K(x) together.

In order to address the case where f has a finite number of poles on jR, we consider the following construction, which avoids the usual -indentation of the Nyquist contour into the rhp. We choose a holomorphic function h on a domain containing C + such that h(s) 0 on C + , lim s→∞ h(s) = 1 on C + , and such that h has a zero of order p at ± jω precisely when F has a pole of order p at ± jω. Let f = f h, D a Nyquist D-contour into the rhp with [jω, jω] ⊂ D containing in its interior all rhp poles of F. Then the modified Nyquist curve f • D has the same winding number as the original Nyquist curve f • D with sufficiently small -indentations.

From the moment onward a controller K(x) has been identified closed-loop stabilizing using the Nyquist test, the nonsmooth optimization method, when considering a trial step K(x + dx) away from the current iterate x, will re-compute the winding number to check stability of the the loop with K(x + dx). In those cases where the number of poles of K(x) is independent of x, this means we simply have to assure that the winding number n p does not change as we go from x to x + dx, which requires preventing the Nyquist curve from crossing the origin. Stability in the H ∞ -sense for the new K(x + dx) will then follow under the proviso that the closed-loop transfer function T (x + dx) remains bounded on jR, which is the case when G, K(x) are bounded on jR, and occurs in particular if these transfer functions are proper. With these preparations the situation for parabolic systems is covered by the following: Theorem 1. Suppose process G(s) and controller K(s) are well-posed, extend meromorphically into a domain containing C + , and satisfy the following conditions: i. F has no zeros on jR and only a finite number N p of poles in C + , n p of which are in C + . ii. There are no pole/zero cancellations on jR. iii. There exist a frequency ω > 0 and α > 0 such that Re f ( jω) > α for all ω ∈ [ω, ∞). iv. G, K are bounded on jR \ [jω, jω]. v. G, K have strongly (exponentially) stabilizable and detectable state-space realizations.

Suppose the modified Nyquist curve f •D winds n p times around the origin in the clockwise sense. Then the closed-loop T (s) is strongly (exponentially) stable.

The theorem was proved in [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF] under the more standard assumptions that G, K are proper and lim s→∞ f (s) 0 on C + exists, but since our present statement concerns only a finite Nyquist contour, f • D, the proof can be adapted with minor changes.

The full force of Theorem 1 is needed when it comes to dealing with hyperbolic systems. Here the situation is complicated because in open loop neutral systems may have infinitely many rhp poles in a strip 0 ≤ Re(s) < α, so condition (iii) may fail, even though G may still be well-posed. When this is the case, it is impossible to use the Nyquist test directly even if an initial stabilizing controller K 0 = K(x 0 ) is given, because the Nyquist curve f (s) = det(I + G(s)K 0 (s)) winds infinitely often around the origin. In this event the method explained in the following section is helpful.

Enabling the Nyquist test

Assuming that is has been verified by some other means that the closed-loop system G 0 = G(I + K(x 0 )G) -1 =: feedback(G, K(x 0 )) is indeed stable. Then a small gain argument tells us that feedback(G 0 , K) remains stable for stable controllers K satisfying K ∞ < 1/ G 0 ∞ . Since feedback(G 0 , K) = feedback(feedback(G, K(x 0 )), K) = feedback(G, K(x 0 ) + K), we see by letting K the difference K(x) -K(x 0 ) that feedback(G, K(x)) is stable for K(x) -K(x 0 ) ∞ < 1/ G 0 , and this can now be verified by applying the Nyquist test to f 0 (s) = det(I + G 0 (s)K(s)), where K = K(x) -K(x 0 ).

Here Theorem 1 applies indeed to G 0 , K, because G 0 has no poles on C + , so that i.-iv. are satisfied for G 0 , K provided iii. was from start satisfied for G, K 0 . We have, however, to recall that despite stability of the loop, f 0 will typically not have a limit as s → ∞ on C + , so we will rely on condition iv., which assures that outside the band [-ω, ω] the Nyquist curve is in no danger to turning around 0.

From here on, we can proceed just as in the previous case for retarded systems, where now N p = n p is the known number of rhp poles of K = K(x) -K(x 0 ). Remark 1. Should our initial stabilizing controller K 0 be unstable, it may be preferable to use controllers of the form K 0 + K(x) for K(x) stable, which gives rise to a modified structure.

The question remains in what sense a state-space representation of T given by (4) will be stable. This is decided by the following Theorem 2. Suppose G, K are well-posed transfer functions which admit strongly (exponentially) stabilizable statespace realizations. Suppose the closed loop transfer function satisfies T ∈ H ∞ . Then the generator of the state-space representation of the closed loop is strongly (exponentially) stable.

Proof. For exponential stability, according to Morris [START_REF] Morris | Justification of input-output methods for systems with unbounded control and observation[END_REF]Theorem 5.2] it suffices to show that the closed loop is exponentially stabilizable and exponentially detectable. By Staffans [START_REF] Staffans | Well-Posed Linear Systems, Encyclopedia of Mathematics and its Applications[END_REF]Lemma 8.2.7] this follows as soon as each of the components G, K is individually exponentially stabilizable and detectable, but the latter is true by hypothesis.

For the statement concerning strong stability we use [8, Lemma 8.2.7] again, which now guarantees that the closed loop T (G, K) is strongly stabilizable and detectable. Since it is H ∞ -stable, we can invoke [START_REF] Staffans | Well-Posed Linear Systems, Encyclopedia of Mathematics and its Applications[END_REF]Theorem 8.2.11 (ii)] to infer that the closed loop is also strongly stable.

⋆ + ⋆ = 0, + = 1

Sampling for the Nyquist test

During optimization, the Nyquist test is applied at every candidate controller K(x) to check stability of the loop feedback(G, K(x)), respectively for neutral systems, of the preliminary stabilized system feedback(G 0 , K(x) -K(x 0 )). For that we have to sample the modified Nyquist curve f = f h at frequencies 0

= ω 0 < • • • < ω N = ω 4
such that the winding number of f • D, based on a typical D-contour extending into the rhp, coincides with the winding number of the closed polygon P f obtained as f ( jω 0 ), . . . , f ( jω N ), f (jω N ), f (jω N-1 ), . . . , f (jω 0 ). The latter can then be computed conveniently using the ray crossing algorithm (see Fig. 1) To begin with, we need a cutoff frequency ω > 0 such that every Nyquist contour D containing the segment [jω, jω] contains all n p unstable poles of G, K in its interior. We then sample on the segment [0, ω], and replace the curved part of the D-contour f

• C by the segment [ f (-jω), f ( jω)].
For the latter to be authorized, the closed curve γ concatenated by this segment and f •

C

has to satisfy ind(γ, 0) = 0 and should not contain any of the n p poles of F in its interior. This is for instance the case if we choose ω as in condition iv. of Theorem 1. Then we sample f on [0, ω] such that none of the curves γ i obtained by concatenating the segment [ f ( jω i+1 ), f ( jω i )] with the piece f ([ jω i , jω i+1 ]) of the modified Nyquist curve encircles 0. This can be assured by the method in [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF]. Call a mapping

L[•, •] : R 2 → R a first-order bound of f if | f ( jω)| ≤ L[ω -, ω + ] for all ω ∈ [ω -, ω + ].
Then we have the following Lemma 1. Suppose that for fixed K the cutoff frequency ω > 0 is as above, and the sampling nodes satisfy the condition

L[ω i , ω i+1 ](ω i+1 -ω i ) < | f ( jω i )| + | f ( jω i+1 )|. ( 5
)
Then the winding numbers of the modified Nyquist curve f • D and the approximating closed polygon P f agree.

For the proof see [12, sect 3.].

Sampling for performance

Sampling for H ∞ -performance was also analyzed in [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF] and can again be based on a first-order bound L[•, •], now for the function φ(ω) = σ (T wz (P( jω), K( jω))). We recall the following Lemma 2. Let γ * = max{φ(ω i ) : i = 1, . . . , N} for a given controller K and a corresponding sampling ω i . Let ϑ > 0 be a user-specified tolerance. If the nodes ω i satisfy

L[ω i , ω i+1 ](ω i+1 -ω i ) < 2γ * + 2ϑ -φ(ω i ) -φ(ω i+1 ), (6)
then the true H ∞ norm is within tolerance ϑ of its estimated value, that is, γ * ≤ T wz (P, K) ∞ ≤ γ * + ϑ.

For the proof see [12, sect. 5]. In that work we have compared sampling for stability via [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF] and sampling for H ∞ -performance based on (6) on a large test bench including finite and infinite dimensional systems. The results fairly consistently show that performance requires at least 10 times more nodes ω i than sampling to assure that the Nyquist test is correct. This leads to the following significative meta-theorem: H ∞ -performance is 10 times more costly than mere stability.

In those cases where a channel T = T wz (P, K) for H 2 -optimization is available, we need a sampling ω i for f (ω) = trace T ( jω) H T ( jω) . Assume that a first-order bound L[ω -, ω + ] for f is available, and let P be the piecewise linear function corresponding to the polygon with nodes (ω i , f

(ω i )), 0 = ω 0 < • • • < ω N = ω with P(ω) = 0 for ω > ω.
Lemma 3. Let ϑ > 0 be a user specified tolerance and suppose the cutoff frequency ω > 0 is such that

e 1 = ∞ ω f (ω)dω ≤ ϑ/2. Suppose the interval [0, ω] is sampled with nodes ω i such that ω 4 (ω i+1 -ω i )L[ω i , ω i+1 ] ≤ ϑ/2, (7)
then the error satisfies

e = ∞ 0 f (ω)dω - ∞ 0 P(ω)dω < ϑ.
Proof. Let e 1 be the error of the high frequency contribution satisfying e 1 < ϑ/2. Now the error of the low frequency part is

e 1 = ω 0 f (ω)dω - ω 0 P(ω)dω ≤ N-1 i=0 1 4 (ω i+1 -ω i ) 2 L[ω i , ω i+1 ] ≤ N-1 i=0 (ω i+1 -ω i )ϑ/2ω = ϑ/2.
Hence altogether e = e 1 + e 2 < ϑ.

Remark 2. It is clear that [START_REF] Curtain | Well posedness of triples of operators (in the sense of linear system theory)[END_REF] is much more binding than [START_REF] Chang | Well-posedness of boundary control systems[END_REF], because sampling to assure the exactness of the maximum value within a tolerance ϑ > 0 has only to be precise at frequencies close to the maximum, whereas approximating an integral requires good approximation on the whole [0, ω]. This suggests avoiding H 2 -optimization if possible. Since robustness requirements further press to avoid H 2 -optimization, we presently seek for workarounds.

Optimization

With a computable test for stability and a method to approximate the objective function T wz ∞ available, we run our nonsmooth optimization method based on [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF][START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], with the interpretation of inexact function and subgradient evaluations as in [START_REF] Noll | Bundle method for non-convex minimization with inexact subgradients and function values, Computational and Analytical Mathematics[END_REF]. For this we have to recall subgradient evaluation as discussed in [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF]. Suppose ω i is one of the sample frequencies where the maximum γ * = φ(ω i ) of the approximation is attained with error γ * ≤ T wz (x) ∞ ≤ γ * + ϑ. Then an approximate subgradient is generated by computing one or several maximum eigenvectors Φ i of T wz (x, ω i ) H T wz (x, ω i ) and using formulas [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF] or ( 13) of [14, Sect. IV]. Inspecting those shows that for the nearby frequency ω where T wz (x) ∞ is in reality attained, the mismatch between the estimated eigenvector Φ i and one of the true maximum eigenvectors Φ of the H ∞ -norm at ω is proportional to the error in the function values, but with proportionality constant depending on the reciprocal of the eigenvalue gap at T wz (x, ω). Since only finitely many frequencies are active, the eigenvalue gap cannot become arbitrarily small. From the same reason it remains bounded in the neighborhood of any of the accumulation points of the sequence of serious iterates generated by the bundle or the bundle trust-region method. That suggests indeed an interpretation of our method as an instance of the inexact bundle trust-region method.

For a recent thorough convergence analysis of the bundle method with inexact function and subgradient evaluations in an infinite-dimensional setting we refer to Hertlein and Ulbrich [START_REF] Hertlein | An inexact bundle algorithm for nonconvex nondifferentiable functions in Hilbert space[END_REF].

Starting the algorithm as presented in [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF][START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF] at the closed-loop stabilizing K 0 = K(x 0 ), the Nyquist test is used at every new iterate K(x + ) to check whether the loop feedback(G, K(x + )), respectively feedback(G 0 , K(x + ) -K(x)) for the case of neutral systems, is stable. If this is not the case, a backtracking step x α = x + α(x +x) for 0 < α < 1 is made such that K(x α ) is still stabilizing, and a repelling cutting plane is included in the bundle, using e.g. the closed-loop sensitivity function S (x) = (I + GK(x)) -1 ∞ as a stability barrier function. See [12, sect. 4] for details.

Remark 3. A special situation occurs if G, K are stable, or if G 0 instead of G is used as described in section 3.1. Here the Nyquist curve f ( jω) does not at all wind around the origin, and as a rule stays outside a conical or parabolic region R α,r = {s ∈ C : Im(s) 2 < α(r -Re(s))} for certain α, r > 0. The constraint Im( f (x, jω)) 2 ≥ α(r -Re( f (x, jω))) can be easily included in the constrained program (3). This is more reliable than just preventing f (x, jω) from crossing the origin. Note that even when the Nyquist test can be applied directly to G, K(x), e.g. for parabolic or first-order hyperbolic systems, it may for numerical reasons be interesting to apply it in just the same way to G 0 , K.

A more conservative but robust way to address the same problem is to add a disk-margin constraint of the form S ∞ ≤ 1/α.

Output-feedback control of a reaction-diffusion equation with input delay

In this section we discuss a one-dimensional reaction diffusion equation with delayed Dirichlet boundary control

x t (ξ, t) = x ξξ (ξ, t) + c(ξ)x(ξ, t), t ≥ 0, ξ ∈ [0, L], x(0, t) = 0, (8)
x(L, t) = u(t -D)

where x(•, t) denotes the state of the system, u(t) the control, D the delay, and where we assume that a finite number of measured outputs y 1 (t) = x(ξ 1 , t), . . . , y m (t) = x(ξ m , t) at sensor positions ξ i ∈ [0, L] are available for control.

A similar control scenario is discussed in Prieur and Trélat [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF] under the assumption of full state measurement. Related work is for instance Sano [START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF], where H ∞ -control of a heat exchanger is discussed, or [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF], where a reactionconvection-diffusion equation with simultaneous boundary and distributed control and a van de Vusse reactor of a coupled system of reaction convection-diffusion equations again with combined boundary and distributed control are discussed without input delay, but with a single point measurement as output. 6

In the present study, we strive to control the system with a finite-dimensional output feedback controller K(x) of simple structure, which could conveniently be implemented, and yet gives satisfactory performances in closed loop.

Performance specifications of the reaction-diffusion PDE are chosen so that responses to non-zero initial conditions show reasonable behavior in terms of damping and settling time. This could be addressed by an H 2 -performance specification, but as we show may also be successfully controlled by way of suitably chosen H ∞ -specification. The latter is advantageous as soon as additional robustness aspects of the design are called for.

Working incrementally, and starting with the case of a single measurement at the mid-point ξ = L/2, our analysis indicates that 5 equidistant measurements are enough to achieve good responses against initial conditions, while mere stability could be assured even on the basis of a single measurement e.g. at ξ = 0.

In our numerical testing we adopt the choices L = 2π, D = 1 and c(x) = 1 2 from [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], where the open-loop F(s) has one unstable pole at s = 1 2 , and an infinity of stable double poles at

s k = 1 2 -k 2 π 2
L 2 following a retarded pattern. This is understood, as the semi-group of the equation is sectorial [8, p. 150]. As indicated in previous chapters, this allows direct application of the Nyquist test to check H ∞ -stability of the closed loop.

Remark 4. Well-posedness of the system (8) in the sense of Salamon-Weiss can be deduced from the functional analytic setting in [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], or from the general approach in [START_REF] Chang | Well-posedness of boundary control systems[END_REF].

The transfer function of ( 8) can be computed analytically as

G(s, ξ) = x(ξ, s) u(s) = e -s e √ s-1 2 ξ -e - √ s-1 2 ξ e √ s-1 2 L -e - √ s-1 2 L
, and for L = 2π the system has one unstable pole.

Proposition 1. Suppose a finite-dimensional structured controller K(x) with m inputs and p = 1 output is found which stabilizes system (8) internally in the H ∞ -sense. Then the closed loop is even exponentially stable.

Proof. By theorem 2 this follows as soon as each of the components G, K is individually exponentially stabilizable and detectable. Since K is finite-dimensional it has clearly an exponentially stabilizable and detectable state-space model. For G exponential stabilizability may be deduced from [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], because the infinite-dimensional state-feedback controller the authors construct has the same control input as [START_REF] Staffans | Well-Posed Linear Systems, Encyclopedia of Mathematics and its Applications[END_REF]. Exponential detectability on the other hand follows from the fact that the differential operator in ( 8) is self adjoint, and that in the adjoint system the five outputs are turned into 5 inputs, one of which is the same as the single input in [START_REF] Staffans | Well-Posed Linear Systems, Encyclopedia of Mathematics and its Applications[END_REF], but now without the delay. Exponential detectability therefore follows from the fact that ( 8) is exponentially stabilizable without the delay.

This still leaves the problem of finding a preliminary stabilizing controller K 0 = K(x 0 ) of the pre-defined structure. As we indicated in previous sections, the latter as a rule requires heuristic methods even for very simple structures. The advantage we have in the case of the present parabolic study [START_REF] Staffans | Well-Posed Linear Systems, Encyclopedia of Mathematics and its Applications[END_REF] is that we can check via the Nyquist test whether a given controller is stabilizing.

Model matching approach

Model matching is a sophisticated control scenario, where specifications are pursued indirectly. It is covered by algorithm 1, and we believe it is particularly suited for PDE-control, where models of different grid scales arise naturally. Here we use model matching to address the reaction of the system to a non-zero initial value. The method consists in two steps (a) and (b).

Non-vanishing initial values may be regarded as disturbances d acting on the system state as in Fig. 2. We now assume that we have to regulate against functions x 0 (ξ) 0 of a given bandwidth. In a first step (a) we therefore compute a reduced state-space model G red (s) of G(s), in which the resolution of the state x r (ξ) reflects the resolution of the potential x 0 (ξ) accurately. In the present study we regulate against initial conditions with resolution comparable to that considered in [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF], which leads us to a finite-difference discretization of (8) with 50 spatial steps or states, complemented by a 3rd-order Padé approximation of the input delay adding 3 more states. This coarse grid model G red is embedded into a plant P red expressing control requirements in terms of damping and settling time in responses to initial conditions. Here this consists in optimizing the root mean-square energy value of the output signal z r = (z 1 , z 2 ) in response to the white noise disturbance d on the state x r :

P red (s) :          ẋr = A r x r + d + B 2 u r z r = C 1 x r + D 12 u r y r = C 2 x r + D 21 d
where x r ∈ R 54 is the reduced state, d is the exogenous input, understood to represent the impulse caused by the non-zero initial value, and z r = (W x x r , W u u) is regulated similar to what is used in LQG-control, with filters W x = I and W u (s) = s 1+s/a , a = 100, the latter adding another state to x r . The reduced output y r ∈ R 5 of G red represents the 5 distributed measurements (see Fig. 5) in the coarse finite-difference discretization. Our testing shows that five equidistant measurements along [0, L] are sufficient to achieve well behaved responses to initial conditions.

B r (sI

-A r ) -1 C r d K(x) - W u z 1 z 2 y r u r 0 Figure 2. Model matching. Preliminary stabilization via structured H 2 -synthesis using P red with channel d → z r = (z 1 , z 2 ) = (W u u r , x r ).
We then solve the structured H 2 -optimization problem for the channel d → z r of P red in Fig. 2, using systune [START_REF]Robust Control Toolbox 5.0[END_REF] based on [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF][START_REF] Apkarian | Nonsmooth optimization for multidisk H ∞ synthesis[END_REF]. Trial and error reveals the interesting fact that fairly low-order controllers K ∈ K 2 of order 2 are adequate, leading to x ∈ R 21 for K: minimize T d→z r (P red , K) 2 subject to K(x) stabilizes G red internally K(x) ∈ K 2 [START_REF] Noll | Bundle method for non-convex minimization with inexact subgradients and function values, Computational and Analytical Mathematics[END_REF] where the closed-loop transfer T d→z r (P red , K) is the lower linear fractional transformation F l (P red , K). The solution of ( 9) is K 0 in Figure 3 right. The difference with traditional H 2 -control is that K does not have observer structure, but the structure K ∈ K 2 we imposed. The resulting controller K 0 ∈ K 2 is obtained as A simulation for initial condition x 0 (ξ) = ξ(Lξ) is shown in Fig. 5 (left). While this produces the expected good results for G red , we now have to check whether K 0 also stabilizes the infinite-dimensional system G(s). It turns out that this is the case, as the Nyquist test reveals, so that we now proceed to the second part (b) of the model matching method, where the controller is further optimized with regard to the full model. As a result of step (a) of the model matching procedure we have so far obtained a reference model, feedback(G red , K 0 ) and a controller K 0 of the desired structure which stabilizes G red . Application of the Nyquist test, in tandem with 8 boundedness of the closed loop transfer function on jR, show that K 0 also stabilizes the infinite dimensional system G(s). Fig. 4 (left) shows one clockwise encirclement, computed by the ray-crossing algorithm, which due to the known single unstable pole in G, and absence of unstable poles in K 0 , confirms the absence of unstable closed-loop poles. Taking into account that G, K 0 are proper shows that the closed loop transfer function is bounded on jR, hence the loop is H ∞ -stable, and by Proposition 1, is exponentially stable.

K 11 0 = 0.
Remark 5. The fact that K 0 stabilizes not only G red , but even G, could be called accidental. However, recall that within most structures K ∈ K , practical methods leading to a stable closed loop are necessarily heuristics, so remain equally accidental. What can be said in favor of our method to obtain K 0 is that it is the result of a local optimization procedure. [START_REF] Apkarian | Non-smooth optimization for robust control of infinite-dimensional systems[END_REF] starting at initial guess K = K 0 leads to overall optimal K = K * and is a special case of (3). In the reference model on the right K 0 obtained in step (a) remains now fixed.

G K K 0 G red r y y r y y r u u r z - + - + - + Figure 3. Model matching step (b). H ∞ -optimization of channel r → y -y r in
In a second step (b) of the model matching procedure, corresponding to step 3 of algorithm 1, the preliminary stabilizing controller K 0 is refined through an H ∞ model matching problem shown in Fig. 3, which takes the true infinite-dimensional dynamics in G(s) accurately into account. In this step, we have to solve an infinite-dimensional structured H ∞ -control problem covered by the general form (3),

minimize (I + G red K 0 ) -1 G red K 0 -(I + GK) -1 GK ∞ subject to K stabilizes G internally K ∈ K 2 (11)
for which we use the bundle algorithm of [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF][START_REF] Noll | Bundle method for non-convex minimization with inexact subgradients and function values, Computational and Analytical Mathematics[END_REF][START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF], initialized at K = K 0 . As a result of optimization the H ∞ -norm of the mismatch channel r → z = yy r is reduced from 1.81 at K 0 to 0.84 at K * . The optimized controller K * obtained is

K * 11 = 0.

Mixed sensitivity approach

A classical approach to improve system performance is via minimization of the sensitivity function S := (I + GK) -1 . Using appropriate weighting filters W e (s), one can achieve better transient and steady-state responses for all references signals of finite energy. Introducing a penalization of the control effort as before leads to a mixed-sensitivity design problem:

minimize W e (s)(I + G(s)K(s)) -1 W u (s)K(s)(I + G(s)K(s)) -1 ∞ subject to K stabilizes G internally K ∈ K 2 (13)
which is a particular case of (3).

Actuation takes place at the edge ξ = L, and its effect propagates with a unit delay e -s along the spatial dimension. The settling time is essentially determined by the slow dynamics that correspond to the states farthest away from the actuation point. On the other end, states closer to ξ = L settle faster, but in turn are affected by much more turbulent transients, as for instance seen in Fig. 6. This suggests shaping response surfaces using weightings which take the distance to ξ = L into account. This leads us to W e = diag(c 1 , c 2 , c 3 , c 4 , c 5 ) with individually adapted c i , where for simplicity static weightings are sought. Note that c 1 corresponds to the edge ξ = 0, while c 5 is associated with ξ = L. As before, penalization of the control effort uses a high-pass filter W u (s) = (s/10)/(1 + s/1e3).

Fig. 6 (left) shows the simulation of G in closed loop with a first optimal controller K * ∈ K 2 given in ( 14) obtained via nonsmooth optimization (13) started at K 0 from (10) and using W (1) e = diag(3, 0, 0, 0, 0). States close to the edge 10 Increasing the cost at ξ = L via W (2) e = diag(1, 0, 0, 0, 0.2), and starting optimization at K * , now leads to the optimal controller K * * given in [START_REF] Hertlein | An inexact bundle algorithm for nonconvex nondifferentiable functions in Hilbert space[END_REF], which removes this undesirable effect. Simulation of the loop with K * * including the control signal is shown in Fig. 6 (right), the Nyquist plot in 6 (right). The optimal controller is obtained as Finally, note that it is possible to obtain even faster responses by accepting more aggressive control signals and therefore obtaining a more academic than practical solution. By proposition 1 all controllers obtained are exponentially stabilizing.

K * * 11 = 0.
Altogether this study shows that by way of program (3) it is possible to conveniently control the reaction-diffusion equation ( 8) with a single input with delay and 5 distributed measurements by synthesizing a finite-dimensional loworder controller such that the result matches the result obtained in [START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF] using full state feedback.

Control of an anti-stable wave equation

In this second study, we discuss the following boundary control system

x tt (ξ, t) = x ξξ (ξ, t), t ≥ 0, ξ ∈ [0, 1] x ξ (0, t) = -qx t (0, t) (16)
x ξ (1, t) = u(t),
where q > 0, q 1. The state of the system is x(•, t), x t (•, t), the control applied at the boundary ξ = 1 is u(t), and we assume that the measured outputs are

y 1 (t) = x(0, t), y 2 (t) = x(1, t) and y 3 (t) = x t (1, t).
The system has been discussed previously in [START_REF] Fridman | Introduction to Time-Delay Systems, Systems and control, foundations and applications[END_REF], and [START_REF] Smyshlyaev | Boundary control of an anti-stable wave equation with anti-damping on the uncontrolled boundary[END_REF][START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF], where potential applications are mentioned. Its wellposedness can be seen from the functional analytic set-up in [START_REF] Smyshlyaev | Boundary control of an anti-stable wave equation with anti-damping on the uncontrolled boundary[END_REF][START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF], and from the general approach to well-posedness of 1D hyperbolic systems in [START_REF] Zwart | Well-posedness and regularity of hyperbolic boundary control systems on a one-dimensional spatial domain[END_REF].

The transfer function of ( 16) is obtained from the elliptic boundary value problems

s 2 x(ξ, s) = x ξξ (ξ, s), s ∈ C, ξ ∈ [0, 1] x ξ (0, s) = -qsx(0, s) (17)
x ξ (1, s) = u(s), which in this particular situation can be solved analytically:

G(ξ, s) = x(ξ, s) u(s) = 1 s • (1 -q)e sξ + (1 + q)e -sξ (1 -q)e s -(1 + q)e -s .
From this general formula the transfer function of ( 17) is obtained as

G(s) = [G(0, s); G(1, s); sG(1, s)] =: [G 1 (s); G 2 (s); G 3 (s)].
The main challenge in the hyperbolic system ( 17) is that along with the unstable pole at s = 0 it exhibits an infinite number of unstable poles on a line Re(s) = σ > 0. This means that the Nyquist test is not directly applicable.

Preliminary stabilization

Following our scheme in algorithm 1, the first step is to provide a preliminary stabilizing controller K 0 = K(x 0 ) of a simple pre-defined structure. We have to stabilize the system

G(s) =                   2e -s /(1-q) s(1-Qe -2s ) 1+Qe -2s s(1-Qe -2s ) 1+Qe -2s 1-Qe -2s                   =              G 1 (s) G 2 (s) G 3 (s)              where Q = (1 + q)/(1 -q).
A first question is whether G can be stabilized by a finite-dimensional controller. Ignoring the input y 2 , which for stabilization is not required, we choose the structure K(x) = [n 1 (s)/d(s), 0, n 3 (s)/d(s)], with x gathering the unknown coefficients of the polynomials n i (s), d(s) with deg(n i ) ≤ deg(d). Stability of the closed loop T (s) = G(s)/(1 + G 1 (s)K 1 (s) + G 3 (s)K 3 (s) leads to testing whether the quasi-polynomial

(1 -q)s(d(s) + n 3 (s)) + (1 -q)sQe -2s (n 3 (s) -d(s)) + 2n 1 (s)e -s
arising in the denominator of T (s) is stable, i.e., has its roots in C -. While there exist general methods to check stability of quasi-polynomials, cf. [START_REF] Pontryagin | On the zeros of some elementary transcendental functions[END_REF], an ad hoc solution is here to choose n 3 = d, whence the quasi-polynomial simplifies to sd(s) + c(s)e -s , where c(s) = n 1 (s)/(1q) and deg

(n 1 ) = deg(c) ≤ deg(d).
If we choose d(s) = s + x 1 and c(s) = x 2 s + x 3 , then stability of the loop is equivalent to stability of the quasi-polynomial

P(s) = A(s) + B(s)e -s , A(s) = s 2 + x 1 s, B(s) = x 2 s + x 3 ,
which is covered by the discussion in [START_REF] Malakhovski | On stability of second-order quasi-polynomials with a single delay[END_REF]. In their terminology we have a 0 = 0,

a 1 = x 1 , b 2 = 0, b 1 = 1, b 0 = 1.
We are then necessarily in the case m = 1, µ 0 = 0 of [START_REF] Malakhovski | On stability of second-order quasi-polynomials with a single delay[END_REF], so the quasi-polynomial P(s) can only be stable if

x 1 > -1.
Moreover the family P h (s) = A(s) + B(s)e -hs is stable for all 0 ≤ h < h σ,0 , where h σ,0 > 0 is determined as follows.

Let ω σ be the positive real solution of

4ω 3 σ -2ω σ (1 -x 2 1 ) = 1 2 5 -2x 2 1 + x 4 1
and let h σ,0 be the smallest positive solution h of

ω σ h = arg - B( jω σ ) A( jω σ ) + 2kπ, k ∈ N,
where arg(•) ∈ [0, 2π). If we let x 1 = 1 > -1, then 4ω 3 σ = 1, ω σ = 4 -3 , and we get

h σ,0 = 4 3 arg - x 2 j4 -3 + x 3 4 -6 + 4 -3 ,
and since our delay is h = 1, this must now be solved for x 2 , x 3 so that h σ,0 > 1. For instance x 2 = -1 and x 3 = -4 -3 gives argument π/4 in the formula, so that h σ,0 = 16π > 1. The leads to the finite-dimensional stabilizing controller K 0 for G:

K 0 = (1 -q)(s + 4 -3 ) s + 1 0 1 . (18)
A second way to seek preliminary stabilization of ( 16) is to stick to the form K = n(s) d(s) 0 1 , but allow n(s), d(s) to be quasi-polynomials, trying to simplify the denominator quasi-polynomial P as much as possible. A very straightforward way is to let d(s) = a(s) + e -s b(s) with a(s), b(s) polynomials, then the denominator quasi-polynomial simplifies to e -s 2(1-q) (sa(s) + (sb(s) + c(s))e -s ), where c(s) = n 1 (s)/(1q). If we now let c(s) = -sb(s), then K will be stabilizing in the H ∞ -sense as soon as a(s) is stable, because the factor s cancels with the factor s in the numerator. If we choose a(s) = s + c 0 , b(s) = -c 0 for some constant c 0 > 0, then we obtain the controller

K = c 0 (1 -q)s s + c 0 (1 -e -s ) 0 1 , (19)
which in [START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF] was obtained using the back-stepping technique. Since only input and output delays along with realrational terms arise, such controllers are implementable, so we are still in line with our general purpose of computing practically useful controllers.

Performance optimization

Let us now discuss a more systematic way which not only leads to preliminary stabilizing K 0 , but also allows performance optimization. In order to compare with [START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF], we optimize again against the effect of non-zero initial values, using the output y ∈ R 3 , and aiming as before at a convenient implementable controller structure.

We start by putting the system G in feedback with the controller K 0 = [0 0 1], which leads to G = G/(1 + G 3 ), where

G(s) =                   2e -s /(1-q) s(1-Qe -2s ) 1+Qe -2s s(1-Qe -2s ) 1+Qe -2s 1-Qe -2s                   , G(s) =                1 s(1-q) 1+Q 2s 1 2                +                 -1-e -s s(1-q) -Q(1-e -2s ) 2s Q 2 e -2s                 .
Re-write this as G = G + Φ, where G is now real-rational and still unstable, while Φ gathers the infinite dimensional part, but is stable. Then we use that stability of the closed loop ( G + Φ, K) is equivalent to stability of the loop ( G, feedback(K, Φ)), as explained in Fig. 8 See also [START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple i/o delays[END_REF].

Then we construct a finite-dimensional structured controller K = K(x) which stabilizes G. The controller K is then recovered from K through the equation K = feedback(K, Φ), which when inverted gives K = feedback( K, -Φ).

The overall controller is then K * = K 0 + K. Construction of K uses systune where we use pole placement via TuningGoal.Poles imposing that closedloop poles have a minimum decay of 0.9, minimum damping of 0.9, and a maximum frequency of 4.0. The controller structure is chosen as static, so that x ∈ R 3 . A simulation with K * is shown in Fig. 9 (right) and some acceleration over the simulation for backstepping controller (left) from the same nonzero initial value is observed. We recall that the controllers obtained in this section stabilize the system exponentially, as follows from the scheme on the right of Fig. 8, where G, K are finite-dimensional and Φ is stable. This shows that G is exponentially stabilizable and detectable, so that every H ∞ -stabilizing controller is also exponentially stabilizing. In particular, this applies retro-actively also to the controllers ( 18) and [START_REF] Smyshlyaev | Boundary control of an anti-stable wave equation with anti-damping on the uncontrolled boundary[END_REF].

In more detail, we have to relate the trajectories of the system on the right of Fig. 8 to the trajectories on the left. Here we can follow [START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple i/o delays[END_REF], using their formulas ( 13), [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF], to the extent that exponential decay of trajectories on the right of Fig. 8 as assured by the finite-dimensional stabilization achieved via systune leads to exponential decay of trajectories on the left of Fig. 8, which uses of course exponential stability of the infinite-dimensional part Φ. While in [START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple i/o delays[END_REF] work with strictly proper G, K, it suffices for our present argument to suppose that all loops are well-posed. This is for instance guaranteed for a proper K, K, since Φ is strictly proper.

Performance with finite-dimensional control

In this section, we show that the anti-stable wave equation ( 16) may be regulated satisfactorily with a simple 3rd-order finite-dimensional controller. We initialize our procedure with the controller K 0 in [START_REF] Apkarian | Nonsmooth optimization for multidisk H ∞ synthesis[END_REF] obtained via the quasipolynomial test. Then we write the desired structure K(x) as K

(x) = K 0 + K 1 (x), where K 1 (x) = [n 1 /d n 2 /d n 3 /d] and deg(n i ) ≤ 2, deg(d) = 2,
which requires 11 variables. This is a subclass of the class of 3rd-order controllers.

According to section 3.1, we consider the pre-stabilized system G 0 = G(I + K 0 G) -1 , build the closed loop feedback(G 0 , K 1 (x)), and find an initial x 0 ∈ R 11 such that K 1 (x 0 ) is stable and K 1 (x 0 ) ∞ < 1/ G 0 ∞ , so that by the Small Gain Theorem the loop T (G 0 , K 1 (x 0 )) is stable. This is achieved e.g. by K 1 (x 0 ) = n 0 /d 0 [1 1 1] with n 0 (s) = 0.3218s + 0.0643, d 0 (s) = s 2 + 100.1s + 10. Since G 0 has zero unstable poles, and since K 1 (x) is not allowed unstable poles, the Nyquist curve 1 + G 0 K 1 (x) turns now zero times around the origin, and this is maintained during optimization.

We now use the mixed sensitivity approach of section 6.2 again, but under the form minimize

W e (I + G 0 K 1 (x)) -1 W u K 1 (x)(I + G 0 K 1 (x)) -1 ∞ subject to K 1 (x) stabilizes G 0 x ∈ R 11 (20)
where we still have to choose the filters. The final H ∞ -norm in (20) was 1.99, with approximately 1000 frequencies for both stability and performance. In [START_REF] Zwart | Well-posedness and regularity of hyperbolic boundary control systems on a one-dimensional spatial domain[END_REF], we have also constrained the controller to have a minimum decay rate of 1e -3 and minimum damping of 0.1 to keep control on the frequency inter-sample behavior [START_REF] Apkarian | Structured H ∞ -control of infinite dimensional systems[END_REF]. Furthermore the constraint |(1 + K 1 (x)G 0 ) -1 | ≤ 1/0.5 stands for a disk margin of 0.5 hence prohibiting any change in the winding number. Simulations are shown in Fig. 10. Top left shows simulation with K 0 + K 1 (x 0 ), bottom left shows the optimized controller K 0 + K 1 (x ), achieving faster convergence and a much smaller smaller steady-state error beyond 4 sec. Simulations of the slices ξ = 0, ξ = 1 and the control signal are displayed in Fig. 11 from top to bottom and confirm the previous analysis.

Gain-scheduling control

Our last study concerns the case where the parameter q ≥ 0 is uncertain or allowed to vary in time with sufficiently slow variations as discussed in [START_REF] Shamma | Analysis of gain scheduled control for nonlinear plants[END_REF]. We assume that a nominal value q 0 > 0 and an uncertain interval [q, q] with q 0 ∈ (q, q) are given. The authors of [START_REF] Bresch-Pietri | Output-feedback adaptive control of a wave pde with boundary anti-damping[END_REF] schedule their controller (19) using an adaptive control scheme, where the scheduling function uses a nonlinear dynamic estimate q(t) ∈ [q, q] of the anti-damping parameter.

Based on the approach in section 7.2 the following scheduling scenarios are possible. (a) Computing a nominal controller K at q 0 as before, and scheduling through Φ(q), which depends explicitly on q, so that K (1) (q) = feedback( K, -Φ(q)). (b) Computing a K(q) which depends already on q, and using K (2) (q) = feedback( K(q), -Φ(q)). (c) Computing a robust controller K rob for the entire interval.

While (a) is directly based on (3) in its finite-dimensional version based on [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF][START_REF] Apkarian | Nonsmooth optimization for multidisk H ∞ synthesis[END_REF], see also [START_REF] Apkarian | Optimization-based control design techniques and tools[END_REF], as available in systune, leading to K (1) (q), we show that one can also apply our approach to case (b). We use the reduction of section 7.2, see Fig. 8, to work in the finite-dimensional system ( G(q), K(q)), where we now have in addition dependency on q, addressed by a parameter-varying design.

For that we have to decide on a parametric form of the controller K(q), which we chose here as K(q, x) = K(q 0 ) + (qq 0 ) K 1 (x) + (qq 0 ) 2 K 2 (x), and where we adopted the simple static form K 1 (x) = [x 1 x 2 x 3 ], K 2 = [x 4 x 5 x 6 ], featuring a total of 6 tunable parameters. The nominal K(q 0 ) is obtained via the synthesis technique in section 7.2. For q 0 = 3 this leads to K(q 0 ) = [-1.049 -1.049 -0.05402], obtained via systune as in section 7.2. With the parametric form K(q, x) fixed, we now use again the feedback system ( G(q), K(q)) in Fig. 8 and design a parametric robust controller using the method of [START_REF] Apkarian | Parametric robust structured control design, Automatic Control[END_REF], which is implemented in the systune package and used by default if an uncertain closed-loop is entered. The tuning goals are chosen as constraints on closed-loo poles including minimum decay of 0.7, minimum damping of 0.9, with maximum frequency 2. The controller obtained is (with q 0 = 3) K(q, x * ) = K(q 0 ) + (qq 0 ) K 1 (x * ) + (qq 0 ) 2 K 2 (x * ), with numerical values K 1 = [-0.1102, -0.1102, -0.1053], K 2 = [0.03901, 0.03901, 0.02855], and we retrieve the final parameter varying controller for the system G(q) as K (2) (q) = K 0 + feedback( K(q, x * ), -Φ(q)).

The methods are compared in simulation in Figs. 12, 13, 14. Comparison of the simulations in Figs. 12, 13, and 14 indicates that the last controller K 3 (q) achieves the best performance for frozen-in-time values q ∈ [2, 4].

Figure 12. Synthesis at nominal q 0 = 3. Simulations of nominal K = K 0 + feedback( K, Φ(3)) for q = 2, 3, 4. Nominal controller is robustly stable over [q, q].

. Method 1. K obtained for nominal q = 3, but scheduled K(q) = K 0 + feedback( K, Φ(q)). Simulations for q = 2 left, q = 3 middle, q = 4 right Figure 14. Method 2. K(q) = K nom + (q -3) K 1 + (q -3) 2 K 2 and K(q) = K 0 + feedback( K(q), Φ(q). Simulations for q = 2, 3, 4

In conclusion, the study of the hyperbolic system [START_REF] Sano | H ∞ -control of a parallel-flow heat exchange process[END_REF] shows that optimization based on the infinite-dimensional program (3) is required to synthesize finite-dimensional controllers for ( 16), while its finite-dimensional counterpart based on [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF] and implemented in systune is sufficient to synthesize infinite-dimensional controllers of the structure covered by Fig. 8. The major difference with parabolic systems or first-order hyperbolic systems (see e.g. [START_REF] Apkarian | Nonsmooth bundle trust-region algorithm with applications to robust stability[END_REF]) is that preliminary structured stabilization, based on a suitable heuristic, cannot be verified using the Nyquist test. A very first stabilizing controller has to be found by way of some other means, but once this is achieved, the Nyquist test can be brought back to serve to control stability of the loop during optimization.

Figure 1 .

 1 Figure 1. Winding number of Nyquist curve (black) and polygon (blue) agree if change of argument of segment [ f ( jω ν ), f ( jω ν+1 )] equals change of argument of corresponding piece of Nyquist curve f ([ω ν , ω ν+1 ]), as guaranteed by condition (5). Winding number of polygon P f is computed by counting ray crossings in red.

Figure 4 .

 4 Figure 4. Model matching. Nyquist curve 1 + GK 0 (left) with initial controller (10) and 1 + GK * (right). Since F(s) has one unstable pole, one counterclockwise encirclement confirms absence of unstable poles in the loop.

Figure 5 .

 5 Figure 5. Model matching. Simulation of optimized reduced-order system (left) with K 0 based on (9) and infinite-dimensional system (right) with H ∞ -optimal controller K * based on (12). The intermediate result on the left leads to the final result on the right. The stars '*' indicate sensor positions.

Figure 6 .

 6 Figure 6. Mixed sensitivity. Simulation with K * ∈ K 2 obtained with weight W(1) e (left). States near ξ = L show undesirable transient oscillations. Simulation with K * * obtained with W(2) e (right) proves satisfactory.

Figure 7 .

 7 Figure 7. Mixed sensitivity. Left. Simulations with controller K * * : top x(L, t), middle x(0, t) and bottom control signal u(t). Right. Nyquist curve of final K * * .

Figure 8 .

 8 Figure 8. Stability of the closed-loop ( G + Φ, K) is equivalent to stability of the closed-loop ( G, feedback(K, Φ)). See also[START_REF] Moelja | Parametrization of stabilizing controllers for systems with multiple i/o delays[END_REF].

Figure 9 .

 9 Figure 9. Wave equation. Simulations with nonzero initial condition for K obtained by backstepping control (left) and K * = K 0 + K obtained by optimizing feedback( G, K) via systune (right). Both controllers are infinite-dimensional, but implementable.

Figure 10 .

 10 Figure 10. Wave equation. Finite-dimensional controllers obtained by mixed sensitivity in (20). Top left initial K 0 . Bottom left optimal K = K 0 + K 1 (x ). Top right, Nyquist plot 1 + (K(x 0 ) -K 0 )G 0 does not encircle origin, bottom right Nyquist plot 1 + (K(x ) -K 0 )G 0 .

Figure 11 .

 11 Figure 11. Wave equation. Simulations of slices ξ = 0, ξ = 1 and control signal from top to bottom. K : solid blue, K 0 : dotted black, K 0 + K 1 (x 0 ): dashed black. Simulations of K 0 and K 0 + K 1 (x 0 ) are nearly indistinguishable due to small gain restriction on K 1 (x 0 ).

  1343s 2 + 0.4535s + 11.34 s 2 + 10.66s + 38.39

	K * 12 =	0.52s 2 + 1.755s + 45.23 s 2 + 10.66s + 38.39
	K * 13 =	0.7443s 2 + 2.621s + 65.23 s 2 + 10.66s + 38.39	(12)
	K * 14 =	0.5976s 2 + 2.036s + 52.82 s 2 + 10.66s + 38.39
	K * 15 =	0.3446s 2 + 2.621s + 20.47 s 2 + 10.66s + 38.39

its simulation is shown in Figs.

5 (right)

and the corresponding Nyquist plot in Fig.

4

(right). The results resemble those obtained by studies based on full-state information, see

[START_REF] Preur | Feedback stabilization of a 1D linear reaction-diffusion equation with delay boundary control[END_REF]

.

  The 3 × 3 filter W e is chosen diagonal where the first entry is a typical low-pass, which corresponds to the output e 1 . The transfer G 2 is non-minimum phase with unstable zeros at the positionslog(1/2)/2 + jkπ, k ∈ Z, which makes the choice of the second filter diagonal element challenging. The above choice turns out to be a good solution, as it forces tight control in the high-frequency range beyond the first unstable zero atlog(1/2)/2 = 0.346. As third weight we choose a simple static gain 0.01. The static filter W u = 0.01 serves to avoid unrealistic control signals.

	W e (s) = diag	0.01s + 0.5002 s + 0.01429	.99s + 0.0007147 s + 0.07941	0.01 ,
	The final controller based on (20) obtained is
	K 1 =	-2.992s 3 -303.5s 2 -104.7s -0.488 s 3 + 102.2s 2 + 101.7s + 0.522
	K 2 =	-0.04494s 2 -4.047s + 0.001097 s 2 + 101.2s + 0.522	(21)
	K 3 =	1.207s 2 + 122.7s + 0.5271 s 2 + 101.2s + 0.522