Huimeng Zhou

Dan Xu

Xiaoyun Shao

Xizhan Ning

Wang Tao

A Robust Linear-Quadratic-Gaussian Controller for the Real-Time Hybrid Simulation on a Benchmark Problem

Keywords: Real-time hybrid simulation, RLQG, LTR, Robustness, Uncertainties

During a real-time hybrid simulation (RTHS), inevitable time delay of actuators when responding to a command will reduce the accuracy of test results and sometimes even cause unstable testing. The inner-loop controller of an actuator is generally capable of eliminating the effects due to small time-delays. However, if a test specimen behaves nonlinearly, accuracy of RTHS results will be impaired. In addition to the uncertainty of test specimens and transfer system, measurement noises of the displacement and force sensors also require a robust external controller for RTHS. In this paper, a robust linear-quadratic-gaussian (LQG) controller with a Loop Transfer Recovery (LTR) procedure and a polynomial-based feedforward prediction (FP) algorithm is proposed to compensate the adverse effects due to time delay and uncertainties within the RTHS testing system. The stability and robustness of the proposed controller are analysed in the frequency domain using the Nyquist curve and the Bode diagrams. Numerical simulations are then carried out on the benchmark problem using both the proposed robust and the conventional LQG controllers and their performance is compared using the nine evaluation criteria. It is demonstrated that the robust LQG (RLQG) controller outperforms the conventional LQG controller in terms of compensating the parameter uncertainties in the testing system and achieving accurate RTHS results.

INTRODUCTION

A structural system is divided into two parts in a real-time hybrid simulation (RTHS): the critical part is taken as an experimental substructure and the remainder of the structural system is simulated using a numerical model in a computer. Utilizing this substructuring method, the RTHS can physically test large-to full-scale specimens in a laboratory, and with relative accuracy it can emulate the dynamic response of a structural system [1-2].

High-quality control is required in RTHS to synchronize the variables at the substructuring interfaces between the experimental and the numerical substructures and to compensate for the additional actuator dynamics, nonlinearities, uncertainties and time-varying parameters in the experimental substructures. A dynamically substructured system (DSS) scheme was developed to obtain whole structural responses through substructural testing, which is the same principle as for RTHS. However, the DDS scheme approaches this principle from a different substructuring framework by synthesizing a synchronizing automatic controller based upon the conventional transfer function method. Various methods were continuously being developed based on the DDS scheme and applied to different types of RTHS experiments, through which the effectiveness and the robustness of the DDS scheme and its variations were verified [START_REF] Wagg | Substructuring of dynamical systems via the adaptive minimal control synthesis algorithm[END_REF][START_REF] Lim | An adaptive control strategy for dynamic substructuring tests[END_REF][START_REF] Bonnet | Real-time hybrid experiments with Newmark integration, MCSmd outer-loop control and multi-tasking strategies[END_REF][START_REF] Neild | Control issues relating to real-time substructuring experiments using a shaking table[END_REF][START_REF] Tu | Testing of dynamically substructured, base-isolated systems using adaptive control techniques[END_REF][START_REF] Stoten | A comparative study and unification of two methods for controlling dynamically substructured systems[END_REF][START_REF] Stoten | Adaptive control of dynamically substructured systems: the single-input single-output case[END_REF][START_REF] Yamaguchi | Synthesised H∞/μ Control Design for Dynamically Substructured Systems[END_REF].

Alternatively, robust RTHS controllers can be achieved when time delay in the actuator is properly compensated using the inverse-dynamics concept. An accurate dynamic model of the actuator is the key to the success when using these compensators, because the inverse of this model forms the basis of the compensator design [START_REF] Carrion | Model-based strategies for real-time hybrid testing[END_REF][START_REF] Jung | Performance of a real-time pseudo dynamic test system considering nonlinear structural response[END_REF][START_REF] Christenson | Large-Scale Experimental Verification of Semiactive Control through Real-time Hybrid Simulation[END_REF]. Among the different types of actuator models, the polynomial model is perhaps the most popular one. Horiuchi et al. [START_REF] Horiuchi | Real-time hybrid experimental system with actuator delay compensation and its application to a piping system with energy absorber[END_REF] firstly adopted a polynomial model and proposed a compensator based on a single-step-forward extrapolation. Since then, several compensators based on polynomial actuator models [START_REF] Darby | Stability and delay compensation for real-time substructure testing[END_REF][START_REF] Wallace | An adaptive polynomial based forward prediction algorithm for multi-actuator real-time dynamic substructuring[END_REF][START_REF] Tu | Modelling and control issues of dynamically substructured systems: adaptive forward prediction taken as a example[END_REF][START_REF] Ahmadizadeh | Compensation of Actuator Delay and Dynamics for Real-Time Hybrid Structural Simulation[END_REF][START_REF] Chen | Dual compensation strategy for real-time hybrid testing[END_REF] were proposed and adaptive features [START_REF] Wallace | An adaptive polynomial based forward prediction algorithm for multi-actuator real-time dynamic substructuring[END_REF][START_REF] Tu | Modelling and control issues of dynamically substructured systems: adaptive forward prediction taken as a example[END_REF] were added to address the time delay variation during RTHS. Another trend in robust RTHS controller design is to compensate the nonlinearity in the testing system, as the actuator's dynamics (including its time delay) changes with the nonlinear responses of the test specimen. Philips et al. [START_REF] Wu | Actuator dynamics compensation based on upper bound delay for real-time hybrid simulation[END_REF] and Chen et al. [START_REF] Phillips | Model-Based Feedforward-Feedback Actuator Control for Real-time Hybrid Simulation[END_REF][START_REF] Chen | Adaptive model-based tracking control for real-time hybrid simulation[END_REF] proposed an improved model-based servo-hydraulic tracking control method with a feedforward-LQG controller to accurately track the desired displacement in real-time. Chen et al. [START_REF] Liu | A Novel Integrated Compensation Method for Actuator Dynamics in Real-time Hybrid Structural Testing[END_REF][START_REF] Chen | Improving the inverse compensation method for real-time hybrid simulation through a dual compensation scheme[END_REF][START_REF] Chen | Improved Adaptive Inverse Compensation Technique for Real-Time Hybrid Simulation[END_REF] proposed a dual compensation scheme with an outer-loop proportional controller to tune the inner-loop inverse compensation based on the actuator tracking error. Wang et al. [START_REF] Wang | An effective online delay estimation method based on a simplified physical system model for real-time hybrid simulation[END_REF] and Chae et al. [START_REF] Chae | Adaptive time series compensator for delay compensation of servo-hydraulic actuator systems for real-time hybrid simulation[END_REF] proposed to use the first-order and second-order Taylor series to calculate the actuators' commands, respectively, during which the optimum Taylor series coefficients were determined using the least squares method.

The main objective of the aforementioned compensators [START_REF] Carrion | Model-based strategies for real-time hybrid testing[END_REF][START_REF] Jung | Performance of a real-time pseudo dynamic test system considering nonlinear structural response[END_REF][START_REF] Christenson | Large-Scale Experimental Verification of Semiactive Control through Real-time Hybrid Simulation[END_REF][START_REF] Horiuchi | Real-time hybrid experimental system with actuator delay compensation and its application to a piping system with energy absorber[END_REF][START_REF] Darby | Stability and delay compensation for real-time substructure testing[END_REF][START_REF] Wallace | An adaptive polynomial based forward prediction algorithm for multi-actuator real-time dynamic substructuring[END_REF][START_REF] Tu | Modelling and control issues of dynamically substructured systems: adaptive forward prediction taken as a example[END_REF][START_REF] Ahmadizadeh | Compensation of Actuator Delay and Dynamics for Real-Time Hybrid Structural Simulation[END_REF][START_REF] Chen | Dual compensation strategy for real-time hybrid testing[END_REF][START_REF] Wu | Actuator dynamics compensation based on upper bound delay for real-time hybrid simulation[END_REF][START_REF] Phillips | Model-Based Feedforward-Feedback Actuator Control for Real-time Hybrid Simulation[END_REF][START_REF] Chen | Adaptive model-based tracking control for real-time hybrid simulation[END_REF][START_REF] Liu | A Novel Integrated Compensation Method for Actuator Dynamics in Real-time Hybrid Structural Testing[END_REF][START_REF] Chen | Improving the inverse compensation method for real-time hybrid simulation through a dual compensation scheme[END_REF][START_REF] Chen | Improved Adaptive Inverse Compensation Technique for Real-Time Hybrid Simulation[END_REF][START_REF] Wang | An effective online delay estimation method based on a simplified physical system model for real-time hybrid simulation[END_REF][START_REF] Chae | Adaptive time series compensator for delay compensation of servo-hydraulic actuator systems for real-time hybrid simulation[END_REF] is to reduce the tracking errors. However, the requirement for the tolerance to uncertainties, and the acceptable range of the control accuracy of transfer systems (i.e., actuator or shaking table) are not clearly clarified. A major performance trade-off facing by the RTHS controller design is between the tracking error reduction and the increasing error caused by the sensor noises in the measurement. This performance trade-off was considered in some RTHS controller design [START_REF] Gao | Real-time hybrid simulation: from dynamic system, motion control to experimental error[END_REF][START_REF] Ou | Robust integrated actuator control: experimental verification and real-time hybrid-simulation implementation[END_REF][START_REF] Ning | Real-time Hybrid Simulation based on Inner-Loop H∞ Control[END_REF][START_REF] Wu | Sliding Mode Controller for Real-time hybrid Test[END_REF], and one example is the H∞ control method [START_REF] Gao | Real-time hybrid simulation: from dynamic system, motion control to experimental error[END_REF][START_REF] Ou | Robust integrated actuator control: experimental verification and real-time hybrid-simulation implementation[END_REF][START_REF] Ning | Real-time Hybrid Simulation based on Inner-Loop H∞ Control[END_REF]. The linear-quadratic-Gaussian (LQG) controller [START_REF] Wu | Sliding Mode Controller for Real-time hybrid Test[END_REF][START_REF] Yang | New optimal control algorithms for structural control[END_REF] provides a systematic procedure to address the aforementioned performance trade-off. The LQG controller was utilized in [START_REF] Yang | New optimal control algorithms for structural control[END_REF][START_REF] Yang | Application of optimal control theory to civil engineering[END_REF] the vibration control of civil engineering structures under random loadings, during which the trade-off between the reduction in structural responses and the increased requirement of control forces were considered. Wu et al. [START_REF] Wu | Numerical model of the servo-hydraulic loading system for rea l-time substructure testing[END_REF] employed the LQG to control the actuator in an RTHS, its superiority over the PID controller is proved by numerical simulations. Philips et al. [START_REF] Wu | Actuator dynamics compensation based on upper bound delay for real-time hybrid simulation[END_REF] utilized a feedback-LQG controller to improve robustness and alleviate the need of online modification of the feedforward controller under the changing specimen conditions. The robustness of the LQG controller can be further improved by using the Loop Transfer Recovery (LTR) procedure, which was proposed by Doyle et al. [START_REF] Doyle | Multivariable Feedback Design: Concepts for a Classical/Modern Synthesis[END_REF] for the general multivariable feedback controller design in aeronautical engineering. Nevertheless, the application of the LQG controller and the LTR procedure to RTHS experiments have to meet the requirements of real-time processing and robustness at the same time.

In this paper, a robust LQG (RLQG) controller is proposed by combining the LQG and the LTR procedure to improve the robustness and accuracy of RTHS experiments. In addition, a polynomial-based forward prediction (FP) is used to compensate the time delay. The theoretical description of this controller is presented in Section 2. The benchmark problem is briefly presented in Section 3. Then, the robustness of the RLQG controller is analysed by the Nyquist curve and the Bode diagrams in Section 4. Numerical simulations are presented in Section 5 to demonstrate the performance of this controller on the benchmark problem as an example.

DESIGN OF THE RLQG CONTROLLER

A detailed description of the RLQG controller design is provided in this section. A block diagram demonstrating the major components of the RLQG controller is shown in Fig. 1, along with the RTHS testing system including both actuator and specimen highlighted in red.

The RLQG controller originated from the linear-quadratic-gaussian (LQG) controller, which includes a linear-quadratic-optimal-regulator (LQR) controller (shown in green in Fig. 1) and a state observer (shown in blue in Fig. 1). The design of the LQR controller and the state observer is presented in the Subsections 2.1 and 2.2, respectively. In Fig. 1, yn is the output of a numerical substructure. yGc is the command signal from the controller to the actuator-specimen system. xm is the measured displacement of the experimental substructure and xcom is the output of FP compensator. The objective of this RLQG controller is to make the measured displacement xm accurately tracking the output of the numerical substructure yn.

In theory, the state feedback gain matrix [K -kI] is determined in the LQR controller to be applied to the states that are obtained directly from the RTHS testing system. However, in practice, it is difficulty to accurately estimate the states because of the measurement noises and the disturbances within the testing system. A state observer is, therefore, introduced to estimate the state to be fed back. However, the conventional state observer will reduce the robustness of the LQG controller and the LRT procedure is therefore adopted to improve the performance of the state observer, as explained in Subsection 2.3. Subsection 2.4 describes the derivation of the RLQG controller (see the orange frame) in the discrete form. To further reduce the time delay effect in the RTHS system, the FP compensation, described in Subsection 2.5, is integrated that extrapolates the output of numerical substructure yn before it is sent to the LQR controller.

LQR controller

A numerical model of the actuator-specimen system must be identified for the LQR controller design. In general, an actuator attached to a nonlinear specimen remains a nonlinear system with nonlinearities from both the specimen and the actuator (i.e., the square root relationship among the oil pressure, flow, and oil leakage effects [2,[START_REF] Jung | Performance of a real-time pseudo dynamic test system considering nonlinear structural response[END_REF]). In this paper, only the uncertainty of the actuator-specimen system model is considered. A linear actuator-specimen transfer function model [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF] is adopted herein for the convenience of the later analysis. This transfer function is expressed as

0 P 5 4 3 2 5 4 3 2 1 0 ( ) B G s A s A s A s A s A s A = + + + + + (1)
where B0, A0, A1, A2, A3, A4, and A5 are the parameters defined in [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF]. This transfer function is further expressed in the form of a differential equation as (2) By defining

1 2 3 4 (4) 5 m m m m m x x x x x x x x x x                 = =                     X & && &&& (3)
the state-space model of the actuator-specimen system is obtained:

 = +   =   X AX B CX & Gc m y x (4) where 0 5 1 5 2 5 3 5 4 5 0 1 0 0 0 0 0 1 0 0 = 0 0 0 1 0 0 0 0 0 1 / / / / / A A A A A A A A A A               - - - - -  

A

, B=[0 0 0 0 B0/A5] T , and C=[1 0 0 0 0].

Because the actuator-specimen system model has no integrator (type-0 plant) [START_REF] Ogata | Modern Control Engineering[END_REF], which will result in steady state error under the step input. To ensure zero steady state error, an integrator is inserted in the feedforward path between the error comparator and the actuator-specimen to build a type-1 [START_REF] Ogata | Modern Control Engineering[END_REF] servo system, as shown in Fig. 1 (see the green frame). For the LQR controller, the actuator command is determined as

=  = - + = -  = - -  KX KX CX & Gc I com m com y k x x x ξ ξ (5)
where

=[ -] Ι k K K
is the feedback gain matrix of the state X of type-1 servo system, and

[ ] = X X T ξ
consists of the state vector of specimen-actuator system X and ξ, which is a scalar representing the output of the integrator (a state variable of the system). Then, substituting Equation [START_REF] Bonnet | Real-time hybrid experiments with Newmark integration, MCSmd outer-loop control and multi-tasking strategies[END_REF] in Equation ( 4), the state-space model of the equivalent regulation problem system becomes:

[ ]

1  = +   =   X AX 0 CX & T com m x x (6) where = 0 Ι -     -   A BK B A C
k , and

[ ] = 0 C C .
The design of the LQR controller can be divided into two steps: the design of the Quadratic performance function J and the determination of the feedback gain matrix K . The Quadratic performance function is usually defined by a linear algebra equation of the system's state and the control input. Thus, the feedback control law, which is combined with the linear part of the regular form of the state-space equation [START_REF] Ogata | Modern Control Engineering[END_REF], can produce a stable and desirable system performance simply by employing linear system design approaches. The Quadratic performance function is defined as

0 1 [ ( ) ( ) ( ) ( )] 2 ∞ = + ∫ X QX T T Gc Gc t J t t y t Ry t dt (7)
Then the feedback gain matrix K is obtained by solving the Riccati function [START_REF] Ogata | Modern Control Engineering[END_REF] when the performance function is given by Equation [START_REF] Tu | Testing of dynamically substructured, base-isolated systems using adaptive control techniques[END_REF]. In Equation [START_REF] Tu | Testing of dynamically substructured, base-isolated systems using adaptive control techniques[END_REF], the weight scalar R and the matrix Q are defined as

χ δ λ β α γ         = == ⋅             Q ; (8)
where ke and me are the stiffness and mass of the experimental substructure [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF], respectively; α, β, γ, δ, χ, and λ are the parameters to be designed. The control effect is not affected if the ratio of γ and χ is constant. Substituting Eq. ( 8) into Equation [START_REF] Tu | Testing of dynamically substructured, base-isolated systems using adaptive control techniques[END_REF], which represents the sum of the kinetic energy, the optimum feedback gain K is then determined according to the LQR algorithm, which will result in the minimum of Equation [START_REF] Tu | Testing of dynamically substructured, base-isolated systems using adaptive control techniques[END_REF]. The detailed derivation process can be found in reference [START_REF] Ogata | Modern Control Engineering[END_REF]. This optimum feedback gain matrix K also can be computed using the command 'lqr' in the MATLAB toolbox [START_REF]MATLAB, version 8.3 (R2014a[END_REF] with the inputs as shown below:

= ( , , , ) lqr R K A B Q (9)
The open-loop transfer function 0 ( ) g jω of the LQR controller based on Equation ( 6) is shown below which will be utilized later in the stability and robustness analysis:

1 0 ( ) ( ) ( ) ( ) - = = - - K I A B Gc m Y s g j j X s ω ω (10)

State observer

The state of the system usually cannot be directly observed whereas a state observer is used to observe the state of the system. LQG controller is the assemble of the LQR controller and the state observer.

When the observability condition is satisfied, the state observers can be designed. As shown in Fig. 1, the state observer estimates the state variables based on the measured displacement of the experimental substructure xm and the command signal from the controller to the actuator-specimen system yGc. Thus, we define the mathematical model of the observer [START_REF] Ogata | Modern Control Engineering[END_REF] to be

E E E ˆˆˆ= + ( ) ( ) + = + - = - + -= - X AX B K CX A K C X B K CX & & Gc m Gc m com m com y x y x x x x ξ ( 11
)
where KE is the observer gain matrix. It is a weighting matrix to the correction term involving the difference between the measured displacement of the experimental substructure xm and the estimated displacement of the experimental substructure = CX m x , X denotes the estimate of the state vector X, and ˆm

x denotes the estimate of the measured displacement of the experimental substructure xm. This is a full state observer. And the observer error and its derivative are defined as = e X X [START_REF] Jung | Performance of a real-time pseudo dynamic test system considering nonlinear structural response[END_REF] whose characteristic equation [START_REF] Ogata | Modern Control Engineering[END_REF] can be obtained as

E =( - ) e A K C e & (13)
The optimum values in matrix KE will make error vector e asymptotically stable and sufficiently quickly converge to obtain the accurate state estimation. The observer gain matrix KE can also be obtained by using the lqe2 function in the MATLAB Toolbox [START_REF]MATLAB, version 8.3 (R2014a[END_REF] with the inputs as shown below:

E = 2( , , , , ) e e lqe Q R K A B C ( 14
)
where Qe is the covariance matrix of the system noises, and Re is the covariance matrix of the measurement noises.

The LQG controller model may now be expressed as Ĝc I y k ξ = -+ KX [START_REF] Darby | Stability and delay compensation for real-time substructure testing[END_REF] Combining Equations ( 15) and [START_REF] Carrion | Model-based strategies for real-time hybrid testing[END_REF]

leads to: [ ] E 1   = + +  = -   X AX 0 K KX & % % % % % T com m Gc x x y ( 16
)
where

E 0 I k - -   =   -   A K C BK B A C % , [ ] E E 0 T = K K % , ˆT ξ   =   X X %
is the assemble of the estimate of the state vector X and an intermediate variable ξ.

The open-loop transfer function of the state observer based on Equation ( 11) is

1 0 E ˆ( ) ( ) = ( ) ( ) - = - C I A K m m X s g j s X s ω ( 17
)
The open-loop transfer function of the LQG controller based on Equation ( 16) is

1 0 E ( ) ( ) = ( ) ( ) - = - - K I A K % % Gc m Y s g j s X s ω ( 18
)
Although the state observer itself can be designed to have a good robustness to the noises and disturbances, a robust state observer may reduce the robustness of the LQG controller [START_REF] Doyle | Multivariable Feedback Design: Concepts for a Classical/Modern Synthesis[END_REF][START_REF] Ogata | Modern Control Engineering[END_REF]. The LTR procedure is therefore adopted herein to improve the controller's robustness as discussed next.

LTR procedure

The open loop transfer function [START_REF] Ahmadizadeh | Compensation of Actuator Delay and Dynamics for Real-Time Hybrid Structural Simulation[END_REF] designed in Section 2.2 can be recovered to target open-loop transfer function [START_REF] Yamaguchi | Synthesised H∞/μ Control Design for Dynamically Substructured Systems[END_REF] by a modified state observe design procedure [START_REF] Doyle | Multivariable Feedback Design: Concepts for a Classical/Modern Synthesis[END_REF]. The covariance matrix Qe in Equation ( 14) is replaced by Q1, which is defined as Q1=qQe, q is a scalar parameter. When q→∞, the open-loop transfer function of the LQG controller will be close to the open-loop transfer function of the LQR controller as

1 1 E lim ( ) ( ) q s s - - →∞ - = - K I A K K I A B % % (19)
Parameters in KE can be calculated using the lqe2 function from the MATALB Toolbox [START_REF]MATLAB, version 8.3 (R2014a[END_REF], similar to how it was used in Equation ( 14):

E 1 = 2( , , , , ) e lqe Q R K A B C (20)
The open-loop transfer function of the RLQG controller is the same as Equation [START_REF] Ahmadizadeh | Compensation of Actuator Delay and Dynamics for Real-Time Hybrid Structural Simulation[END_REF].

Implementing the RLQG controller

The RLQG controller proposed in this paper is based on the LQG controller and added with the LTR procedure to improve its robustness. The RLQG controller can be rewritten from Equation ( 16)

  = +  =   X AX Bu CX & % % % % % % % Gc y (21)
where

E = 0 1       K 0 B % , = - C K % and [ ] = T m com x x u % .
The discrete time model [START_REF] Åström | Computer-controlled systems theory and design[END_REF] of the LQG controller is written as

* * * ( +1)= ( ) ( ) ( ) ( ) ( ) d d d d Gc d d k k k k y k k + = X A X B u C X % % % % (22)
where A * d, B * d, and C * d are determined using the MATLAB command 'c2d' [START_REF] Åström | Computer-controlled systems theory and design[END_REF], which converts a model from continuous time (defined by A % , B % , and C % ) to discrete time by the Tustin's method [START_REF] Åström | Computer-controlled systems theory and design[END_REF].

( +1) d k X %
is the controller state vector at time t=(k+1)Δt (k is the discrete index, and Δt is the simulation sampling duration), ( )=[ ( )

( ) ] T m com k x k x k u %
, where ( )

com
x k is the output of FP compensator, and ( ) m x k is the measured displacement of the experimental substructure; both are at time t=kΔt.

FP compensation method

The RLQG controller is aimed to cancel the adverse effects to the RTHS testing system due to uncertainties. To compensate the residual time delay, the polynomial-based forward prediction (FP) compensation algorithm [START_REF] Wallace | Real-time dynamic substructuring for mechanical and aerospace applications; control techniques and experimental methods[END_REF] is applied as shown in Fig. 1 (the purple frame), which determines the output of FP compensator before it is sent to the RLQG controller , +1 1 1 ( )

+ +   = ∆ ∆ =   a L N com k n , k P x P t P t y ( 23
)
where P is the forward prediction step, a is the coefficient matrix determined using the least square fitting based on the output vector of numerical substructure of the previous n steps.

1 , 1 0 , 1 1 , 2 1 0 0 1 ( ) 1 ( 1) ( ( 1) ) - + - -+             -∆ -∆       = = =             --∆ --∆         a X y L L M M M M M M L n k N n k n N n k n N y a y a t t y a n t n t (24)
where [•] -1 is the pseudo-inverse of the matrix, the parameter n must be larger than or equal to the order of the polynomial N+1, the upper/top limit of n depends on the parameters N and the maximum signal frequency fmax, whose relationship is discussed in details in the reference [START_REF] Wallace | Real-time dynamic substructuring for mechanical and aerospace applications; control techniques and experimental methods[END_REF]. Increasing n will smooth the command history and decrease the noises' effect. The larger the value of N, the better fitting of the polynomial. However, the FP algorithm tends to be unstable when the value of N is too large. Therefore, N is usually set between 2~4. The parameter P is the prediction step, which is constant during a test and can be determined by dividing the time delay τ by the simulation sampling duration Δt.

BENCHMARK PROBLEM

The actuator-specimen model of the benchmark problem is defined in Equation (1) whose parameters, and their associated uncertainties are described in reference [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF]. Numerical simulations of RTHS (i.e., vRTHS) were conducted in the time domain using the vRTHS tools programmed in Matlab TM , which are part of the benchmark problem provided as well.

A laboratory-size model of a typical three-story steel frame structure is the prototype structure in the benchmark problem. For RTHS, the reference model is partitioned into numerical and experimental substructures. Five partitioning cases are defined in the benchmark problem by varying the parameters of the reference structure, as summarized in Table 1. Based on the predictive stability indicator [START_REF] Maghareh | Predictive stability indicator: A novel approach to configuring a real-time hybrid simulation[END_REF], partitioning cases 1-4 are slightly sensitive to the experiment setup to desynchronization, and partitioning case 5 is moderately sensitive. Case 5 is specifically designed and added here to demonstrate the robustness of the proposed RLQG+FP controller. The response of the entire structural model is numerically simulated and the power spectral of the reference displacement of the first floor x1 is determined and shown in Fig. 2. It can be seen that the main frequency content of this reference displacement is from 0~20 Hz. This frequency range is referred in the bandwidth design of the controller system, as discussed in Sections 4.2 and 4.3.

STABILITY AND ROBUSTNESS ANALYSIS OF THE RLQG CONTROLLER

The stability and robustness analysis of the LQR, LQG controller and the proposed RLQG controllers are presented in this section, whose designs are explained in Subsections 2.1, 2.2, and 2.3, respectively. The parameters used herein for the three controllers are determined based on the parameters of the benchmark problem model, as summarized in

Stability and robustness theorem

When design the LQR, LQG and RLQG controllers, the frequency condition for the optimum tuning system must be satisfied [START_REF] Zhen | Linear system theory (Second Edition)[END_REF]. This frequency domain condition can also be expressed in a geometrical form as: Theorem Ⅰ [START_REF] Zhen | Linear system theory (Second Edition)[END_REF]:

The curve of open-loop frequency response 0 ( ) g jω from ω=0 to ω=∞ in the complex plane must not intersect the unit circle Γ(-1,j0), and the curve g0(jω) has a finite tangent point with the unit circle. The optimum tuning system must satisfy the following two requirements according to theorem Ⅱ: Theorem Ⅱ [START_REF] Zhen | Linear system theory (Second Edition)[END_REF]: (ⅰ) with less than ±60° phase margin; (ⅱ) with (1/2,∞)gain margin, where the phase margin and the gain margin are defined as: Gain margin=1/|g0(jωα)|, when ∠g0(jωα)=180° (25)

Phase margin=θ, when |g0(jωα)|=1

(26)

Stability and robustness analysis of the LQR and LQG controllers

Using the parameters defined in Table 2 of the two controllers, the Nyquist curves and the Bode diagrams of the LQR and the LQG controllers for the benchmark problem are shown in Fig. 3 based on their respective open-loop transfer function defined in Equations ( 10) and [START_REF] Ahmadizadeh | Compensation of Actuator Delay and Dynamics for Real-Time Hybrid Structural Simulation[END_REF]. In addition, the Nyquist curve and the Bode diagram of the state observer are also included in Fig. 3 based on its transfer function defined in Equation ( 17). The Nyquist curve of the LQR controller shown in Fig. 3(a) in red has only two tangent points within the unit circle Γ(-1,j0). This curve defined by g0(jω) satisfies the robustness requirements according to Theorems I, and Ⅱ. From Fig. 3(b), the phase margin Pm is greater than 60°, and the gain margin Gm is infinite. The slope of the LQR controller curve around the crossing frequency is approximately -20dB, which implies good robustness [START_REF] Anderson | Optimal Control, Linear Quadratic Methods[END_REF]. The controller bandwidth is 474 rad/s, which is approximately 4 times of 125.6rad/s (i.e., 20Hz, the largest frequency of input signal). The Nyquist curves of the LQG controller (blue) and the state observer (green) are also shown in Fig. 3(a). Both curves have only one tangent point within the unit circle Γ(-1, j0) at the origin and infinite phase margin, which can also be observed in Fig. 3(b). This infinite phase margin of the LQG controller does not meet the robustness requirement according to the Theorem Ⅱ (ⅱ) [START_REF] Zhen | Linear system theory (Second Edition)[END_REF], although the state observer has infinite gain and phase margin, which meet the robustness requirement based on the same theorem. It can be seen from Fig. 3(b) that the LQG g0(jω) curve decreases more rapidly than the LQR curve and its gain margin Gm is less than that of the LQR controller. This reveals that the LQG controller has poor robustness even though it is designed based on the combination of a robust LQR controller and a robust state observer [START_REF] Anderson | Optimal Control, Linear Quadratic Methods[END_REF]. The RLQG controller is, therefore, proposed for the RTHS benchmark problem to address this conflict in the robustness performance of the LQG controller, whose stability and robustness analysis is presented next.

Stability and robustness analysis of the RLQG controller

The Nyquist curve and the Bode diagram of the RLQG controller are shown in Fig. 4, where the blue and red lines represent the RLQG and LQR controllers, respectively. From Fig. 4(a), it can be seen that the RLQG controller has no tangent point within the unit circle Γ(-1,j0), so it has an infinity gain margin which can also be seen in Fig. 4(b). This infinite gain margin of the RLQG controller after adopting the LTR procedure satisfies the robustness requirement according to Theorem Ⅱ(i). Fig. 4(c) shows that the RLQG g0(jω) curve decreases more rapidly than the LQR curve and the slope of the RLQG curve around the crossing frequency is also approximately -20 dB, the gain margin Gm of the RLQG controller is recovered to infinity as for the LQR controller. The infinite gain margin of the RLQG controller shows that it has the same robustness as the LQR controller. The system bandwidth is 383rad/s, which is approximately 3 times 125.6rad/s (i.e., 20Hz, the largest frequency of the input signal). Although, the RLQG controller has a slightly smaller phase margin than the LQR controller, which implies that it has less tolerance to time delay than the LQR and LQG controllers, a fixed time delay compensation such as the FP algorithm can address this issue. To this end, the RLQG controller proposed in this paper is considered to have adequate robustness for the benchmark problem.

NUMERICAL SIMULATION

Numerical simulations of the benchmark problem using the LQG controller and the proposed RLQG controller combined with the FP compensator were conducted. The simulation time step Δt=1/4096 sec. The same parameters listed in Table 2 for the LQG and RLQG controllers are used in the numerical simulation. The parameters of the FP compensator are the same when applied to the two controllers, which are N=3, n=10. The time delay is 0.01 s, when only the RLQG controller is used, which satisfies the maximum time delay of 10ms required in the benchmark problem (Section 3.4 in [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF]). The prediction step P is therefore set to 0.01/(1/4096)=40.96. The simulation results of the virtual RTHS (vRTHS) using the LQG+FP and the RLQG+FP controllers are presented and compared below.

Displacement response and tracking error time histories

Displacement time histories of the reference model and the RTHS are shown in Fig. 5 for all five partitioning cases. As can be seen from the response comparison in the range of 9.54-9.6 seconds, the RTHS responses using the RLQG+FP controller follow the reference model responses better than that using the LQG+FP controller. The entire time histories of the reposes are not shown in Fig. 5 because similar observations can be made in other segments. The tracking error between 5-10 seconds of the two controllers are shown in Fig. 5 for comparison. The absolute values of the tracking errors of the RLQG+FP controller are less than those of the LQG+FP controller in all 5 cases. And the RLQG errors experience less fluctuation than the LQG errors. The larger fluctuation in the LQG errors explains the larger difference between the responses obtained from the vRTHS using the LQG+FP and the reference model when compared to the RLQG+FP controller. Therefore, both the time history responses and the tracking error comparison demonstrate that the RLQG+FP controller is more robust than the LQG+FP controller for the five cases defined in Table 2.

Nine Evaluation Criteria Values

To consistently compare the controller's performance proposed for the benchmark problem, nine evaluation criteria are defined in [START_REF] Silvaa | Benchmark Control Problem for Real-Time Hybrid Simulation[END_REF]. The values of these nine evaluation criteria were computed based on the vRTHS responses and are listed in Table 3, where white cells show the results of the RLQG+FP controllers and grey cells show the results of the LQG+FP controllers for easy comparison. The detail data of these nine evaluation criteria are listed as A1-A10 in Appendix A.

As can be seen from Table 3, the J1 values are zeros for all 5 cases for the RLQG controller demonstrating its capability to achieve accurate control result even when a fixed FP compensator is used. In addition, J2 and J3 values of all 5 cases are less than 3%, which echo the good match between the reference model responses and the vRTHS responses shown in Fig. 5. The variances (Std. dev.) of the J2 and J3 values are very small revealing the robustness of the RLQG controller when it is applied to five partitioning cases. The variance values of J4-J9 are relatively larger than those of J2 and J3. This maybe be attributed to the uncertainties of the experimental substructure stiffness ke that not only affects the control accuracy, but also affects the restoring force feedback to the numerical substructure, which in turn results in relatively larger errors in the global response of the reference model. Nevertheless, the mean values of these criteria (J4-J9) are generally less than 5%, within the acceptance level for most vRTHS.

Comparing the J1-J9 values of the RLQG+FP and the LQG+FP controllers in Table 3, it can be seem that generally the nominal values of RLQG+FP are smaller than that of the LQG+FP controller, while the variances (Std. dev.) of the J1-J9 values using the RLQG+FP are much smaller than of those of the LQG+FP controller (usually less than 1/10). Both the small nominal values and the variances of the RLQG+FP controller demonstrate its superior robustness when compared to the LQG+FP controller. To facilitate easy visual comparison of the evaluation criteria values among different partitioning cases when using the RLQG+FP controller, Fig. 6 is provided. It can be seen that J4-J9 values of case 1 (red) are smaller than those of case 4 (black), which indicates that the lightly damped structure is more sensitive to the parameter uncertainties compared to highly damped structure. J4-J9 values of case 5 (cyan) are the largest among the five cases because the reference model of this case has the smallest damping, further proves that damping has a great effect on the accuracy of the RTHS results. It is interesting to see that J2-J3 values of case 4 are smaller than those of case 1 and J2-J3 values of case 5 are the smallest. This phenomenon can deduce that for the same tracking error of different displacement command yn value (as is seen from Fig. 6), the lower the damping of the reference model, the smaller the values of J2 and J3. From Fig. 6, J2-J7 values of case 3 (green) are smaller than those of case 4, which indicates that the larger the mass ratio between the experimental substructure and the numerical substructure, the more sensitive of the RTHS results to the parameter uncertainties. These results also agree with the predictive stability indicator (PSI) analysis, in which PSI of case 3 and case 1 are smaller than that of case 4 and PSI of case 5 is the smallest in Table 2. Figure 6. vRTHS evaluation criteria for cases 1, 2, 3, 4, and 5.

CONCLUSION

A robust LQG (RLQG) controller is proposed in this paper to accurately control the RTHS system, and it is implemented in the benchmark problem with parameter uncertainties. The design of the RLQG controller is presented focusing on its integration of several existing control methods and strategies. The stability and robustness of the RLQG controller are analyzed using the Nyquist curve and Bode diagram with gain and phase margins. The stability and robustness analysis show that the RLQG controller has a high level of robustness with an infinite gain margin. Numerical simulations of RTHS are performed on the benchmark problem, considering five partitioning cases. The simulation results demonstrate that the RLQG controller is more robust than the than LQG controller and when applied in RTHS of the benchmark problem, RLQG yields more accurate response histories than the LQG when compared to the responses of the reference model. The sensitivity of the RTHS results to the damping and the mass ratio between the experimental substructure and the numerical substructure is observed, which may help future RTHS partition design to achieve higher accuracy.

XFigure 1

 1 Figure 1 Block diagram of the RLQG controller

Figure 2

 2 Figure 2 Frequency spectrum of reference signal x1.

Figure 3

 3 Figure 3 Nyquist curve and Bode diagram of the LQG controller and plant.

FrequencyFigure 4

 4 Figure 4 Nyquist curve and Bode diagram of the RLQG controller and plant.From Fig.4(a), it can be seen that the RLQG controller has no tangent point within the unit circle Γ(-1,j0), so it has an infinity gain margin which can also be seen in Fig.4(b). This infinite gain margin of the RLQG controller after adopting the LTR procedure satisfies the robustness requirement according to Theorem Ⅱ(i). Fig.4(c) shows that the RLQG g0(jω) curve decreases more rapidly than the LQR curve and the slope of the RLQG curve around the crossing frequency is also approximately -20 dB, the gain margin Gm of the RLQG controller is recovered to infinity as for the LQR controller. The infinite gain margin of the RLQG controller shows that it has the same robustness as the LQR controller. The system

Figure 5 .

 5 vRTHS response comparison of the 1st floor displacement under El Centro excitation, cases 1, 2, 3, 4, and 5.

ACKNOWLEDGEMENT

  The authors are grateful to Professor Wu of Wuhan University of Technology and Dr. Wang of Harbin Institute of Technology for their suggestions with the development of the RLQG controller. The research is financially supported by the Scientific Research Fund of Institute of Engineering Mechanics, China Earthquake Administration (Grant No. 2018D10, 2018D03, 2016B09, 2017A2), National Science Foundation of China (Grant No. 51408565, 51378478, 51878630, 51878525). REFERENCE 1. Nakashima M, Kato H, Takaoka E. Development of real-time pseudodynamic testing.

Table 1 :

 1 RTHS partitioning cases.

	Partitioning	Reference	Reference	Predictive	Sensitivity of
	configuration	floor mass	modal damping	stability indicator	experiment setup to
		(kg)	(%)	[42]	desynchronization
	Case 1	1000	5	1.05~1.1	slightly sensitive
	Case 2	1100	4	0.9~1	slightly sensitive
	Case 3	1300	3	0.9	slightly sensitive
	Case 4	1000	3	0.8~0.9	slightly sensitive
	Case 5	950	2	0.6~0.7	moderately sensitive

Table 2

 2

	below.

Table 3 :

 3 vRTHS evaluation criteria for cases 1, 2, 3, 4, and 5 for El Centro earthquake excitation

		Case (#)	J1 (msec)	J2 (%)	J3 (%)	J4 (%)	J5 (%)	J6 (%)	J7 (%)	J8 (%)	J9 (%)	Method
		Nominal value	0	2.89	2.344 3.409 2.621 1.799 1.806 0.976 0.981	RLQG +
	1	Std. dev. Nominal value	0 0.414	0.026 0.056 0.058 0.043 3.798 3.941 8.787 5.66	0.13 7.225 7.249 3.803 3.799 0.128 0.069 0.069	FP LQG +
		Std. dev.	0.484	0.567 0.771 2.281 1.219 2.238 2.247 1.114 1.115	FP
	2	Nominal value Std. dev. Nominal value	0 0 0.295	2.711 2.256 3.157 2.594 1.629 1.635 0.895 0.888 RLQG+ FP 0.021 0.042 0.036 0.045 0.088 0.089 0.048 0.05 LQG 4.069 4.406 8.789 5.615 7.418 7.436 4.141 4.134 +
		Std. dev.	0.413	0.854 1.107	1.85	1.005 1.625	1.63	0.95	0.95	FP
	3	Nominal	0	1.752 1.596 3.057 2.616 2.501 2.501 1.742 1.741	RLQG

Table A -

 A Table A-4: vRTHS evaluation criteria for cases 4 for El Centro earthquake excitation, using a RLQG controller improved with a FP compensator. 5: vRTHS evaluation criteria for cases 5 for El Centro earthquake excitation, using a RLQG controller improved with a FP compensator.

	a RLQG controller improved with a FP compensator. 14 0 1.75 1.59 3.22 2.73 5 0 1.63 1.84 5.56 4.09	2.71 5.31	2.71 5.31	1.88 2.96	1.88 2.95
	15 6	0 0	1.75 1.65	1.6 1.90	3.02 5.18	2.59 3.92	2.45 4.88	2.45 4.88	1.71 2.73	1.71 2.71
	Number 16 7	J1 0 0	J2 1.75 1.66	J3 1.58 1.91	J4 3.15 5.28	J5 2.73 3.92	J6 2.59 5.05	J7 2.59 5.05	J8 1.8 2.82	J9 1.8 2.81
	(#) 17 8	(msec) 0 0	(%) 1.76 1.64	(%) 1.61 1.87	(%) 2.99 5.55	(%) 2.58 4.08	(%) 2.4 5.31	(%) 2.4 5.31	(%) 1.67 2.97	(%) 1.67 2.95
	1 18 9	0 0 0	2.7 1.75 1.61	2.25 1.6 1.76	3.16 3.14 6.33	2.59 2.68 4.55	1.65 2.62 6.09	1.65 2.62 6.10	0.91 1.82 3.40	0.90 1.82 3.39
	2 19 10	0 0 0	2.71 1.77 1.62	2.28 1.61 1.80	3.17 2.9 6.20	2.64 2.45 4.47	1.67 2.32 5.99	1.68 2.32 5.99	0.92 1.62 3.34	0.91 1.62 3.33
	3 20 11	0 0 0	2.68 1.76 1.66	2.18 1.6 1.91	3.24 3.03 5.13	2.63 2.59 3.84	1.79 2.46 4.86	1.8 2.46 4.86	0.99 1.72 2.71	0.98 1.72 2.70
	4 21 12	0 0 0	2.71 1.75 1.62	2.26 1.59 1.79	3.17 3.1 6.10	2.61 2.62 4.43	1.64 2.58 5.86	1.65 2.58 5.86	0.90 1.8 3.27	0.90 1.79 3.25
	5 Table A-3: vRTHS evaluation criteria for cases 3 for El Centro earthquake excitation, using 0 2.7 2.21 3.2 2.65 1.73 1.73 0.95 0.94 13 0 1.62 1.82 5.78 4.21 5.51 5.52 3.08 3.06
	6 a RLQG controller improved with a FP compensator. 0 2.7 2.23 3.17 2.61 14 0 1.65 1.86 5.38 3.99	1.67 5.12	1.68 5.13	0.92 2.86	0.92 2.85
	7 15	0 0	2.71 1.64	2.26 1.87	3.14 5.52	2.55 4.05	1.6 5.29	1.61 5.29	0.88 2.95	0.87 2.94
	8 Number 16	0 J1 0	2.7 J2 1.64	2.27 J3 1.88	3.17 J4 5.47	2.63 J5 4.05	1.66 J6 5.23	1.66 J7 5.23	0.91 J8 2.92	0.90 J9 2.90
	9 (#) 17	0 (msec) 0	2.73 (%) 1.65	2.3 (%) 1.88	3.12 (%) 5.36	2.6 (%) 3.97	1.56 (%) 5.12	1.57 (%) 5.12	0.86 (%) 2.86	0.85 (%) 2.84
	10 1 18	0 0 0	2.7 2.24 1.62	2.24 2.12 1.81	3.19 3.71 5.97	2.6 2.8 4.32	1.7 2.97 5.72	1.71 2.97 5.72	0.93 1.5 3.19	0.93 1.5 3.18
	11 2 19	0 0 0	2.72 2.22 1.63	2.25 2.08 1.83	3.13 3.87 5.79	2.52 2.92 4.24	1.54 3.16 5.55	1.54 3.17 5.55	0.85 1.59 3.10	0.84 1.59 3.08
	12 3 20	0 0 0	2.76 2.25 1.61	2.36 2.11 1.79	3.12 3.7 6.03	2.59 2.83 4.37	1.51 2.95 5.77	1.51 2.96 5.77	0.83 1.49 3.22	0.82 1.49 3.20
	13 4 21	0 0 0	2.69 2.25 1.66	2.2 2.11 1.87	3.14 3.6 5.45	2.55 2.75 4.03	1.61 2.79 5.23	1.61 2.8 5.23	0.88 1.41 2.92	0.88 1.41 2.90
	14 5	0 0	2.75 2.25	2.29 2.13	3.09 3.65	2.48 2.76	1.43 2.9	1.44 2.91	0.79 1.47	0.78 1.47
	15 6	0 0	2.73 2.27	2.26 2.16	3.18 3.57	2.67 2.74	1.67 2.78	1.67 2.79	0.91 1.41	0.90 1.41
	Earthquake Engng Struct. Dyn. 1992; 21:779-792. 2. Williams DM, Williams MS, Blakeborough A. Numerical modeling of a servo-hydraulic 16 0 2.69 2.23 3.17 2.63 1.7 1.7 0.93 0.93 7 0 2.22 2.04 3.96 3.02 3.27 3.28 1.65 1.64 J9 17 0 2.73 2.3 3.13 2.56 1.54 1.55 0.85 0.84 8 0 2.25 2.12 3.67 2.81 2.91 2.92 1.47 1.47 Number J1 J2 J3 J4 J5 J6 J7 J8 J9 (%) 1 0 2.89 2.35 3.38 2.6 1.75 1.75 0.95 18 0 2.68 2.19 3.2 2.61 1.75 1.76 0.96 0.96 9 0 2.23 2.08 3.81 2.87 3.11 3.12 1.58 1.57 (#) (msec) (%) (%) (%) (%) (%) (%) (%) (%) 0.95 2 0 2.89 2.34 3.43 2.65 1.86 1.86 1.01 19 0 2.73 2.3 3.1 2.56 1.52 1.52 0.83 0.82 10 0 2.24 2.12 3.72 2.83 2.94 2.95 1.48 1.48 1 0.2 3.56 3.55 7.54 5.09 5.95 5.97 3.16 3.16 1.01 3 0 2.92 2.41 3.31 2.54 1.57 1.58 0.85 20 0 2.7 2.27 3.15 2.6 1.62 1.63 0.89 0.88 11 0 2.23 2.07 3.82 2.89 3.12 3.13 1.57 1.57 2 0 4.36 4.69 9.89 6.42 8 8.02 4.3 4.3 0.86 4 0 2.86 2.26 3.44 2.62 1.86 1.87 1.01 21 0 2.7 2.25 3.16 2.59 1.65 1.66 0.91 0.90 12 0 2.24 2.1 3.67 2.78 2.91 2.92 1.47 1.47 3 0.5 2.75 2.4 5.28 3.51 4.16 4.17 2.19 2.18 1.01 5 0 2.9 2.38 3.38 2.59 1.73 1.74 0.94 Table A-2: vRTHS evaluation criteria for cases 2 for El Centro earthquake excitation, using 13 0 2.26 2.14 3.58 2.73 2.8 2.8 1.41 1.41 4 0.2 3.94 4.05 8.37 5.66 6.63 6.65 3.52 3.52 0.95 6 0 2.94 2.44 3.33 2.57 1.62 1.63 0.88 a RLQG controller improved with a FP compensator. 14 0 2.24 2.09 3.61 2.75 2.8 2.81 1.41 1.41 5 0.2 3.58 3.62 7.72 5.12 6.15 6.17 3.26 3.26 0.88 7 0 2.9 2.37 3.37 2.62 1.73 1.74 0.93 15 0 2.23 2.1 3.76 2.85 3.01 3.02 1.52 1.52 6 0.2 3.14 2.94 6.44 4.27 5.03 5.05 2.67 2.66 0.94 8 0 2.88 2.32 3.41 2.6 1.8 1.8 0.98 Number J1 J2 J3 J4 J5 J6 J7 J8 J9 16 0 2.25 2.16 3.64 2.77 2.88 2.89 1.46 1.46 7 1.7 4.43 4.45 14.13 8.35 12.94 12.99 6.47 6.47 0.98 9 0 2.85 2.26 3.54 2.71 2.08 2.08 1.12 (#) (msec) (%) (%) (%) (%) (%) (%) (%) (%) 17 0 2.24 2.14 3.68 2.81 2.92 2.93 1.47 1.47 8 0.7 3.58 3.83 9.48 5.81 8.17 8.2 4.25 4.24 Number J1 J2 J3 J4 J5 J6 J7 J8 J9 1.13 10 0 2.87 2.27 3.49 2.69 1.98 1.98 1.07 1 0 1.75 1.58 3.04 2.62 2.48 2.48 1.73 1.73 18 0 2.24 2.1 3.84 2.92 3.14 3.15 1.59 1.59 9 0 3.78 3.91 7.59 5.29 5.83 5.85 3.12 3.12 (#) (msec) (%) (%) (%) (%) (%) (%) (%) (%) 1.07 11 0 2.92 2.41 3.34 2.55 1.62 1.63 0.88 2 0 1.75 1.6 3.06 2.6 2.52 2.52 1.76 1.76 19 0 2.26 2.12 3.68 2.76 2.96 2.97 1.5 1.5 10 0 4.89 5.2 12.47 7.76 10.48 10.51 5.61 5.61 1 0.2 2.98 3.42 10.1 7.78 9.47 9.47 6.62 6.63 0.89 12 0 2.88 2.3 3.47 2.66 1.94 1.95 1.05 3 0 1.75 1.6 2.99 2.55 2.43 2.43 1.69 1.69 20 0 2.24 2.11 3.72 2.82 2.99 3 1.51 1.51 11 0 3.9 4.28 7.73 5.43 5.88 5.9 3.2 3.2 2 0 2.54 2.77 6.64 5.5 5.89 5.89 4.09 4.1 1.06 13 0 2.9 2.4 3.39 2.6 1.76 1.77 0.96 4 0 1.75 1.6 3.06 2.63 2.49 2.49 1.74 1.73 21 0 2.24 2.05 3.85 2.89 3.17 3.18 1.6 1.6 12 0.7 4 4.22 9.67 5.98 8.09 8.12 4.24 4.23 3 0.7 3.38 3.77 12.8 9.53 12.25 12.25 8.6 8.62 0.96 14 0 2.9 2.36 3.39 2.61 1.75 1.76 0.95 5 0 1.76 1.6 3.01 2.55 2.46 2.46 1.71 1.71 13 0.2 2.81 2.5 5.51 3.73 4.33 4.34 2.28 2.27 4 0 3.05 3.48 10.3 7.86 9.65 9.65 6.74 6.76 0.96 15 0 2.85 2.27 3.46 2.65 1.9 1.91 1.03 6 0 1.75 1.59 3.12 2.67 2.58 2.58 1.8 1.8 14 0 4.16 4.61 8.58 5.92 6.66 6.68 3.61 3.61 5 0 2.49 2.77 7.84 6.15 7.18 7.18 5.01 5.02 1.03 16 0 2.93 2.4 3.37 2.6 1.68 1.69 0.92 7 0 1.76 1.61 3.11 2.69 2.56 2.56 1.78 1.78 15 0.7 3.86 4.1 9.32 5.81 7.82 7.84 4.09 4.08 6 2.2 5.13 5.48 24.52 17.93 24.08 24.09 17.1 17.16 0.92 17 0 2.91 2.4 3.36 2.6 1.7 1.71 0.92 8 0 1.73 1.56 3.23 2.81 2.69 2.69 1.87 1.87 Number J1 J2 J3 J4 J5 J6 J7 J8 J9 16 0.5 3.26 3.28 7.64 4.82 6.37 6.39 3.34 3.33 7 0.5 3.48 3.98 12.35 9.34 11.73 11.73 8.23 8.25 0.93 18 0 2.85 2.29 3.47 2.65 1.94 1.95 1.06 9 0 1.75 1.59 3 2.55 2.42 2.42 1.69 1.69 (#) (msec) (%) (%) (%) (%) (%) (%) (%) (%) 17 0.7 3.45 3.59 8.4 5.25 7.05 7.08 3.69 3.68 8 0 3.66 4.14 12.26 9.2 11.54 11.54 8.06 8.08 1.06 19 0 2.87 2.31 3.45 2.64 1.88 1.89 1.02 10 0 1.74 1.59 3.13 2.71 2.57 2.57 1.79 1.79 1 0 1.65 1.89 5.40 4.01 5.15 5.16 2.88 2.87 18 0 4.38 4.9 9.48 6.49 7.53 7.56 4.1 4.1 9 0 2.88 3.17 7.75 6.31 7 7 4.86 4.87 1.03 20 0 2.88 2.32 3.42 2.66 1.86 1.87 1.01 11 0 1.77 1.63 2.89 2.45 2.31 2.31 1.61 1.61 2 0 1.63 1.83 5.77 4.20 5.54 5.54 3.10 3.08 19 0.7 3.63 3.88 8.97 5.54 7.6 7.62 3.96 3.95 10 0 3.17 3.54 9.3 7.22 8.56 8.56 5.96 5.97 1.02 21 0 2.9 2.37 3.38 2.63 1.76 1.77 0.96 12 0 1.75 1.59 2.96 2.53 2.38 2.38 1.66 1.66 3 0 1.67 1.91 4.92 3.73 4.63 4.63 2.59 2.57 20 0 3.58 3.61 7.04 4.94 5.35 5.37 2.86 2.86 11 0 3.76 4.21 11.61 9 10.85 10.86 7.55 7.56 0.96 Table A-1: vRTHS evaluation criteria for cases 1 for El Centro earthquake excitation, using 13 0 1.75 1.6 3.05 2.6 2.5 2.5 1.74 1.74 4 0 1.64 1.86 5.62 4.12 5.40 5.41 3.02 3 21 1.5 4.71 5.16 13.28 7.68 11.7 11.74 5.95 5.95 12 0 3.37 3.88 9.89 7.85 9.19 9.19 6.35 6.37
					25


Olivier Pujol
email: olivier.pujol@univ-lille.fr

Andrew Jensen
email: ajensen@northland.edu

Cloud-rain predator-prey interactions: analyzing some properties of the Koren-Feingold model and introduction of a new species-competition bulk system with a Hopf bifurcation

Keywords: Cloud-rain interactions, Species competition dynamics, Predator-prey model, Cloud physics

This paper deals with some properties of predator-prey cloud-rain models. Focus is put on scaling and on some mathematical features such as stability and limit cycles. Precisely, the Koren-Feingold delay differential equation model is first investigated and it is shown that it has no limit cycles. Then, by considering another point of view (i.e. species competition dynamics) for parametrizing cloud-rain processes, a system of ordinary differential equations to model these processes is formulated. Some examples are given to illustrate that this model reproduces in a realistic way the essential macroscopic behavior of a cloud-rain system. The model has a Hopf bifurcation at which certain properties of cloud-rain interactions in the model are represented. This is an important point to prepare for further examination of cloud synchronization in a cloud field by Kuramoto model, for instance.

Introduction

Clouds are complex nonlinear dynamic systems with many degrees of freedom and interactions across a vast range of spatio-temporal scales. Despite this complexity, certain non-trivial aspects of their macrobehavior are predictable without a consideration of the full complexity of the dynamic system [START_REF] Feingold | Precipitationgenerated oscillations in open cellular cloud fields[END_REF][START_REF] Koren | Aerosol-cloud-precipitation system as a predator-prey problem[END_REF][START_REF] Feingold | A model of cloupled oscillators applied to the aerosol-cloud-precipitation system[END_REF][START_REF] Koren | Exploring the nonlinear cloud and rain equation[END_REF][START_REF] Mülmenstädt | The Radiative Forcing of Aerosol-Cloud Interactions in Liquid Clouds: Wrestling and Embracing Uncertainty[END_REF]. This implies that models that capture the essential physics have the potential to contribute to our understanding of clouds and their interactions. Such models are thus always welcome, especially in conjunction with more detailed models. This is of importance in the context of the climate system (e.g. [START_REF] Glansdorff | The gap between simulation and understanding in climate modeling[END_REF] since, as is well known, clouds play a major role in the climate system (IPCC, 2013).

An example of macrobehavior that can be tackled by nonlinear dynamics is the interaction between cloud droplets and rain. Recently, [START_REF] Koren | Aerosol-cloud-precipitation system as a predator-prey problem[END_REF] (hereafter KF11) analyzed cloud-rain coupling by means of the following non-linear system of delay differential equations (DDEs):

         dH dt = H 0 -H τ 1 - αH 2 (t -T ) c 1 N d (t -T ) dN d dt = N 0 -N d τ 2 -αc 2 H 3 (t -T ) (1)
with c 1 = 2 × 10 -6 mm m -2 , c 2 = 0.3 m -1 , and α = 2 mm m -6 d -1 (see also

Table 1 in appendix). Here, H is cloud depth (in m) and N d is cloud droplet concentration (in cm -3 ). These are the two macroscopic degrees of freedom.

In addition, t is time and T is an arbitrary constant delay (both in minutes).

The constants τ 1 and τ 2 are timescales of the order of ten minutes. H 0 is a height that represents the full environmental potential for cloud development (i.e. the maximum cloud depth possible) and N 0 is the background concentration of aerosol (i.e. the maximum concentration of cloud droplets that can be reached). According to KF11, this system exhibits predator-prey like behavior. In particular, two examples of cloud behavior are provided which represent (i ) an oscillator limit cycle and (ii ) a damped oscillator.

In the first part of this paper, we explore some physical and mathematical properties of Eq.1 that we consider important from the perspective of nonlinear dynamics that have thus far not been addressed. More precisely, our analysis, which is presented in Section 2, concentrates on scale analysis, form of the equations, stability analysis, and the question of the existence of limit cycles by means of the Busenberg theorem applied to DDEs.

In the second part of this paper, we introduce a bulk model which follows the initial idea of modeling the macrobehavior of a cloud-rain system and has properties in common with species competition dynamics. In contrast to the KF11 model, it is a system of ODEs rather than DDEs. However, the model presented here reproduces the behavior of the KF11 model, and is derived from a physically-based parametrization. This model differs from the KF11 model not only in its structure, but in the variety of its dynamics which we examine through the lense of nonlinear dynamics, e.g. linear stability analysis and bifurcation theory. This exploration yields regions of parameter space in which clouds grow at the expense of rain and vice versa.

It must be emphasized that this model does not consider the full complexity of cloud-precipitation interactions; rather, we examine certain predictable elements of the macrobehavior of the full complex dynamical system.

In particular, the model we propose incorporates population dynamics in a natural way, which includes exchanges with the surrounding environment.

Our model is presented in Section 3 and examples are given in order to show that it is physically realistic and able to reproduce the macroscopic behavior of a cloud-rain system. Section 4 provides a discussion and conclusions.

2. Analysis of some physical and mathematical properties of the KF11 model

Scale or order-of-magnitude analysis

The system (1) can be solved numerically by means of the pydelay package of Python (http://pydelay.sourceforge.net/) or MATLAB's dde23 algorithm. The PyDDE solver can also be useful. These methods are based on the Bogacki-Shampine method which is a 3(2) Runge-Kutta scheme adapted to DDEs. However, before diving into numerical computations, we analyze the magnitudes of the different terms in the system of equations.

In order to have a consistent set of equations in terms of units, we replace α and c 2 in Eq. ( 1) by α = 10 -3 α/(60 × 24) and c 2 = 10 -6 c 2 (see appendix). Typically, for cloud height, H 0 ∼ 500 m and so, by definition of H 0 , 0 < H < H 0 . For cloud droplet concentration, N 0 ∼ 400 cm -3 and so, by definition of N 0 , N d < N 0 . A reasonable minimum for N d (N d,min ) is about 50 cm -3 . Hence, considering the values of the constants involved in Eq. (1), i.e. τ 1 = τ 2 ∼ 60 min, we find:

H 0 -H τ 1 < H 0 τ 1 ≈ 8 m min -1 α H 2 c 1 N d < α H 2 0 c 1 N d,min ≈ 3 500 m min -1
and

N 0 -N d τ 2 < N 0 τ 2 ≈ 7 cm -3 min -1 α c 2 H 3 < α c 2 H 3 0 ≈ 0.5 cm -3 min -1
It follows that, since H and N d are physically expected to maintain orders or magnitude as mentioned above during the evolution of the cloud rain system, that α H 2 /(c 1 N d ) and (N 0 -N d )/τ 2 are the dominant terms in Eq. ( 1). This is particularly clear for the equation governing the evolution of H. Concerning the equation for dN d /dt, having a ratio close to 1 implies that c 2 be 10 4 or 10 5 higher than the value given in Table 1. Consequently, the evolution of N d is either an increasing or a decreasing exponential function according to the sign (positive or negative, respectively)

of N 0 -N d .
The examples chosen by KF11 have positive signs. Precisely, of any constant delay in the equation. Some numerical computations with the numerical scheme mentioned above have confirmed the evolution described above (Fig. 1). Nonetheless, we have obtained some oscillatory behavior in rare cases by varying arbitrarily some constants in Eq. (1) as the constant c 1 around 10 -3 mm m -1 and the constant c 2 around 10 4 m -1 or 10 5 m -1 . We have observed that these oscillations are, in addition, extremely sensitive to the value of c 2 chosen since they disappear if c 2 is changed even slightly (we recover the evolution described above).

N d (t) -N 0 = [N d (0) -N 0 ] exp(-t/τ 2 ),

Form of the system of equations

It is worth noting that the system (1) of two equations is not formally similar to a predator-prey system, stricto sensu. Indeed, in its basic form, the Lokta-Volterra system is the following (x and y are the two degrees of freedom):

       dx dt = ax -bxy dy dt = -cy + dxy (2)
with {a, b, c, d} a set of four positive constants. The signature (+, -| -, +), or sometimes (+, -| + , -), is also present in more elaborate versions of Lokta-Volterra models (e.g. [START_REF] Lipowski | Nonequilibrium phase transition in a lattice prey-predator system[END_REF] since it is a characteristic of species-competition models, with or without predation. However, the system 1 has the form: dH/dt = -AH-BH

H τ 1 + α H 2 c 1 N d > H 0 τ 1 ∼ 10 m min -1 and N d τ 2 + α c 2 H 3 > N 0 τ 2 ∼ 1 -10 cm -3 min -1
So, with regard to the typical orders of magnitude of the two degrees of freedom (H ∼ 400 m, N d ∼ 100 cm -3 ), care is needed in the choice of the constants (i.e. in the model parametrization) in order to avoid always decreasing functions H(t) and N d (t).

Stability analysis

To determine further some of the properties of the system (1), we briefly discuss the stability of the model (with no delay) near its equilibrium points. These are independent of any delays, so that a system of DDEs has the same equilibrium points as the corresponding system (zero delay) of ODEs. However, investigating the stability for DDEs is a more complex task (e.g. Engelborghs et al., 2000, and references therein) which is beyond the scope of this paper (some elements can be found in Koren et al. (2017, Section IV)). Nonetheless, it is still instructive to analyze stability in the simplest case of no delay. A quick calculation yields that the KF11 system has an equilibrium point given by

N e d = N 0 -α c 2 τ 2 (H e ) 3
and

H e = c 1 N e d 2α τ 1   -1 + Ç 1 + 4α τ 1 H 0 c 1 N e d å 1/2   (3)
Physically, because N e d (and H e ) must be positive, the value of α c 2 is constrained: α c 2 < N 0 /[τ 2 (H e ) 3 ]. The behavior of Eq. ( 1) near this equilibrium point can be obtained by evaluating the Jacobian matrix J at (H e , N e d ) and finding its eigenvalues. These are given by the well-known formula:

λ ± = 1 2 ï T J ± Ä T 2 J -4∆ J ä 1/2 ò
where T J and ∆ J are respectively the trace and the determinant of [J]. A straightforward calculation yields:

T J = - Ç 1 τ 1 + 1 τ 2 + 2α H e c 1 N e d å and ∆ J = 1 τ 2 Ç 1 τ 1 + 2α H e c 1 N e d å + 3α 2 c 2 (H e ) 4 c 1 (N e d
) 2 Since T J < 0 and ∆ J > 0, the equilibrium point (H e , N e d ) is a stable node or a stable focus (spiral) according to the sign of (T 2 J -4∆ J ). We note that if the value of α c 2 does not satisfy the above constraint, then T J can be negative, which means that Eq. ( 1) describes an unstable dynamical system.

In addition, since the sign of T J never changes, there is no local bifurcation.

A similar analysis that corroborates the above result has been performed

recently by [START_REF] Jiang | Aerosol Replenishment and Cloud Morphology: A VOCALS Example[END_REF].

Existence of limit cycles

One of the important points of this paper is the possibility of limit cycle solutions of the KF11 delay differential equations (DDEs). This needs rigorous examination. First, we note that the KF11 equations with constant delay T do not have any limit cycle solutions with period T . To see this, we consider either the Bendixson-Dulac theorem or the special case known as

Bendixson's criterion (e.g. [START_REF] Minorsky | Nonlinear Oscillations[END_REF][START_REF] Glansdorff | The gap between simulation and understanding in climate modeling[END_REF].

To illustrate the use and consequences of this theorem to the problem at hand, we first rewrite the KF11 system: dH

/dt = f [H(t -T ), N d (t -T )] and dN d /dt = g[H(t -T ), N d (t -T )],
where the functions f and g are given by the rhs (right hand sides) in Eq. ( 1). The Bendixson-Dulac theorem may be applied to this system by considering the sign of the expression

∂ H (φf ) + ∂ N d (φg), where φ(H, N d ) = 1 (this case is Bendixson's criterion). Now, suppose that this set of equations has a non-constant T -periodic solution [H(t), N d (t)]. Then, it follows that dH/dt = f [H(t), N d (t)] and dN d /dt = g[H(t), N d (t)]. A simple calculation shows that ∂ H (f ) + ∂ N d (g)
is always negative. We can therefore conclude by the Bendixson-Dulac theorem that no solutions of period T exist in the simply connected region defined by positive values of H and N d .

Since the KF11 equations are DDEs, a customary application of the Bendixson-Dulac theorem to prove the non-existence of limit cycles is inadequate. However, this theorem has been generalized (Busenberg and van den Driessche, 1993, Section 4) to higher dimensional ODEs, including DDEs like those proposed by KF11, for instance. This generalized "Bendixson-Dulac" theorem extends the concept of limit cycle to that of simple loop solutions,

i.e. any continuous solution of the system of equations whose orbit contains a closed curve. Here, we apply this generalized theorem to the KF11 DDEs in order to demonstrate the absence of loop solutions.

To illustrate the use and consequences of this generalized theorem to the problem at hand, especially theorem 4.1 in Busenberg and van den Driessche (1993), we first rewrite the KF11 system:

       dH dt = F 1 (H, N d ) h t (H t , N d,t ) + l 1 (H, N d ) dN d dt = F 2 (H, N d ) k t (H t , N d,t ) + l 2 (H, N d )
where [START_REF] Busenberg | A method for proving the nonexistence of limit cycles[END_REF] gives conditions under which the KF11 DDEs have no simple loop solution. To verify the first condition, we define the vector function g(H, N d ) = (g 1 , g 2 , 0). Here,

F 1 (H, N d ) = -α/c 1 , F 2 (H, N d ) = -αc 2 , h t (H t , N d,t ) = H 2 t /N d,t and k t (H t , N d,t ) = H 3 t are time-delayed functions (the subscript "t" stands for "time-delayed"), l 1 (H, N d ) = (H 0 -H)/τ 1 , and l 2 (H d , N d ) = (N 0 -N d )/τ 2 . Theorem 4.1 in
g 1 = α c 2 A 3 + (N 0 -N d )/τ 2 , g 2 = -α A 2 /(c 1 B) -(H 0 -H)/τ 1 ,
A and B being constants. For notational simplicity, let (u, v) represent the coordinates of the curve in Theorem 4.1 of [START_REF] Busenberg | A method for proving the nonexistence of limit cycles[END_REF]. Then:

g(H, N d ) • Ç -αu 2 c 1 v + l 1 , -αc 2 u 3 + l 2 , 0 å = α 2 c 2 A 2 u 2 c 1 Ç u B - A v å + αc 2 l 1 Ä u 3 + A 3 ä - αl 2 c 1 Ç A 2 B + u 2 v å ≤ 0
This follows from the following observations. The third term always contributes negatively, the second positively, to the sum. The first contributes negatively whenever uv < AB. Because c 1 is small, the first and third terms dominate the second so that, whenever uv < AB, the entire sum is negative.

Since A and B can be taken as large as wished, the region is finally R 2 + . To verify the second condition, we calculate

∂ H g 2 -∂ N d g 1 = 1/τ 1 + 1/τ 2 > 0 for
any solution (H, N d ) in the region specified above. So, an application of the Busenberg-van den Driessche theorem implies that there is no simple loop solution traversed in the clockwise sense (the analogue for DDEs of the limit cycle of ODEs) for the positive values of H and N d . A similar argument with vector function -g works in the opposite sense as well. In other words, there is no solution of the KF11 system which enclose a smooth oriented simple closed curve.

A new species-competition bulk model

Description

In order to retain the simplicity of the KF11 model with ODEs instead of DDEs, we propose an approach which incorporates population competition, interaction, and carrying capacities in a natural way. Here, the populations are given by cloud water content L c (g cm -3 ), rain water content L r (g cm -3 ), and cloud droplet number concentration N d (cm -3 ). We emphasize that these variables are considered here as macroscopic averages and that our model Other terms can be added to the equations, such as self-collection for instance, but to retain some measure of simplicity we use only the terms already mentioned and which we describe further. This parametrization is a straightforward consequence of the stochastic collection equation (STE) with polynomial kernel. For this reason, it seems to us that this approach is more rigorous than many other methods which consist in using heuristic (or empirical) parametrizations. Details on the STE can be found in [START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF] and [START_REF] Pruppacher | Microphysics of Clouds and Precipitation[END_REF] and references therein. It is worth mentioning that the goal of [START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF] has been to fill a gap between heuristic parametrizations and detailed microphysical (computationally expensive) schemes. The terms representing autoconversion and accretion are given by the following expressions (for a comprehensive overview see the references):

-Autoconversion: the contribution to dL c /dt and dL r /dt is given by

dL c dt auto = -KL 4 c N -2 d = - dL r dt auto (4)
-Accretion: the contribution to dL c /dt and dL r /dt is given by

dL c dt accret = -k r L c L r = - dL r dt accret . (5)
We note that these expressions are adapted from the [START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF] parametrization by considering the variables as bulk, cloud-scale variables. The various parameters in the equations above are defined as follows:

K = [k c /(20x * )] (ν + 2)(ν + 4)(ν + 1) -2
, where x * is a cloud drop mass separating the cloud droplets from raindrops, k c (k r ) is a constant from the cloud (rain) water kernel, and ν is the shape parameter of the gamma distribution. See [START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF] for complete details. By combining these equations, neglecting spatial dependence and adding source/sink terms for cloud and rain water contents, we arrive at the following system of bulk coupled first-order differential equations for L c , L r and N d :

               dL c dt = A c L c -B c (N d )L 4 c -k r L c L r dL r dt = -A r L r + B c (N d )L 4 c + k r L c L r dN d dt = f (N d ,L c ,L r ) (6)
Here, B c (N d ) = KN -2 d , and A c and A r are two positive constants that can be chosen so that the source/sink terms represent the background meteorological conditions. In particular, A r represents the rain out process. These two constants give the timescales, τ c = 1/A c and τ r = 1/|A r |, of cloud water content and rain water content evolutions respectively. In addition, f (N d , L c , L r ) is a function which can be defined according to the question being investigated. Here, we choose to define f as follows:

dN d dt = A c (N 0 -N d ) - 4 3 k c L 2 c -k r L r N d . (7)
On the rhs, the first term represents a supply of cloud droplets from the surroundings, N 0 being the background aerosol concentration that feeds the system (nucleation). Considered alone, this term causes an exponential increase (decrease) of N d with N 0 as horizontal asymptote if initially N d is lower (higher) than N 0 . This represents cloud droplet concentration and when considered alone it tends toward a constant background aerosol population concentration. Moreover, this term (or similar ones) should be dominant in the very early stages of the cloud system evolution. Indeed, at the initial time (and close to t = 0), N d is zero or very small (clear air condition) and must increase enough in order that the cloud can further evolve and produce rain.

The two other terms in Eq. ( 7) come from Seifert and Beheng (2001, Eqs.

A-5, A-6, and A-9) and represents respectively cloud droplet self-collection and accretion.

It has to be recalled that the state variables considered (L c , L r , and N d ) in this bulk model are macroscopic averages over the whole cloud. They can also be seen as local variables over a volume element of the cloud for which they are more or less uniform. A further step would be to introduce in the equation an explicit inside-cloud location (e.g. altitude z) dependence.

Nonetheless, as showed below, the present model describes some macroscale features of could-precipitation coupling.

Comment on the form of the system of equations

We note that the first two equations in (6) constitute a predator-prey system modified to include autoconversion. To see this, observe that the constants {a, b, c, d} in the Lotka-Volterra equations (Eq. 2) have the following formal correspondence with the model presented here:

{a = A c , b = -d = k r , c = A r }.
Hence, the source/sink and accretion terms form a standard predator-prey system when considered alone. The quartic terms not only represent autoconversion of cloud droplets to rain given in the Seifert and Beheng ( 2001) parametrization, but also ensure that the cloud droplet population does not exhibit unbounded growth in the absence of competition.

Thus, it plays a similar role to the carrying capacity term used in population dynamics to modify some of the unrealistic features of the original Lotka-Volterra system.

The entire system of equations is not rigorously equivalent to a threespecies competition model. As already discussed, the first two equations constitute a modified Lotka-Volterra system. The autoconversion terms which modify the traditional predator-prey system have coefficients driven by the third equation, viz. that for N d . The interpretation of the varying coefficients is that the carrying capacities of the water contents vary according to the state of N d , which is more realistic for clouds.

It is worth mentioning that L c and L r are, for physical reasons, necessarily bounded. So, in virtue of the Poincaré-Bendixson theorem, the first two equations of system (6) have a solution that either converges towards a limit or presents an asymptotic behaviour that can take the form of a limit cycle.

Boundedness has not been proved in a mathematical sense but, rather, we rely on numerical evidence (see some examples below, in Subsection 3.4) and on the fact that the present system of equations has much in common with other predator-prey models which are known to have bounded solutions. Obviously, if k r and B c are zero, the evolution of L c is exponentially divergent.

The terms involving non-zero values of k r and B c prevent unbounded evolution within certain parameter ranges, as is the case in predator-prey models.

It is far beyond the scope of this paper to address rigorously (i.e. mathematically speaking) boundedness of solutions, but that could be formulated and explored deeply in a future specific work.

Stability analysis

Equilibrium points

The system

       dL c dt = A c L c -B c (N d )L 4 c -k r L c L r dL r dt = -A r L r + B c (N d )L 4 c + k r L c L r (8)
has two equilibrium points X e = {L (e) c , L (e) r }. In this section we treat N d as a parameter. The first one is trivial: X (1) e = 0. There is no cloud nor a fortiori rain. Using Eq. ( 7), N (e) d = N 0 , which means that the atmospheric volume considered has a constant CCN (cloud condensation nuclei) loading (that does not allow cloud formation). The second equilibrium point X (2) e is such that:

A c L (e) c = |A r |L (e) r and B c (N d ) |A r | 3 A 4 c Ä L (e) r ä 3 + k r A c L (e) r -1 = 0 (9)
Since the two first terms of this third order polynomial are positive, there is only one real positive solution L (e) r -the other two are complex conjugates. The expression of N (e) d follows from those of L (e) c and L (e) r using Eq. ( 7) Plugging L (e) c and L (e) r into Eq. ( 7) yields:

N (e) d = Ç A c + k r A c |A r | L (e) c å -1 ñ A c N 0 - 4 3 k c Ä L (e) c ä 2 ô (10)
The eigenvalues of the jacobian matrix J of the system (6) are:

λ ± = 1 2 ï T J ± Ä T 2 J -4∆ J ä 1/2 ò
where the trace T J and the determinant ∆ J of [J] can be written as follows:

T J = Tr(J) = A c -|A r | + k r Ä L (e) c -L (e) r ä -4B c (N d ) Ä L (e) c ä 3 (11)
and

∆ J = det(J) = -A c |A r | + k r Ä A c L (e) c + |A r |L (e) r ä + 4|A r |B c (N d ) Ä L (e) c ä 3 (12)
They must be evaluated for each of the two equilibrium points.

Stability analysis around X (1)

e

In this case,

T J = A c -|A r | and ∆ J = -A c |A r | < 0. The two eigenvalues have opposite sign, λ + λ -< 0, with λ + = A c > 0 and λ -= -|A r | < 0. The state X (1)
e is thus a saddle point.

Stability analysis around X (2)

e

In this case

T J = (3A c + |A r |) k r L (e) r A c -1 and ∆ J = -2A c |A r | k r L (e) r A c - 3 2
It follows that the equilibrium point

X (2) e = Ä L (e) c , L (e) r ä is stable if: k r L (e) r A c < 1 i.e. L (e) r < A c k r and L (e) c < |A r | k r (13) Otherwise, for k r L (e) r /A c > 1, X (2)
e is an unstable equilibrium point.

The identity k r L (e) r /A c = 1, which separates a stable regime from an unstable one, has the straightforward consequence that B c N The other characteristic value is k r L (e) r /A c = 1.5. Below (above) this value, ∆ J is positive (negative). So

X (2) e is a saddle point if L (e) r = 1.5A c /k r .
The nature (node, focus or spiral) of X (2) e is provided by the sign of T 2 J -4∆ J :

T 2 J -4∆ J =   (3A c + |A r |) 2 k r L (e) r A c -1 2   + 4|A r |A c 2 k r L (e) r A c -3 If k r L (e) r /A c = 1, T 2 J -4∆ J = -4|A r |A c < 0 which means that X (2)
e is a focus.

For k r L (e) r /A c = 1, the sign of T 2 J -4∆ J depends on the number N τ = |A r |/A c through the second order polynomial equation:

(3A c + |A r |) 2 k r L (e) r A c 2 -2 k r L (e) r A c î (3A c + |A r |) 2 -4A c |A r | ó +(|A r | -3A c ) 2
The two solutions are easy to find:

Z ± = (3A c + |A r |) 2 -4|A r |A c (3A c + |A r |) 2 ± 2 [(9A 2 c + |A r |A c ) (|A r | 2 + A c |A r |)] 1/2 (3A c + |A r |) 2
Three cases can be distinguished at first sight:

-|A r | = A c : so Z ± = (3 ± √ 5)/4 > 0. Between these two values, T 2 J - 4∆ J < 0, so that X (2)
e is a focus. Outside of this interval, we have a node

-|A r | A c : so Z ± = 1 ± 2 (A c /|A r |) 1/2 > 0. Same as previously. -|A r | A c : so Z ± = 1 ± (2/3) (|A r |/A c ) 1/2 > 0.
Same as previously.

Andronov-Hopf bifurcation

It follows from what precedes that the eigeinvalues of the jacobian matrix cross the imaginary axis for B c = 0. At this point, they are pure imaginary numbers, λ ± = i(A c |A r |) 1/2 . There, the system has a limit cycle. Just after

(B c < 0), the spiral is unstable whilst it is stable just before (B c > 0).
This behavior is precisely a Hopf bifurcation: an unstable focus gives birth to a stable focus (and vice versa) through a limit cycle. The existence of such a bifurcation is an important point for further investigations of cloud organization in a coupled-cloud field at the mesoscale [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF].

Examples

We now illustrate the behavior of this set of three equations with some examples. The values chosen for k r and k c are those given in Seifert and Beheng ( 2001) which follow from approximations to the collection kernel in [START_REF] Long | Solutions to the Droplet Collection Equation for Polynomial Kernels[END_REF], i.e. k c ≈ 5.66×10 11 cm 3 g -2 min -1 and k r ≈ 3.47×10 5 cm 3 g -1 min -1 .

The value of x * = 2.6 × 10 -7 g, following [START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF]. We take ν = 2, which is a typical value for warm clouds (e.g. , Pruppacher and Klett, 1997, chap. 2). With these values, the parameter K ≈ 2.905 × 10 17 cm 3 g -3 min -1 .

In addition, we take A c = 0.02 min -1 and A r = 0.1 min -1 . So, cloud water content and rain water content are expected to evolve on the timescales τ c = 50 min and τ r = 10 min. The latter is lower than the former since raining out is a rapid process compared to the evolution of cloud water content. These values are chosen for illustrative purposes.

First example: periodic behavior

We choose as initial conditions N d (0) = 10 cm -3 , L c (0) = 10 -9 g cm -3 , L r (0) = 0, and N 0 = 50 cm -3 . The system exhibits oscillatory behavior after approximately 300 min, with a period of about 200 min, which is physically realistic (Fig. 2a). Both N d and L c increase initially (the former faster than the latter) until rain appears. Soon thereafter, rain water content peaks and droplet number and mass concentration decrease drastically in a short interval of time. Minimum values of L c are reached and rain water content decreases to zero. Then N d and L c start to increase again, and so on. This periodicity appears clearly as a limit cycle on the plot (L c , L r ) displayed in Fig. 2b. It is interesting to mention that peaks in rain water content are delayed from L c peaks by about 20 -30 min, which is a very realistic value for rain to form from a cloud droplet by collection (Fig. 2c). We also note that similar oscillations and limit cycles are obtained for other initial conditions (not shown). Note that the model presented here represents well a timescale of precipitation production that has to be explicity and arbitrarily accounted for via delays in the more complicated DDE framework.

Second example: damped behavior

If N 0 is lower than N d (0), the feeding of the system is not sufficient to maintain the cloud-rain system in a harmonic (or quasi) oscillation and a damping is observed (Fig. 3): N d , L c , and L r still oscillate but their respective amplitudes decrease progressively until a steady state is reached. Initially N d decreases exponentially as expected, and then decreases more rapidly when L c and L r peak for the first time.

Damping can also be obtained by introducing in the model a decreasing evolution for the CCN supply, for instance N 0 = 50 exp(-0.002 t). This means that the initial CCN burden is divided by around 2.72 in 500 min, i.e. 8 hours. Keeping the same initial conditions as previously, we get the behavior displayed in Figure 4, which is typical of a damped oscillator: N d increases until a certain time t (approximately 120 min) at which the exponential function becomes significant. Then N d decreases dramatically. The same behavior is seen for L c and L r which increase until t = 500 min and then exhibit decaying amplitude oscillations that eventually (t > 1 500 min) die out. Again, we note that the peaks of L r are delayed with respect to those of L c by about 20 -30 min.

The same damping behavior is observed (not shown) for different initial conditions, in particular if L r = 0, provided the values are realistic.

Discussion and conclusion

Certain aspects of a system's bulk behavior arise out of many smallscale interactions. These aspects may be investigated with simple models that faithfully represent the bulk behavior, but do not consider directly the full complexity of the system. In conjunction with more detailed models and simulations, models which represent bulk behavior can also enhance our understanding of complex dynamical systems. The model presented here, which has aspects of species competition in biology, is a tool with which to examine the complexity of cloud-rain interactions in a way that highlights certain aspects of their behavior in a simple but realistic manner. This model is not meant to supplant more detailed simulations of cloud-precipitation interactions but to represent certain predictable elements of the macrobehavior of the full complex dynamical system.

In the first part of the paper, we examined some important mathematical properties of the KF11 DDEs model which is a similar model devoted to analyze the bulk behavior of a cloud-rain system. The focus has been put on scaling and limit cycles or loop solutions. In particular, it has been shown, through the Busenberg theorem that this model has no proper limit cycle solutions. The model we have introduced in this paper (second part) is another point of view which retains the idea of modeling cloud-precipitation interaction through population dynamics and exhibits predator-prey behavior with rain as predator and cloud droplets as prey. The stability analysis has shown that our model has limit cycles and a Hopf bifurcation.

In particular, the model formulated here is a system of ODEs rather than DDEs. It is derived from a physically-based parametrization. Several examples illustrate that our model reproduces realistically the behaviour of an unique cloud-rain system. In particular, rain water content peaks are delayed by about 20 -30 min with respect to cloud water content peaks.

Some of the novel properties of the model include regions of parameter space in which clouds grow at the expense of rain and vice versa as determined from bifurcation analysis.

Further properties of our model such as ability to represent various physical situations, and sensitivity tests to parametrization schemes similar to those performed by [START_REF] Wacker | Competition of precipitation particles in a model with parametrized cloud microphysics[END_REF] for mixed and ice clouds, will be explored in more detail in future work. Due to its properties, this cloud-rain model is intended to be further developed for examination of cloud organization at larger scales in a coupled-cloud field.

and not [L] -1 (cf. Table 1). Indeed, the well-established relationship is given by LW P (g m -2 ) = 0.5 c w H 2 , where c w = 2 × 10 -6 kg m -4 (e.g. [START_REF] Geoffroy | Relationship between drizzle rate, liquid water path and droplet concentration at the scale of a stratocumulus cloud system[END_REF]. Actually, c w = ρ w c 1 , where ρ w ≈ 10 3 kg m -3 is liquid water density and c 1 = 2 × 10 -9 m -1 -in Tab. 1, c 1 = 2 × 10 -6 mm m -2 , i.e. c 1 (m -1 ) = 10 -3 c 1 (mm m -2 ). The relationship between LW P and H we get is thus:

LW P (g m -2 ) = 1 2 c w (kg m -4 ) H 2 (m 2 ) = 1 2 ρ w (kg m -3 ) c 1 (m -1 ) H 2 (m 2 ) ( 14
)
If we look now at Eq [5] and the units chosen, the decreasing rate of LW P due to rain cannot equal exactly the rainrate R (in mm d -1 ). In fact, dLW P/dt is proportional to the precipitation flux density (or rain current J r = -ρ w R).

The coefficient of proportionality is determined by the choice of units. If, as in KF11, we express time in minutes (it is the most natural timescale for cloud macroscopic physics), and keep LW P in g m -2 as above, we obtain:

dLW P dt (g m -2 min -1 ) = -ρ w (g m -3 ) R (m min -1 ) = -10 3 ρ w (kg m -3 ) R (m min -1 ) = - 10 3 60 × 24 ρ w (kg m -3 ) R (m d -1 ) = - 1 60 × 24 ρ w (kg m -3 ) R (mm d -1 ) (15)
since R (m min -1 ) = 10 -3 R (mm d -1 )/(60 × 24). Combining Eqs. ( 14) and (15) yields:

dH dt (m min -1 ) = dLW P/dt (g m -2 min -1 ) c w (kg m -4 ) H (m) = - ρ w (kg m -3 ) R (mm d -1 ) (60 × 24) c w (kg m -4 ) H (m) = - 1 60 × 24 R (mm d -1 ) c 1 (m -1 ) H (m) = - 10 3 60 × 24 R (mm d -1 ) c 1 (mm m -2 ) H (m) (16)
It is possible, although not common, to divide LW P by ρ w and then to de-370 fine a volumetric liquid water path: LW P v (m) = 10 -3 LW P (g m -2 )/ρ w (kg m -3 ).

371

In this case, Eqs. [2] and [5] of Table 1 are correct provided that, for Eq. [5], R be expressed in m min -1 , or equivalently, considering 10 -3 R (mm d -1 )/(60× 24). However, we have not found any explicit mention of LW P in m or mm in KF11 and further papers; there LW P is explicitly expressed in g m -2 several times. The KF11 model also employs two parametrizations which relate R, N d , and dN d /dt (see Eqs. [4] and [9] in Table 1). The units of the latter two quantities are naturally cm -3 and cm -3 min -1 , respectively. According to the units used in KF11, the rhs of Eq. [4] is in mm m -3 cm 3 d -1 , while the lhs (left hand side) is in mm d -1 . Although Eq. [4] is dimensionally correct ([L][T ] -1 , where [T ] is the time dimension), the units of the rhs and lhs are not the same. Indeed, the relation R = αH 3 /N d should have N d in m -3 when H is in m. If we choose to express N d in cm -3 , then we must multiply H 3 /N d by 10 -6 . With the unit of α, R is thus in mm d -1 : R (mm d -1 ) = 10 -6 α (mm m -6 d -1 ) H 3 (m 3 ) N -1 d (cm 3 ) (17)

Thus, for Eq. ( 16), we obtain: dH dt (m min -1 ) = -10 -3 60 × 24 α (mm m -6 d -1 ) H 2 (m 2 ) c 1 (mm m -2 ) N d (cm -3 ) (18)

In addition, Wood (2006) (ref. 20 in KF11) gives a value of c 2 which does not correspond to that used in KF11. Aside from the conflation of rain current and rainrate in [START_REF] Wood | Rate of loss of cloud droplets by coalescence in warm clouds[END_REF] -this article seems to use P as rainrate and divides it by liquid water density (ρ w ), whereas P should be the rain current and P/ρ w the rainrate -, the coefficient of proportionality between dN d /dt and R is 3E 0 /4 = 3 × 10 3 m -1 since E 0 = 4 × 10 3 m -1 . This value is different from the 0.3 m -1 used in KF11 (assuming that P , in [START_REF] Wood | Rate of loss of cloud droplets by coalescence in warm clouds[END_REF], is rainrate, and not rain current, which is very few probable due to a dimensional analysis of the Eq. 14 of [START_REF] Wood | Rate of loss of cloud droplets by coalescence in warm clouds[END_REF], c 2 would be 3 m -1 , i.e. the previous value of c 2 divided by ρ w , which is again different from the value of Table 1). Because, in Eq. [9], dN d /dt has to be in cm -3 min -1 , we have:

dN d (cm -3 min -1 ) dt = c 2 N d (cm -3 ) R (m min -1 ) (19)
with c 2 = 3 × 10 3 m -1 .

Then, because R (m min -1 ) = 10 -3 R (mm d -1 )/(60 × 24), we obtain, by substituting the rainrate of Eq. ( 17) into Eq. ( 19): dN d (cm -3 min -1 ) dt = 10 -9 60 × 24 c 2 (m -1 ) α (mm m -6 d -1 ) H 3 (m 3 ) (20)

  where N d (0) = N d (t = 0). It ensues that N d evolves exponentially towards N 0 . With respect to dH/dt, the dominant term is always negative and, unless N d becomes very small, it will remain the dominant term in the evolution of H. Hence, H is a decreasing function which must stop at H = 0 by definition. This scale analysis, based only on the physical realistic orders of magnitude of H and N d is independent

  captures some of the macroscale properties of cloud-precipitation interactions by means of three relatively simple equations in three variables, rather than a fine-scale description with complete microphysical details. The model consists of three coupled differential equations, one each for L c , L r , and N d . The system includes three main processes that contribute to the cloud budgets of L c , L r , and N d , namely: (1) sources/sinks which represent exchanges with the surroundings, and two internal processes, (2) autoconversion of cloud water to rain water, and (3) accretion of cloud water by rain water.The various terms in the equations were suggested by the[START_REF] Seifert | A double-moment parameterization for simulating autoconversion, accretion and selfcollection[END_REF] parametrization for simulating autoconversion and accretion.

  .e. B c = 0 since L (e) c = 0. Physically, this critical value means that the net rate of the autoconversion process is zero: formation of a raindrop from cloud droplets only is compensated over a given time interval by the formation of cloud droplets from a raindrop only. Thus, the system maintains constant rain and cloud water contents whose values -see Eq. (13)depend only on the exchange rates with the surroundings (A c and A r ) and on the cloud-to-rain reaction constant (k r ). Below the critical value (B c < 0), the autoconversion process results in net production of cloud droplets from raindrops. It constitutes a nonlinear amplification term for cloud droplets and a nonlinear damping term for raindrops, i.e. disturbances from equilibrium grow. In other words, the cloud is growing such that cloud droplets become more numerous while the number of raindrops vanishes. This behavior is reversed above the critical value (B c > 0): the non linear term compensates any deviation from the equilibrium point (stability). In the case B c < 0, without any cloud-rain interaction (k r = 0), cloud (rain) water content is an increasing (decreasing) function.

Fig. 1 -Fig. 2 -Fig. 3 -Fig. 4 -

 1234 Fig. 1 -Evolution of H and N d according to system (1) with H 0 = 530 m, N 0 = 180 cm -3 , τ 1 = τ 2 = 60 min, and T = 12 min. Initial conditions are H(t = 0) = 200 m and N d (t = 0) = 50 cm -3 . We have used the pydelay package of Python.

  2 /N d and dN d /dt = -CN d -DH 3 with {A, B, C, D} a set of four positive constants. So, system (1) does not have the general form of a Lokta-Volterra system since the signature is fully negative, (-, -| -, -). It follows that Ḣ and Ṅd are negative if (H, N d ) always satisfy the two inequalities:

Acknowledgements

We are very grateful to the Fulbright Committee which, by means of funding for United States-France scientific researcher exchanges, facilitates our collaboration and permits our investigations in Statistical Physics for Atmospheric Cloud Systems (SPACS group).

Appendix 1: Unit analysis

The system (1) is derived from four empirical equations which we have summarized in Table 1 together with their numbering in KF11.

Eq. number Equation Constants in KF11

[2] LW P = 0.5 c 1 H 2 c 1 = 2 × 10 -6 mm m -2 [5] dLW

Tab. 1 -Equations and constants in KF11 discussed here. In KF11, LW P and R are liquid water path and rainrate, respectively. The corresponding units are explicitely g m -2 and mm d -1 .

Eq. [2] connects liquid water path (LW P , in g m -2 ) with cloud depth (H, in m). The dimension of LW P is [M ][L] -2 (where [M ] and [L] are mass and length dimensions), which implies that the dimension of

Finally, with this set of units, i.e. with H and N d in m and cm -3 respectively, and time in minutes, the consistent (in terms of units) system of equations we get is:

consider that LW P is a volumetric liquid water path, as it might be assumed as suggested above, and R in mm d -1 , the factor 10 -3 is replaced by 10 3 . Tab. 2 -Values of the constants in [START_REF] Geoffroy | Relationship between drizzle rate, liquid water path and droplet concentration at the scale of a stratocumulus cloud system[END_REF] When using this parametrization, it is necessary to first use Eq. ( 14) to change LW P into H.

(iii ) Finally, a parametrization of the delay term dN d /dt by taking

This follows from Eq. ( 14). The equation for dH r /dt can then be substituted into this equation to complete the parametrization. This equation with delay 20 is finally substituted for the delay term in the second equation in the KF11 system.

When using these different parametrizations, a set of coupled first-order DDEs similar to Eq. ( 1) is obtained. The units must be made consistent when the various empirical relationships are employed. Namely, H is in m, N d is in cm -3 , and time is expressed in minutes. For the numerical tets we have performed, these parametrizations have given results similar to those in Figure 1. However, we have not made a complete test of sensitivity of the model to changes in parametrization since it is not the purpose of this paper. Our goal here is to suggest other possibilities to parametrize cloudrain interaction processes that could be investigated more deeply in future research.