Said Ouala

Ananda Pascual

Ronan Fablet

RESIDUAL INTEGRATION NEURAL NETWORK

Keywords: Dynamical systems, Data-driven models, Neural networks, Forecasting, Runge-Kutta methods

In this work, we investigate residual neural network representations for the identification and forecasting of dynamical systems. We propose a novel architecture that jointly learns the dynamical model and the associated Runge-Kutta integration scheme. We demonstrate the relevance of the proposed architecture with respect to learning-based state-of-the-art approaches in the identification and forecasting of chaotic dynamics when provided with training data with low temporal sampling rates.

PROBLEM STATEMENT AND RELATED WORK

The modeling of physical dynamics is a critical issue. While model-driven strategies based on the definition of ordinary differential equations (ODE) governing the observable phenomena are the classic frameworks to address such a problem [START_REF] Evensen | Data Assimilation[END_REF]. Limitations in terms of numerical complexity [START_REF] Hong | Next-generation numerical weather prediction : Bridging parameterization, explicit clouds, and large eddies[END_REF] as well as the ability to better relate models and observation data for poorly-resolved processes [START_REF] Van Leeuwen | Nonlinear data assimilation in geosciences : an extremely efficient particle filter[END_REF] open the venue for datadriven representations as an appealing alternative [START_REF] Tandeo | The analog data assimilation : application to 20 years of altimetric data[END_REF], where one can determine a representation of processes of interest directly from data, especially with the increasing availability of representative data collections.

Regarding the data-driven identification of dynamical systems, one may distinguish representations based on physical priors and machine learning. The first category comprises polynomial representations [START_REF] Paduart | Identification of nonlinear systems using Polynomial Nonlinear State Space models[END_REF] and sparse regression frameworks [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. They are particularly appealing in explicitly relating the learnt representation and the associated physical operator. Such methods may, however, fail in representing complex non-linear systems for which strong priors are not available. Machine learning methods can greatly broaden modeling capabilities often however at the expense of a lack of interpretability. A large interest has for instance recently emerged in analog methods, which are based on nearest neighbours [START_REF] Lguensat | The Analog Data Assimilation[END_REF], in ocean and atmosphere science. Though leading to significant advances for simulation and reconstruction issues, the lack of physical interpretability of such approaches advocates for frameworks bridging the physical paradigm and the statistical paradigm underlying machine learning. In this respect, neural networks, especially Residual Networks (ResNet), are of key interest [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. More specifically, the discretized numerical integration schemes of an ODE can be stated as a ResNet [START_REF] Wang | Runge-Kutta neural network for identification of dynamical systems in high accuracy[END_REF][START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF], which allows for the data-driven identification of the dynamical operator governing a dynamical system of interest. To our knowledge, previous works mostly focused on the parameterization of the dynamical operator using or not physical priors (e.g., bilinear setting [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF], physics-informed parameterization [START_REF] De Bezenac | Deep learning for physical processes : Incorporating prior scientific knowledge[END_REF]). Besides the identification of the dynamical operator, the selection of the numerical integration scheme (e.g., explicit Eurler and Runge-Kutta schemes) may significantly affect modeling and forecasting performance, especially when the data available for training may not involve fine-scale time sampling with respect to the characteristics time scales of the considered processes, as shown in a previous work [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF]. Such issues appear critical when addressing the identification of dynamical systems from observation data (e.g., satellite earth observation data).

In this work, we address the joint data-driven identification of the dynamical operator governing a dynamical process of interest and the associated numerical integration scheme. We propose a residual integration neural network which jointly learns a dynamical model and an explicit Runge-Kutta integration scheme with an arbitrary number of stages. From an insight on high order numerical integration schemes, we demonstrate the relevance of the proposed architecture for identification and forecasting purposes when considering large integration step. We make explicit the relationship between the number of residual blocks and the order of a given integration scheme in terms of truncation error and show that increasing the number of residual layers in our architecture results in a behaviour similar to the integration of an ODE depending on the integration time-step. Overall, our key contributions are three-fold : i) we propose a new neural network architecture for the joint identification of dynamical systems and their corresponding integration scheme ; ii) we make explicit the link between the considered residual architecture and high-order integration schemes in terms of truncation error, iii) we demonstrate the relevance of the proposed architecture for Lorenz-63 and Lorenz-96, which are representative of chaotic geophysical dynamics.

NUMERICAL INTEGRATION AND RUNGE-KUTTA METHODS

This section briefly introduces Runge-Kutta numerical integration schemes, that will provide the basis for the definition and analysis of the proposed residual integration architecture.

Let us consider a dynamical system, whose time-varying state X t is governed by the following ordinary differential equation (ODE) :

dXt dt = F (t, Xt) (1)
where F is the dynamical operator. Most of the time, this ODE cannot be solved analytically and numerical integration techniques using discrete approximations are implemented.

Assuming that we are provided with an initial condition X t0 , we aim to solve the ODE for an interval t ∈ [t 0 , t f ]. Given a discretization of the interval using a time-step h > 0 as h = t f -t0 N and t n = t 0 + nh, where 0 < n < N an integer and N is the number of grid points, it comes to approximate the value of variable X t at each grid point : X t1 , ..., X tn , ..., X t N . Explicit and implicit numerical integration schemes may be considered [START_REF] Fried | Numerical Solution of Differential Equations[END_REF]. In this work, we focus on explicit integration schemes. A one-step explicit integration scheme is defined as :

Xt n+1 = Xt n + hΦ(tn, Xt n , h) (2)
with Φ(t n , X tn , h) a numerical integration operator. Here, we aim to learn a prediction operator based on Φ so that the forecasting error (typically, a one-step-ahead error) is minimized.

From an integration point of view, one may rather consider the truncation error to characterize the numerical integration scheme. The truncation error is defined with respect to the true analytic solution X T t as follows :

en = X T t n+1 -Xt n+1 (tn, X T tn , h) (3)
en = X T t n+1 -X T tn -hΦ(tn, X T tn , h) (4)
A p-order numerical resolution method can be derived, using the Taylor development of the analytic solution X T t up to the order p + 1. Assuming that F is a C p function (p times derivable with a continuous p th derivative), we can write the Taylor expansion as :

X T t n+1 = X T tn + p k=1 h k 1 k! F k-1 (tn, X T tn ) + h p+1 1 (p + 1)! F p (tn, X T tn ) + O(h p+2 ) (5)
The corresponding p-order numerical integration scheme can be derived by replacing Φ(t n , X tn , h) in Equation 2 such as :

Φ(tn, Xt n , h) = p k=1 h k-1 1 k! F k-1 (tn, Xt n ) (6)
The corresponding truncation error of the p-order method can be deduced by neglecting the term O(h p+2 ) in equation 5 as :

en = h p+1 1 (p + 1)! F p (tn, Xt n ) (7)
The explicit Euler method corresponds to p = 1 and its truncation error is proportional to h 2 . To use a first-order method like Euler, the integration time step should be small enough which is not always possible for complex systems due to computational issues. Higher-order techniques are more robust to the integration time step [START_REF] Fried | Numerical Solution of Differential Equations[END_REF]. However, the computation of high-order derivatives becomes quickly expensive which may limit their use in practice. Runge-Kutta integration schemes were introduced as an efficient trade-off between high-order approximations and computational complexity. It relies on the following recurrent update :

Xt n+1 = Xt n + s i=1 βiki ( 8
)
where s is the number of stages of the method,

k i = F (t n + c i h, X tn + h( i-1 j=1 α i,j k j )) with 0 < j < i ≤ s and s i=1 β i = 1, 0 < c i < 1, i-1 j=1 α i,j = c i . When s = 1
, it simply corresponds to the explicit Euler method. For a given number of stages s, Runge-Kutta method coefficients need to satisfy some extra conditions (by matching it to the corresponding Taylor series) to get a given order p [START_REF] Butcher | Coefficients for the study of runge-kutta integration processes[END_REF]. Formally, the Runge-Kutta method order p is always inferior or equal to the number of stages s. For s = 4, we retrieve the well-known Runge-Kutta-4 method. For p > 4, we need more integration stages s to truly reach a given error order p.

RESIDUAL INTEGRATION NETWORK

This section introduces the proposed residual integration neural network (RINN) framework. We first introduce the proposed architecture and the associated learning scheme.

We then analyze the characteristics of RINNs in terms of forecasting error.

Proposed Neural Network architecture

Let us assume we are provided with representative time series {X} with a given time sampling rate h, which are governed by an unknown ODE. For the sake of simplicity, we consider below a single time series of length N + 1, {X 0 , X 2 , . . . , X N }. The same applies for a dataset formed by different time series possibly of varying lengths. We aim to identify the unknown dynamical operator F (Equation ( 1)) from time series {X} when sampling rate h may be high. As illustrated in the reported experiments, in such situations, Euler-based learning schemes [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF] may fail in providing relevant forecasts.

Motivated by the effectiveness of high-order integration schemes in solving differential equations with relatively high time-step h, we propose a novel architecture based on residual networks and Runge-Kutta to effectively identify dynamical systems when provided with observations with low time sampling rates. The proposed architecture involves a residual neural network architecture. A residual block F N N is shared upon all the residual layers up to the predefined stage S. This residual block is the neural-network parameterization of the dynamical operator F in Equation (1). Our architecture mimics a Runge-Kutta numerical integration scheme with S stages, which imposes the following constraints on weighing parameters {β i } i , {α i,j } i and {c i } i :

ŝ i=1 βi = 1, ∀i, 0 < ci < 1 and i-1 j=1 αi,j = ci (9)
Overall, two main components need to be defined to specify a RINN :

-The parametrization chosen for the residual block F N N approximating our true dynamical model F in terms of neural network structures. It may rely on physicsinformed parameterizations [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Raissi | Multistep Neural Networks for Data-driven Discovery of Nonlinear Dynamical Systems[END_REF][START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF][START_REF] De Bezenac | Deep learning for physical processes : Incorporating prior scientific knowledge[END_REF]. ;

-The number of stages S of our residual integration network.

The learning procedure is stated as the minimization of the forecasting error subject to constrain (9) :

min θ N N ,c,β,α N n=1 X T tn -Ψ(X T t n-1 , θNN , c, β, α)
subject to(9) [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF] where Ψ is the output of the RINN obtained by applying the Runge-Kutta recursion (Equation 8) based on the approximate model F N N and the weights c, α and β introduced above. θ N N are the parameters of operator F N N . We implement the considered architecture and learning criterion under pytorch framework. As optimization solver, we consider the ADAM algorithm [START_REF] Diederik | Adam : A Method for Stochastic Optimization[END_REF]. The constrained optimization is solved by clipping the integration weights after each training epoch.

Performance of the learnt RINN

Assuming that the S-stage RINN corresponds to an porder numerical integration scheme, the loss function of the RINN relates to the truncation error of the learnt integration scheme :

ên 2 = (X T t n+1 -Xt n+1 ) 2 ( 11
)
where Xtn+1 is the output of our RINN.

Using the Taylor expansion given by equation 5 over the true state X T tn+1 up to the order p + 1, the training error of the learnt p-order numerical integration scheme is given by :

ên 2 = ( p k=1 h k 1 k! F k-1 (tn, X T tn ) + h p+1 1 (p + 1)! F p (tn, X T tn ) - p k=1 h k 1 k! F k-1 (tn, X T tn )) 2 (12)
This squared truncation error depends on two learnt parameters θ N N (i.e., the parameters of dynamical operator F N N ) and p. It reaches a minimum for F N N = F and p = p. Hence, a theoretic lower bound of the training loss function of the RINN is given by the truncation error of our true dynamical model :

ên 2 > (h p+1 1 (p + 1)! F p (tn, X T tn ))) 2 (13)
Equations 12 and 13 illustrate two main characteristics about learning neural network representations of dynamical models : (i) one cannot expect a training error lower then a theoretical lower bound represented by the truncation error of the true dynamical model (assuming p to be high enough to properly integrate the true dynamical system so that terms proportional to h p+2 or higher are negligible), (ii) we may jointly tune F N N and p in the RINN architecture to lower the training loss function. Assuming that the integration time step h is set by the temporal sampling of our training data, one may improve the approximation F N N of the true dynamical model F as mostly studied in the data-driven community [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF]. One may also decrease the training loss function through a greater order p of the integration scheme reproduced by the RINN, parameterized through the number of stages S. This clearly motivates the development of residual networks with several residual layers. Our previous work [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF] can be regarded as an illustration of this aspect. We showed that a residual dynamical model reproducing the Runge-Kutta-4 scheme outperforms the model reproducing the Euler setting and this even with a low integration time-step. Those results clearly relate to Equation 12.

NUMERICAL EXPERIMENTS

In this section, we evaluate the proposed framework and demonstrate its relevance to identify and forecast dynamical systems governed by an unknown ODE when only provided with data with a low sampling rate. As case studies, we consider two models widely studied in geophysics as examples of chaotic patterns in ocean-atmosphere science [START_REF] Edward | Deterministic Nonperiodic Flow[END_REF].

Case studies

Lorenz-63 system : The Lorenz 63 dynamical system is a 3-dimensional model governed by the following ODE :

     dX t,1 dt = σ (Xt,2 -Xt,2) dX t,2 dt = ρXt,1 -Xt,2 -Xt,1Xt,3 dX t,3 dt = Xt,1Xt,2 -βXt,3 (14)
Under parameterization σ = 10, ρ = 28 and β = 8/3, this system involves chaotic dynamics with two attractors [START_REF] Edward | Deterministic Nonperiodic Flow[END_REF]. Lorenz-96 system : The Lorenz 96 dynamical system is a 40-dimensional system. It involves propagation-like dynamics governed by :

dXt,i dt = (Xt,i+1 -Xt,i-2)Xt,i-1 + A ( 15
)
with periodic boundary conditions (i.e. X t,-1 = X t,40 and X t,41 = X t,1 ) and A = 8. Lorenz-96 system provides means to demonstrate the relevance of the proposed framework for higher-dimensional systems.

We simulate Lorenz-63 (resp. Lorenz-96) state sequences using the LOSDA ODE solver [START_REF] Hindmarsh | ODEPACK, a systematized collection of ODE solvers[END_REF] with an integration step of 0.01 (resp. 0.05). We then subsample the simulated sequences to different timesteps while making sure we don't exceed the characteristic time-scale.

Results

In this section, we compare several Residual Integration neural network performances in predicting the Lorenz dynamics from a given initial state. For benchmarking purpose, the following models were tested :

-Sparse regression model [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] (SR) : This model computes a sparse regression over an augmented states vector based on second order polynomial representations of the Lorenz states. The learnt dynamical model is then integrated to compute forecasts using the LOSDA ODE solver [START_REF] Hindmarsh | ODEPACK, a systematized collection of ODE solvers[END_REF]. -Residual Integration Neural Network 1 (RINN1) : the proposed residual architecture with a number of stages equal to one. This corresponds to the first order Euler integration method.

-Residual Integration Neural Network 4 (RINN4) : the proposed residual architecture with a number of stages equal to four. This comprises the fourth-order Runge-Kutta 4 integration technique with integration parameters {β i } i , {α i,j } i and {c i } i set to the true Runge Kutta 4 parameters.

-Residual Integration Neural Network 11 (RINN11) : Proposed residual architecture with a number of stages equal to 11. In this architecture, the weights of the integration scheme are learnt as explained in section 3.1.

In all these reported experiments, the parameterization used for the neural-network approximation F N N of the dynamical operation F is a bilinear architecture as proposed in [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF]. This bilinear architecture ensures that the true model truly lies within the space of possible model parameterizations.

We report the forecasting performances in Tab. 1 and 2. Figure 1 illustrates the trajectories generated using the trained data-driven models on the Lorenz-63 system with h = 0.4. The proposed residual integration neural network with 11 residual layers leads to the best performances comparing to identification techniques based on low order integration schemes. This gain clearly motivates the investigation of such representations in data-driven dynamical modeling.

CONCLUSION

In this work, we demonstrate the relevance of the residual integration neural network in the identification of dynamical systems. Through the representation of residual networks as high order numerical integration schemes, we prove that high order residual networks can allow the learning of dynamical models even when provided with training data with low temporal sampling. Further works could investigate a relationship between the number of stages in our residual network and the corresponding order of the learnt numerical integration scheme.

Fig. 1 :

 1 Fig.1: Generated time series of the proposed models. Generated time series of the proposed models for h = 0.4. For visualization purpose, the time series were interpolated to an h = 0.01 grid using a cubic interpolation.

Table 1 :

 1 Forecasting performance of data-driven models for Lorenz-63 dynamical model : mean RMSE for different forecasting time steps.

	Model	h=0.3	h=0.4	h=0.5
	SR	t 0 + h t 0 + 4h	11.10 12.56 9.64 12.51 57.90 7.48
	RINN1	t 0 + h t 0 + 4h	11.36 10.64 8.09 9.96	3.83 8.60
	RINN4	t 0 + h t 0 + 4h	2.24 8.33	7.66 10.79	2.80 8.64
	RINN11	t 0 + h t 0 + 4h	0.23 1.01	0.54 2.22	0.41 2.06
	Model	h=0.3 h=0.4 h=0.5
	RINN4	t 0 + h t 0 + 4h	1.30 2.69	2.89 3.31	2.76 3.05
	RINN11	t 0 + h t 0 + 4h	0.02 0.09	0.68 2.48	2.06 3.34

Table 2 :

 2 Forecasting performance of data-driven models for Lorenz-96 dynamical model : mean RMSE for different forecasting time steps.

This work was supported by GERONIMO project (ANR-13-JS03-0002), Labex Cominlabs (grant SEACS), Region Bretagne, CNES (grant OSTST-MANATEE), Microsoft (AI EU Ocean awards) and by MESR, FE-DER, Région Bretagne, Conseil Général du Finistère, Brest Métropole and Institut Mines Télécom in the framework of the VIGISAT program managed by "Groupement Bretagne Télédétection" (BreTel).


Said Ouala

Ronan Fablet

Cédric Herzet

Bertrand Chapron

Ananda Pascual

Fabrice Collard

Lucile Gaultier

LEARNING STOCHASTIC REPRESENTATIONS OF GEOPHYSICAL DYNAMICS

Keywords: Probabilistic modeling, Dynamical systems, Neural networks, Kalman filter

In the last years, Neural Networks have enriched the stateof-the-art in probabilistic modeling. This is principally due to the advances in deep learning which allow a better understanding of complex systems. However, the stochastic representation of spatio-temporal fields is still an open challenge that may benefit from the recent advances in probabilistic modelization. In this work, we explore neural network to derive a stochastic representation of spatio-temporal dynamical systems based on ensemble forecasting. Trough the implementation of our stochastic model in a classical Kalman filtering scheme, we demonstrate the relevance of the proposed architecture in the reconstruction of geophysical fields with respect to the state-of-the-art approaches.

PROBLEM STATEMENT AND RELATED WORK

The high-resolution monitoring of sea surface geophysical parameters is one of the significant challenges in oceanography. Producing high resolution gridded spatio-temporal products of physical variables such as sea surface temperature, sea surface height and sea surface salinity is of key interest for a variety of scientific fields and applications, from numerical weather prediction to maritime operations [START_REF] Hardman-Mountford | Relating sardine recruitment in the Northern Benguela to satellite-derived sea surface height using a neural network pattern recognition approach[END_REF][START_REF] Le Traon | Satellites and operational oceanography[END_REF][START_REF] Von Schuckmann | The copernicus marine environment monitoring service ocean state report[END_REF].

Data assimilation is known as the state-of-the-art techniques used for the reconstruction of such geophysical fields. They are present in two main categories : variational and stochastic data assimilation. While variational data assimilation relies on the minimization of a cost function and requires the computation of the cost function gradient through the adjoint This work was supported by GERONIMO project (ANR-13-JS03-0002), Labex Cominlabs (grant SEACS), Region Bretagne, CNES (grant OSTST-MANATEE), Microsoft (AI EU Ocean awards) and by MESR, FE-DER, Région Bretagne, Conseil Général du Finistère, Brest Métropole and Institut Mines Télécom in the framework of the VIGISAT program managed by "Groupement Bretagne Télédétection" (BreTel). of the dynamical model [START_REF] Lorenc | The Met. Office global three-dimensional variational data assimilation scheme[END_REF]. Stochastic data assimilation involves Monte Carlo strategies through an ensemble formulation to explore the state space of our physical model [START_REF] Evensen | Data Assimilation[END_REF].

The Key feature of both variational and stochastic data assimilation techniques is the mathematical formulation of the dynamical model describing the temporal evolution of our physical variables. The definition of this dynamical model dramatically limits the use of data assimilation in the case of geophysical fields interpolation. First, analytical modeling techniques may result in complex representations that can not be used in assimilation schemes due to their numerical complexity [START_REF] Hong | Next-generation numerical weather prediction : Bridging parameterization, explicit clouds, and large eddies[END_REF]. Complex dynamics may also involve uncertain parameterizations, whereas simplified parameterizations may only hold for specific dynamical modes [START_REF] Jose | Estimating model parameters with ensemble-based data assimilation : A review[END_REF].

Data-driven strategies in the other hand arise as a relevant alternative to model-driven techniques [START_REF] Tandeo | The analog data assimilation : application to 20 years of altimetric data[END_REF]. several representations have shown a great success in forecasting physical dynamics in several domains [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Tandeo | The analog data assimilation : application to 20 years of altimetric data[END_REF][START_REF] Ouala | Sea surface temperature prediction and reconstruction using patch-level neural network representations[END_REF]. However, when dealing with chaotic systems such as oceanic tracers, modeling the variability in terms of probabilistic representation is still an open challenge due to the complexity of such dynamics.

Neural networks, especially deep models, has become the state-of-the-art in a range of machine learning applications. Among others, neural networks with stochastic layers and generative models have significantly advanced the field of probabilistic modeling for complex processes and signals [START_REF] He | Deep Residual Learning for Image Recognition[END_REF][START_REF] Fraccaro | Sequential Neural Models with Stochastic Layers[END_REF].

In this study, we are specifically interested in the probabilistic modeling of spatio-temporal fields for simulation and assimilation issues. We investigate stochastic data-driven representations of spatio-temporal dynamics based on neural networks. The key idea of our work is to design a Gaussian neural-network representation for the propagation of spatiotemporal dynamics. Inspired by the variational approximation based techniques [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], the proposed approach can be regarded as a neural-network representation of the ensemble forecasting based scheme classically used in Bayesian filters (e.g., particle filters, Kalman filters, Ensemble Kalman fil-ters) to propagate probabilistic representations of non-linear dynamics [START_REF] Evensen | Data Assimilation[END_REF]. The proposed neural network representation is learnt from data. It can be combined with a classic Kalman recursion to address assimilation issues. Overall, our key contributions are three-fold : i) we propose a novel Gaussian NNbased representation of high-dimensional geophysical dynamics, ii) we derive the associated NN-based Kalman filtering scheme for spatio-temporal interpolation issues, iii) We demonstrate the relevance of these contributions with respect to state-of-the-art approaches [START_REF] Donlon | The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system[END_REF][START_REF] Ping | An Improved DINEOF Algorithm for Filling Missing Values in Spatio-Temporal Sea Surface Temperature Data[END_REF][START_REF] Fablet | Data-Driven Models for the Spatio-Temporal Interpolation of Satellite-Derived SST Fields[END_REF] for an application to the spatio-temporal interpolation of satellite-derived data.

This paper is organized as follows. Section 2 describes the proposed stochastic neural-network dynamical model and its implementation in a Kalman filtering scheme. Section 3 presents the results of the numerical experiments. We further discuss our contributions in Section 4.

LEARNING GAUSSIAN DYNAMICAL MODELS

Proposed stochastic dynamical model

We aim to model the propagation of a Gaussian approximation (i.e., a mean and covariance) of the true dynamics of the system of interest. Contrary to parametric techniques [START_REF] Pannekoucke | Modelling of local length-scale dynamics and isotropizing deformations[END_REF][START_REF] Pannekoucke | Parametric Kalman filter for chemical transport models[END_REF] based on an analytic PDE formulation of the propagation of a Gaussian dynamical model using specific dynamical priors, we aim to benefit from the flexibility and genericity of neural-network architectures to learn our Gaussian dynamical model from a representative training dataset.

Formally, assuming that we are provided with a dataset X 0:T of our physical system, our objective is to learn a probabilistic representation of the dynamical model through a parametric model of the conditional distribution of state variable X t given X t-1 . Considering a Gaussian approximation, it comes to the following model :

Xt|Xt-1 ∼ N (Xt = F(Xt-1), Σt = FΣ(Xt-1, Σt)) (1)
where functions F, F Σ are neural networks with parameter vectors θ = (θ µ , θ Σ ).

Following our previous works on analog data assimilation [START_REF] Lguensat | Data-driven Interpolation of Sea Level Anomalies using Analog Data Assimilation[END_REF][START_REF] Fablet | Data-Driven Models for the Spatio-Temporal Interpolation of Satellite-Derived SST Fields[END_REF], we consider a patch-based representation as a mean to decompose space-time scales as well as to provide an explicit relationship between the global and local (patch-level) representations. Regarding the mean model F, the proposed architecture proceeds as follows :

-At a given time t, the first layer of the network, which is parameter-free in terms of training, comes to decompose an input field x t into a collection of N p P × P patches x Ps,t , where P is the width and height of each patch and s ∈ [1, ..., N p ] the patch location in the global field. Each patch is decomposed onto an EOF basis B according to :

z Ps,t = x Ps,t B T (2)
with z Ps,t the EOF decomposition of the patch x Ps,t . The EOF decomposition matrix B is trained offline as a preprocessing step ; -The second layer implements a numerical integration scheme (typically, an Euler or 4th-order Runge-Kutta scheme) using a patch-level dynamical model F Ps to predict z Ps,t+1 . For patch-level models F Ps , we consider residual architectures [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] with a bilinear parameterization [START_REF] Fablet | Bilinear residual Neural Network for the identification and forecasting of dynamical systems[END_REF] ; -The third layer is a reconstruction network F r . It combines the predicted patches x Ps,t+1 = z Ps,t+1 B to reconstruct the output field x t+1 . This reconstruction network F r involves a convolution neural network [START_REF] Lecun | Object Recognition with Gradient-Based Learning[END_REF]. The details of the considered parameterizations for the second and third layers are given in Section 3. To train the mean dynamical model F, we apply a two-step procedure. We first learn the local dynamical models F Ps , s ∈ [1, ..., N p ] based on the minimization of the EOF-patch based forecasting error. The reconstruction network F r is then optimized using the same criterion over the global field.

Due to dimensionality issues, one cannot expect to model the full covariance of the field. We benefit from the patchlevel representation to investigate a block-diagonal approximation of the covariance structure. It comes to consider a diagonal covariance parameterization of the patch-level covariance matrix in the EOF space. Formally, the block-diagonal covariance model is learnt as follows :

-For the t th state X t in our training set, a patch-based projection is first performed to compute the states in the EOF space z Ps,t according to equation 2. -For each s ∈ [1, ..., N p ], an ensemble E Ps,t = {z

Ps,t , ..., z

(M )
Ps,t } of size M is sampled from a Gaussian distribution with mean z Ps,t and covariance matrix Σ 0 computed as the covariance matrix of our training data in the EOF space.

-From the ensemble data E Ps,t , we build a forecasting ensemble F Ps (E Ps,t ) by applying to each sample the trained model F Ps . A neural network covariance model F Ps D,Σ is then trained to forecast a diagonal covariance matrix in the EOF space Σ EOF Ps,t+1 = F Ps D,Σ (z Ps,t , Σ 0 ) using the ensemble dataset {E Ps,t-1 , F Ps (E Ps,t-1 )} according to a ML criterion.

-Ones our EOF covariance model learned, the patchbased covariance model F Ps Σ is deduced according to :

F Ps Σ (x Ps,t+1 , Σ Ps,t+1 ) =B t Ψ(Σ Ps,t , Σ 0 ) × F Ps D,Σ (z Ps,t , Σ 0 ) • B (3)
where Ψ a scaling function.

-The global covariance model F Σ is given by the recombination of the N p local patch-level covariance models F Ps Σ .

The parametrization of the patch-based covariance model in the EOF-space allows us learning only the variance vector of each patch sins the EOF components are supposed to be orthogonal (resulting in insignificant covariance components in the EOF-space). However, and to illustrate the relevance of the proposed patch-based parametrization of the covariance matrix in the EOF-space, we also investigate a diagonal covariance matrix model in the patch space. according to :

F Ps Σ (xP s ,t+1, ΣP s ,t+1) = Ψ(ΣP s,t , Σ0)•F Ps D,,Σ (xP s,t , Σ0) (4)
The aim of using an ensemble likelihood optimization criteria is that geophysical dynamical models usually present complex, possibly chaotic, dynamical behaviours. Therefore, using a single forecast to characterize the variability of the dynamical model can result in an underestimation (or an overestimation) of the covariance matrix. The proposed Monte Carlo ML criterion explores the one-step prediction variability resulting in a parametric Gaussian approximation of an ensemble forecasting.

Application to data assimilation

Given this neural-network representation of the conditional distribution X t |X t-1 , we can derive an associated Kalman-based filtering under the assumption that the observation model is linear and Gaussian.

In contrast to classical non linear filtering frameworks illustrated for instance by Extended Kalman filtering and Ensemble Kalman filtering. The proposed Kalman-based filtering technique directly propagates the Gaussian pdf (see equation 5) at a given time step (using the Gaussian transition model). It simply comes to apply the Kalman recursion with the mean and covariance propagated by the trained neural networks :

(X - t+1 , Σ - t+1 ) = (F (X + t ), F Σ (X + t , Σ + t )) (5)
X + t+1 = X - t+1 + K t+1 [Y t+1 -H t+1 X - t+1 ] (6)
K t+1 = Σ - t+1 H t t+1 [H t+1 Σ - t+1 H t t+1 + Rt] -1 (7)
where H t+1 corresponds to some linear observation models. The superscript (-) refers to the forecasting of the mean of the state variable x - t+1 and of its covariance matrix Σ - t+1 , the superscript (+) refers in the other hand to the assimilated mean x + t+1 and covariance Σ + t+1 . R t is a Gaussian noise covariance matrix.

As explained above, we consider two different parameterizations of the covariance model F Σ . The key aspect of this model is the structure of the patch based covariance operator F Ps Σ . Different assimilation techniques could be envisaged depending on its parametrization :

-Modeling the patch based covariance in the EOF space and projecting it back to the patch space which results in a block-diagonal covariance matrix (illustrated by equation 3) : PB-NNKF-EOF ;

-Directly modeling the variance vector of each patch resulting on a diagonal covariance matrix (illustrated by equation 4) : PB-NNKF.

NUMERICAL EXPERIMENTS

Dataset description

As a case-study, we address the spatio-temporal interpolation of satellite-derived sea surface temperature (SST) fields associated with infrared sensors. The time series is delivered by the UK Met Office [START_REF] Donlon | The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system[END_REF] from January 2008 to December 2015. The spatial resolution of our SST field is 0.05 • and the temporal resolution h = 1 day. The data from 2008 to 2014 were used as training data and we tested our approach on the 2015 data. To perform a quantitative evaluation, we simulated realistic spatio-temporal cloud patterns over the test set using the METOP-AVHRR masks. This sensor is highly sensitive to the cloud cover and results in very high missing data rates.

As case-study region, we select an area off South Africa (from 2.5 • E, 38.75 • S to 32.5 • E, 58.75 • S). This region involves complex fine-scale SST dynamics (e.g., fronts, filaments). It makes it relevant for the considered quantitative evaluation.

Experimental setting

The proposed stochastic dynamical model exploits patchlevel representations with non-overlapping 20 × 20 patches. Each patch P s is then projected into an EOF basis learnt from the training data. We keep the first 50 EOF components, which amount to encode 95% of the total variance. For the patch-level NN model F Ps , we use a bilinear residual neural network architecture as proposed in [START_REF] Ouala | Sea surface temperature prediction and reconstruction using patch-level neural network representations[END_REF] with 60 linear neurons, 100 bilinear neurons and 10 fully-connected layers with a Relu activation. The reconstruction model F r is a convolutional neural network with 3 convolutional layers. The first two layers comprise 64 filters of size 3 × 3 with a Relu activation and the last layer is a linear convolutional layer with one filter. Regarding covariance model F Ps D , we consider a diagonal covariance model within each patch. Each element of diagonal involves a 3-layer MLP with 4 neurons and Relu activation functions on the hidden layers and a softplus activation in the output layer. With a view to evaluating the EOF-based covariance parameterization, we consider both PB-NNKF-EOF and PB-NNKF schemes. A constant scaling function Ψ() = 1 in equations 3 and 4 led to the best performances in our experiments.

We perform a quantitative analysis of the interpolation performance of the proposed schemes with respect to an optimal interpolation, the analog data assimilation [START_REF] Fablet | Data-Driven Models for the Spatio-Temporal Interpolation of Satellite-Derived SST Fields[END_REF] and the EOF based interpolation method VE-DINEOF. The considered parameter setting is as follows :

-Optimal interpolation (OI) : We use a Gaussian kernel with a spatial correlation length of 100km and a temporal resolution length of 3 days. These parameters were empirically tuned for the considered dataset using a cross-validation experiment. -Analog data assimilation (LAF-EnKF) : We test the local analog data assimilation scheme [START_REF] Fablet | Data-Driven Models for the Spatio-Temporal Interpolation of Satellite-Derived SST Fields[END_REF][START_REF] Lguensat | The Analog Data Assimilation[END_REF]. Similarly to the proposed architecture, we consider 20×20 patches and 50-dimensional EOF decomposition with an overlapping of 10 pixels. We let the reader refer to [START_REF] Fablet | Data-Driven Models for the Spatio-Temporal Interpolation of Satellite-Derived SST Fields[END_REF][START_REF] Lguensat | The Analog Data Assimilation[END_REF] for a detailed description of this data-driven approach, which relies on nearest-neighbor regression techniques. -EOF based reconstruction (PB-VE-DINEOF) : We also compare our approach to the state-of-the-art interpolation scheme based on the projection of our observations with missing data on an EOF basis [START_REF] Ping | An Improved DINEOF Algorithm for Filling Missing Values in Spatio-Temporal Sea Surface Temperature Data[END_REF]. The SST field is here decomposed as described in the analog data assimilation application into a collection of 20 × 20 patches with 10 pixels overlapping. Each patch is then reconstructed using the VE-DINEOF method.

Results

We report the mean interpolation performance in Tab. 1. Figure 1 illustrates an interpolation example of the tested methods. The proposed NN-based scheme (PB-NNKF-EOF) leads to very significant improvements with respect to the optimal interpolation in terms of RMSE and correlation coefficients, which emphasizes fine-scale structures (e.g., relative improvement of the RMSE above 50% for missing data areas). A clear gain is also exhibited w.r.t. analog data assimilation and PB-VE-DINEOF schemes with a relative gain greater than 20% in terms of RMSE. Although the considered NN-based representation exploits non-overlapping patches, we still come up with significant improvements w.r. between patches. This clearly illustrates the relevance of NN-based representation, which fully embeds the direct and inverse mappings between the SST field and its patch-level representation. Interestingly, Tab.1 also reveals the importance of the EOF-based parameterization of the NN-based covariance model (Equation 3) in the improvement of interpolation results w.r.t. AnDA schemes.

CONCLUSION

In this work, we address stochastic data driven representations for spatio-temporal fields identification. we rely on neural networks and ensemble forecasting to derive a neural network Gaussian dynamical model. From our numerical experiments, combining our probabilistic representation with Kalman based data assimilation techniques outperforms classical state-of-the-art techniques in spatio-temporal fields reconstruction.

Fig. 1 :

 1 Fig. 1: Interpolation of the SST field on July 19 2015 : first row, the reference SST and the observation with missing data (here, 82% of missing data) ; second row, interpolation results using respectively OI, PB-VE-DINEOF, LAN-ENKF, PB-NN-NNKF, PB-NN-NNKF-EOF.

Table 1 :

 1 t AnDA scheme which involve a 50% overlapping rate SST interpolation experiment : Mean reconstruction correlation coefficient and RMSE over the SST time series.

	Model		Entire map Missing data areas
	OI	RMSE Correlation	0.76 99.35%	0.75 99.37%
	PB-VE-DINEOF	RMSE Correlation	0.54 99.68%	0.54 99.66%
	LAF-EnKF	RMSE Correlation	0.43 99.79%	0.42 99.77%
	PB-NNKF	RMSE Correlation	0.51 99.75%	0.51 99.71%
	PB-NNKF-EOF	RMSE Correlation	0.33 99.87%	0.35 99.85%


J Lerbet
email: jean.lerbet@ibisc.univ-evry.fr

M Aldowaji

N Challamel

F Nicot

F Prunier

F Darve

P-positive de nite matrices and stability of nonconservative systems

Keywords: stability, kinematic constraints, second-order work criterion

Key words p-positive definite matrices, nonconservative systems

Introduction

This paper investigates the linear static stability of constrained nonconservative mechanical systems. More precisely, the systems studied are elastic systems subjected to nonconservative positional forces (or circulatory forces) and we investigate the loss of stability by divergence but the similar mathematical problem is met in geomechanics for non associated plastic materials. It is well known that for such systems, there is no general inclusion between the domain of instability by divergence and the domain of instability by flutter (see for example [START_REF] Kirillov | Collapse of the Keldysh chains and stability of continuous non-conservative systems[END_REF][START_REF] Seyranian | Bifurcation diagrams and stability boundaries of circulatory systems[END_REF] for precisions). It is also well known that such systems may present paradoxical behaviors if some damping is introduced in the models: damping may increase the instability domain (see [START_REF] Bolotin | Non-Conservative Problems of the Theory of Elastic Stability[END_REF][START_REF] Ziegler | Die Stabilitätskriterien der Elastomechanik[END_REF] for example and most recently see also [START_REF] Kirillov | Paradoxes of dissipation-induced destabilization or who opened Withney's umbrella?[END_REF]). It is, however, less reported that other paradoxical effects may be met for additional constraints. This problem has recently been investigated (see [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF][START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF][START_REF] Nicot | Bifurcation and generalized mixed loading conditions in geomaterials[END_REF] for example) and remarkable new results obtained for one additional kinematic constraint: the additional constraint may destabilize the system and preventing the instability by divergence of the constrained system (ie for any kinematic constraint) leads to the second order work criterion: the symmetric part of the stiffness matrix must be positive definite. In the present paper, the general mathematical framework of the problem is highlighted, which leads to the new linear algebra concept of p-positive definite matrices (0 < p ≤ n) of size n and the previous result is generalized for any family of p independent kinematic constraints. The paper is organized as follows. In the first part, as introduction to later developments, we start by illustrating on the example of a three degree of freedom Ziegler's column how the known approach is performing and how it does not lead to appropriate results. In the second part, the general mechanical background is given, which leads to the central mathematical problem. In the third part, the new concept is defined with its elementary properties and the main result of the paper is proved. Various other interesting mathematical problems related to the new concept are also suggested. In the last part, we return to mechanics and apply the results to the problem of stability. The initial example of a 3 d.o.f. Ziegler pendulum is again used as an illustration of general results.

A motivating example

Let us consider the following 3 d.o.f. Ziegler column as in Fig. 1. The mechanical system Σ consists of three bars OA, AB, BC with OA = AB = AC = linked by three elastic springs of the same stiffness k. The nonconservative external action (the circulatory force) is the follower force P . The kinetic energy T of Σ reads (suppose an uniform distribution of mass):

T = m 2 2 7 3 θ2 1 + 4 3 θ2 2 + 1 3 θ2 3 + 3 θ1 θ2 cos(θ 1 -θ 2 ) + θ1 θ3 cos(θ 1 -θ 3 ) + θ2 θ3 cos(θ 2 -θ 3 ) = 1 2 θT M (θ) θ with θ = ⎛ ⎜ ⎝ θ 1 θ 2 θ 3 ⎞ ⎟ ⎠ , M(θ) = m 2 ⎛ ⎜ ⎝ 7 3 3 2 (cos θ 2 -θ 1 ) 1 2 (cos θ 3 -θ 1 ) 3 2 (cos θ 2 -θ 1 ) 4 3 1 2 (cos θ 2 -θ 3 ) 1 2 (cos θ 3 -θ 1 ) 1 2 (cos θ 3 -θ 1 ) 1 3 ⎞ ⎟ ⎠
the elastic energy of the springs is

U = k 2 ((θ 2 1 + (θ 1 -θ 2 ) 2 + (θ 2 -θ 3 ) 2 ) = k 2 (2θ 2 1 + 2θ 2 1 + θ 2 3 -2θ 1 θ 2 -2θ 2 θ 3 )
and the virtual power of P in any configuration θ = (θ 1 , θ 2 , θ 3 ) reads (P > 0 in compression): Dimensions of the parameters m, , k are independent. Hence one can introduce units of measurement in such a way that all these parameters would be equal to one. Let us do it. Then it remains the only dimensionless parameter β = P k which is thet loading parameter. The chosen equilibrium position is the vertical one θ eq = 0 = (0, 0, 0). The stiffness matrix reads:

P * P = P (sin(θ 3 -θ 1 )θ * 1 + sin(θ 3 -θ 2 )θ * 2 ).
K(β) = ⎛ ⎜ ⎝ 2 -β -1 β -1 2 -β -1 + β 0 -1 1 ⎞ ⎟ ⎠ .
As it is well known but nevertheless paradoxical, det(K(β)) = 1 is independent of β: the system cannot be unstable by divergence for any value of the loading parameter β (or with the notation of the paper β s = ∞): to investigate the stability of the system, it is necessary to take into account the inertia of the system and to make a spectral analysis of the classical dynamical system (see (2) below). As consequence, flutter is here the unique mode of linear instability. With the chosen mass distribution, flutter occurs for the lowest value β f such that the spectral equation P (x) = det(M (0)x + K(β)) = 0 with x = μ 2 has a multiple root. For the chosen mass distribution (uniform)

P (x) = (13/108)x 3 + (131/36 -(19/18)β)x 2 + (12 -(41/6)β 2 )x + 1.
This equation possesses a multiple root for β f ≈ 1.483 such that the discriminant

Δ = (1/46656)(45548392 -89906052β + 70924849β 2 -29092268β 3 + 6603748β 4 -790848β 5 + 39168β 6 )
vanishes.

However, as it is also known but still paradoxical, adding a kinematical constraint may destabilize the system. For example, consider the constraint θ 3 = 0. For this new two dof system Σ 1 , the reduced stiffness matrix reads: (3β) which leads to the value β = 1 for the instability by divergence of the constrained system.

K 1 (β) = 2 -β -1 -1 2 -β and det(K 1 (β)) = (1 -β)
This simple example leads to the following six questions:

• What about the (in)stability by divergence for another single kinematic constraint?

• What about the (in)stability by divergence for any other single kinematic constraint?

• Is it possible to preview the result with only the matrix K of the system without constraint ie without calculating the reduced stiffness matrices and solving the corresponding problems of divergence instability? • What about the (in)stability by divergence for two kinematic constraints? • What about the (in)stability by divergence for any system of two kinematic constraints?

• Is it possible to preview the result with only the matrix K of the system without constraint ie without calculating the reduced stiffness matrices and solving the corresponding problem of divergence instability?

In [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF][START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF][START_REF] Nicot | Bifurcation and generalized mixed loading conditions in geomaterials[END_REF], the answer of the three first questions leads to the beautiful following result. Let Σ be any n dof mechanical system the stiffness matrix of which is K(β) not necessary symmetric. Suppose Σ divergence stable for the range of variations of β. As long as the symmetric part K s (β) of K(β) is positive definite, no single kinematic constraint may destabilize Σ. For λ m the lowest positive singular value of K s (β) (root of det(K s (β)) = 0), there exists a kinematic constraint destabilizing the system and we may find this constraint on the cone of the quadratic form x → t xK s (β m )x. The main goal of this paper is to generalize this result to any family of p kinematic constraints (0 < p < n). In the last part, we will come back to this example and we will apply to it these results for solving the last three previous questions.

Structural mechanics

General framework

Let us come back now to more general considerations. Stability of mechanical systems (body systems, continuum, discretized continuum) is considered to be a significant issue and one of the fundamental ones in mechanics and science in general. Therefore, in mechanics today, two main approaches may be identified. The first one called the "static approach" in the present paper concerns evolutions in which the inertia is not taken into account, either because it is negligible or because it does not contribute to the stability of the system. The second one can be linked to the stability in the sense of Lyapunov and mainly concerns a system's dynamic evolution. Mechanically, this implies that inertial effects are considered, this called here the "dynamic" or spectral approach.

However, it is clear that the relation between the two approaches is essential for a proper understanding of the phenomena concerned. More specifically, we assume that the studied system Σ is evolving under the action of a parameter β (loading parameter) and we study the stability of an equilibrium state of Σ, a state defined by the position of each particle and the external and internal forces allowing this equilibrium to be reached. Each approach ("static" or "dynamic") in general leads to two critical different values, β s and β d of β. These two values do not highlight exactly the same physical phenomena. The "static" approach makes it possible to understand the bifurcations of equilibrium of Σ and the "dynamic" approach the instability of the equilibrium considered (evolution of the system after a small kinematical (positions and/or velocities) disturbance of the equilibrium), the remarkable property consisting in the identity of these two phenomena for conservative systems.

Let us now consider a mechanical system Σ and an equilibrium reached under the loading β. Let M and K(β) be, respectively, the (symmetric definite positive) matrix of mass and the matrix of stiffness of Σ at equilibrium. The framework is linear mechanics and the equation of bifurcation or of "static" instability (or also of the so-called divergence) is:

det(K(β)) = 0. (1)
The (linear) "dynamic" stability analysis is traditionally done through the spectral analysis of the dynamical system:

M Ẍ + K(β)X = 0. (2)
Equation ( 1) leads to a value β s of static critical load by divergence (and sometimes it may lead to any critical load; see [START_REF] Bolotin | Non-Conservative Problems of the Theory of Elastic Stability[END_REF] for example or the previous motivating example) and ( 2) leads to a value β d of the dynamic critical load. As previously noted, if the system is conservative, then β s = β d , but for elastic nonconservative systems considering β f as the flutter critical load so that β d = Min{β s , β f }, both cases β s < β f or β s > β f may occur. Such nonconservative systems -as first reported by Ziegler [START_REF] Ziegler | Die Stabilitätskriterien der Elastomechanik[END_REF] or later by Bolotin (see [START_REF] Bolotin | Non-Conservative Problems of the Theory of Elastic Stability[END_REF] for example) -possess "counterintuitive" behaviors. For example, the introduction of damping in such systems may be unfavorable for the systems stability (see [START_REF] Bottema | The Routh-Hurwitz condition for the biquadratic equation[END_REF] for the first explanation). Another remarkable point, but less reported, concerns additional kinematic constraints to the system. While, for Rayleigh, it is conventionally accepted that the addition of such constraints reduced the extent of the spectrum and therefore supports to the stabilization (when the equilibrium is the minimum of the potential), some authors such as Tarnai ( [12], for example) have called this principle into question even in the conservative case. However, in the given examples, the addition of kinematical constraints also changes the equilibrium position itself and are not within the Rayleigh theory field of applications. On the other hand, as it is showed in the motivating example of Sect. 2, the constrained system may be less stable than the unconstrained system, this counterintuitive result illustrating the nonconservativity of Σ as well.

The goal of this paper is to highlight a number of original mathematical developments of this problem of constrained nonconservative systems.

Statement of the problem

By [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF], we consider the study of stability of a configuration of an elastic structure under any load, that is to say that the linearized equations around the equilibrium configuration studied are [START_REF] Bottema | The Routh-Hurwitz condition for the biquadratic equation[END_REF]. The linearized equations of statics become:

K(β)X = 0 (3)
a solution X = 0 indicating a "static" instability by loss of uniqueness.

It is then a solution of the equation det(K(β)) = 0 and this loss of stability is called divergence. If the system is also subjected to nonconservative actions, K(β) is any (nonsymmetric) matrix and it will be assumed in the following that, during the evolution of the system (variations of β), this matrix remains invertible: instability by divergence of the unconstrained system is assumed to be impossible. Now suppose the same system is subjected to additional kinematical constraints such that the equilibrium configuration that we may consider as 0 is not disturbed. If there are p independent constraints f (q

) = f (q 1 , . . . , q n ) = (f 1 (q), . . . , f p (q)) = 0 I R p and if A ∈ M np (I R) means the Jacobian matrix of f (A i,j = ∂fj (0)
∂qi ), then the column vector X of infinitesimal displacement must satisfy the following system:

K(β)X + AΛ = 0, t AX = 0, (4)
where Λ = t (λ 1 , . . . , λ p ) is the Lagrange multipliers associated with the constraints in the equilibrium's configuration. Following [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF], we can rewrite the previous system in a condensed manner. Let

K A (β) = K(β) A t A 0 , Y = X Λ .
So ( 4) is equivalent to:

K A (β)Y = 0. (5)
Note that:

K(β) A t A 0 K -1 (β) -K -1 (β)A 0 I pp = I nn 0 t AK -1 (β) -t AK -1 (β)A .
Thus:

det(K A (β)) det(K -1 (β)) = (-1) p det( t AK -1 (β)A), so: det(K A (β)) = (-1) p det(K(β)) det( t AK -1 (β)A). (6)
Let us write à = K -1 (β)A; it follows that:

det(K A (β)) = (-1) p det(K(β)) det( t ÃK(β) Ã) (7)
because:

det( t AK -1 (β)A) = det( t A t K -1 (β) t K(β)K -1 (β)A) = det( t Ãt K(β) Ã) = det( t ( t ÃK(β) Ã)) = det( t ÃK(β) Ã).
A part of these transformations are actually the Schur Complement formula and for more precisions see [START_REF] Zhang | The Schur Complement and Its Applications[END_REF]. The analysis of singularities of the matrix K A (β) highlighting the instability by divergence of the constrained system is brought back to the singularities of the matrix t ÃK(β) Ã, because A → Ã = K -1 (β)A is a bijection in the set of matrices M np (I R) of maximal rank p. We denote G np this set.

In the next section, we introduce the new concept of p-positive definite matrices.

p-positive de nite matrices 4.1 De nition and immediate properties: symmetric case

Although the concept of a definite positive matrix is often used for symmetric matrices, it may be defined for any square matrix:

De nition 1 Let K ∈ M nn be a square matrix, K is said to be positive de nite if and only if

t XKX > 0 ∀X ∈ M n1 X = 0
and that criterion concerns only the symmetric part

K s = K+ t K 2 of K because t XKX = t XK s X.
Any positive definite matrix K is invertible (with a strictly positive real determinant), and

K -1 is positive definite also. Moreover, let us recall that if K is symmetric positive definite, it is diagonalizable in an orthonormal basis (Y k ) k=1,...,n of I R n for the canonical scalar product I R n : t Y k Y l = δ k,l , and K defines a scalar product of I R n .
We now define a sort of generalization of the useful concept of a definite positive matrix.

De nition 2

Let K ∈ M n (I R) be a square matrix, 1 ≤ p ≤ n an integer. K is said to be p-positive de nite if, for all matrices A ∈ M np of rank p, then det( t AKA) > 0. The subset of M n (I R) constituted by p-positive de nite matrices is noted D p,n .

For p = 1 it is obvious that 1-positive definite is equivalent to positive definite.

Remark 1 Let φ K (A) = det( t AKA) ∀A ∈ M np , 1 ≤ p ≤ n.
Thus we have:

φ -K (A) = det( t A(-K)A) = (-1) p φ K (A). (8)
Remark 2 This de nition may be extended to p = 0 by K is 0-de nite positive if det(K) > 0 and the inclusion D 0,n ⊂ D 1,n is another way to formulate the well-known Bendixon theorem (sometimes also called the Bromwich theorem), which states that the smallest eigenvalue of the symmetric part K s of any square matrix K is lower than any real part of the eigenvalues of K ( [START_REF] Ishaq | Sur les spectres des matrices, Séminaire Dubreil[END_REF] for example).

Central result: 1-positive de nite implies p-positive de nite

We begin now by demonstrating a partial result:

Proposition 1: if K is (1-)positive de nite then K is 2-positive de nite. In other words, D 1,n ⊂ D 2,n . P r o o f. Let A ∈ M n2 of rank 2, A = (x 1 x 2 ) with x i ∈ M n1 for i = 1, 2.
We will sometimes identify M n1 and I R n .

det t AKA = det t x 1 Kx 1 t x 1 Kx 2 t x 2 Kx 1 t x 2 Kx 2 = det t x 1 K s x 1 t x 1 Kx 2 t x 2 Kx 1 t x 2 K s x 2 ,
where K = K s + K a with K s (resp. K a ) the symmetric part (resp. skew-symmetric part) of K.

det t AKA = t x 1 K s x 1 t x 2 K s x 2 -t x 1 Kx 2 t x 2 Kx 1 = t x 1 K s x 1 t x 2 K s x 2 -t x 1 K s x 2 t x 2 K s x 1 -t x 1 K s x 2 t x 2 K a x 1 -t x 1 K a x 2 t x 2 K s x 1 -t x 1 K a x 2 t x 2 K a x 1 we have t K s = K s , t K a = -K a ,
so by a straightforward calculation, we obtain:

det t AKA = det t AK s A + t x 1 K a x 2 2 . ( 9
)
If we now consider the quadratic form associated to K s , it defines a scalar product for which the Cauchy-Schwarz inequality gives:

det t AK s A = t x 1 K s x 1 t x 2 K s x 2 -t x 1 K s x 2 t x 2 K s x 1 > 0,
where the inequality is strict, because A being of rank 2, the vectors x 1 , x 2 are not collinear. We deduce that det

( t AK s A) = t x 1 K s x 1 t x 2 K s x 2 -t x 1 K s x 2 t x 2 K s x 1 > 0,
and therefore according to (9), φ K (A) > 0 and K is 2-positive definite. This proposition is actually general and constitutes the main result of the paper: every 1-definite positive matrix is p-definite positive for any 1 ≤ p ≤ n. The demonstration of this general result is straightforward if K is symmetric: Proposition 2: if K is symmetric and 1-positive de nite, then for any 1 ≤ p < n, K is p-positive de nite. P r o o f. Let us suppose that K is symmetric (1-)definite positive and let A ∈ M np be a matrix of rank p. Let us consider S the square root of K. S is also symmetric (1-)definite positive and t SS = S 2 = K. Thus:

φ K (A) = det( t AKA) = det( t A t SSA) = det( t RR),
where R = t (AS)AS ∈ M pp is a symmetric positive definite matrix, because if u ∈ M p1 we have t uRu = t u t SSu = t (Su)(Su) = Su 2 ≥ 0 and t uRu = 0, implies Su = 0. It follows that u = 0 because rank(S)= p. We then deduce that K is p-positive definite.

Let us recall that for nonconservative mechanical systems, we are concerned by nonsymmetric matrices. The general theorem holds for any matrix K: Theorem 1 If K is 1-positive de nite then for any 1 ≤ p < n, K is p-positive de nite. In other words, D 1,n ⊂ D p,n . P r o o f. Let K be positive definite. Then K s is symmetric positive definite, so there is a symmetric matrix (the square root of

K s ) denoted M ∈ M nn such that K s = M 2 .
Let A ∈ M np of rank p, then the following equalities hold:

det t AKA = det t A (K s + K a ) A = det t A M 2 + K a A = det t AM I n + M -1 K a M -1 M A = det t A t M I n + M -1 K a M -1 M A = det t AM I n + M -1 K a M -1 M A = det t (AM ) I n + M -1 K a M -1 (M A) Let U np = B ∈ M np | rank (B) = p}, then ψ : A -→ B = ψ(A) = M A is a bijection from U np into itself, we have A = M -1 B with B ∈ U np .
Then:

det t AKA = det t B I n + M -1 K a M -1 B = det t B I n + Ka B = χ B, Ka (-1),
where Ka = M -1 K a M -1 is also skew-symmetric as K a because M and M -1 are symmetric, and where

χ B, Ka (x) = det t B Ka -xI n B
could be called a generalized characteristic polynomial of Ka relative to B. λ is a root (possibly complex) of χ B, Ka if and only if the matrix t B Ka -λI n B is not invertible, this means if there is U ∈ M np (C) nonzero such that:

t B Ka -λI n BU = 0, then:

t B Ka BU = λ t BBU.
Since the matrices are real and the matrix Ka skew-symmetric, by taking the conjugate transpose (or adjoint), this gives:

-t Ū t B Ka B = λt Ū t BB.
By multiplying on the right by U :

-t Ū t B Ka BU = λt Ū t BBU so:
λ t Ū t BBU = λt Ū t BBU then:

-λ BU 2 = λ BU 2 .
Consequently, λ = -λ, meaning that the no null roots are purely imaginary and the only real root possible of χ B, K is then 0.

The roots of χ B, K are complex conjugates, whereas the polynomial is real. So let us group the conjugate roots together, let us separate the null root assumed to be of multiplicity r and the other conjugate roots assumed to be of number s and written as

λ k = iα k with α k > 0. Since the dominant coefficient of χ B, K is (-1) p det( t BB): χ B, Ka (x) = (-1) p det( t BB)x r s k=1 x 2 + α 2 k (10)
and then:

χ B, Ka (-1) = (-1) r+p det( t BB) s k=1 1 + α 2 k . ( 11
)
We can deduce that the sign of χ B, Ka (-1) is that of (-1) r+p det( t BB). But r = p -2s, so (-1) r+p = (-1) 2p-2s = 1 > 0, and det( t BB) > 0 because t BB is a symmetric real matrix, so it is diagonalizable and has real eigenvalues μ i for i = 1, . . . , p.

Let X i , i = 1, . . . , p be the eigenvectors corresponding to μ i , then

t BB X i = μ i X i , thus t X t i BB X i = μ i t X i X i ⇒ BX i 2 = μ i X i 2 ⇒ μ i > 0, ∀i = 1, . . . , p.
We deduce that χ B, Ka (-1) > 0, i.e., φ K (A) = det t AKA > 0.

The reciprocal result

First of all, let us note that the converse of Theorem 1 is false. More precisely, the converse of Proposition 1 is itself false, which implies that the converse of Theorem 1 is false: Proposition 3: For any n-square matrix K, K 2de nite positive does not imply K (1)-de nite positive. In other words,

D 2,n ⊂ D 1,n . P r o o f. Let K = -I n . -K = I n is 1-definite positive. So by Proposition 2, -K is 2-definite positive. But according to relation (8) of remark 1, φ -K (A) = φ K (A) for any A ∈ M n2 .
Thus K is 2-definite positive as well and it is not (1)-definite positive: K is actually definite negative! It is interesting to highlight the reasons why any (invertible) matrix is not a (p)-definite positive matrix. Indeed, returning to the problem of structural stability, a value of loading parameter λ corresponding for an (invertible) stiffness matrix K(λ) to pass from (1)-positive definite to not (1)-positive definite corresponds exactly to a loss of stability by divergence of an adequate constrained system. This very important result has already been proved in [START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF] and the object of the following section is to investigate the general case for any 1 ≤ p < n. We now investigate the reciprocal result of Theorem 1 by supposing that K is not definite positive but with K s invertible and by investigating the resulting properties about φ K . The following result holds:

Theorem 2 If K is any invertible (neither positive de nite nor negative de nite) matrix, then φ K (A) = det( t AKA) is not strictly positive for all A ∈ M np for 1 ≤ p < n and in particular for p = 2. P r o o f. To rank K S = r = n, the signature of the quadratic form q is (k, nk) and none of these numbers is zero. Let B = (e 1 , . . . , e k , e k+1 , . . . , e n ) be an orthogonal basis for q. We can also assume that:

q(e i ) = 1 if i = 1, . . . , k q(e i ) = -1 if i = k + 1, . . . , n.
Let F + = V ec(e 1 , . . . , e k ) and F -= V ec(e k+1 , . . . , e n ). So the sum I R n = F + ⊕F -is orthogonal for q, i.e, t e i K s e j = 0 if i = j.

Let us now consider P is the matrix passage from the canonical basis to the basis B.

t P K s P = I k 0 0 -I n-k .
Moreover, let us consider x ∈ I R n , and let us write x = x + + x -the only decomposition in I R n , so t (x + )K s x -= 0.

Let us assume now that p = 2 and let A = t (x 1 x 2 ) ∈ M n2 be a matrix of rank 2. Let us also choose t x 1 and t x 2 in I R n (identification I R n and M 1n ) orthogonal for q or K s . Then:

• t x 1 K s x 2 = 0, • t x 1 K a x 1 = 0, t x 2 K a x 2 = 0.
Therefore:

φ K (A) = det q(x 1 ) t x 1 K a x 2 t x 2 K a x 1 q(x 2 ) = q(x 1 )q(x 2 ) + ( t x 1 K a x 2 ) 2 . Thus, φ K (A) ≤ 0 ⇐⇒ ( t x 1 K a x 2 ) 2 ≤ -q(x 1 )q(x 2 ). By homogeneity (x i → xi q(xi) ), we can choose q(x 1 ) = 1, q(x 2 ) = -1. It follows that φ K (A) = -1 + ( t x 1 K a x 2 ) 2 .
1. If K a is singular (recall that the only real eigenvalue possible of K a is 0), there is x = 0 such that K a x = 0, and let us write x = x + + x -the only decomposition in I R n , so we have K a x + = -K a x -.

• If x + = 0, we consider x 2 = x -, it follows that K a x -= 0, and then t x 1 K a x 2 = 0, so φ K (A) = -1 for all x 1 ∈ F + .

• If x -= 0, we consider x 1 = x + , it follows that K a x + = 0, and then t x 1 K a x 2 = 0, so φ K (A) = -1 for all x 2 ∈ F -. • If neither x + nor x -are not zero, then we take x 1 = x + , x 2 = x -. Then (x 1 , x 2 ) are independent (they are not zero and in two supplementary vector subspaces). The calculation then gives

t x 1 K a x 2 = t (x + )K a x -= -t (x + )K a x + = 0 because K a x + = -K a x -and K a is skew-symmetric. It follows that φ K (A) = -1.
2. If K a is nonsingular, we will assume for example that k > nk, i.e., dim(F + ) > dim(F -), so dim(K a (F -)) = dim(F -) = nk, because K a is a bijective. Let us consider H the orthogonal of K a (F -) for the canonical scalar product (x|y) = t xy, then we have

dim(H) = n -(n -k) = k. Therefore dim(F + + H) = dim(F + ) + dim(H) - dim(F + ∩ H) ≤ n and thus dim(F + ∩ H) ≥ dim(F + ) + dim(H) -n = k + k -n = 2k -n > 0. Consider then 0 = x 1 ∈ F + ∩ H and 0 = x 2 ∈ F -, so K a x 2 ∈ K a (F -) and x 1 ∈ H, and thus (x 1 |K a x 2 ) = t x 1 K a x 2 = 0.
The demonstration in the opposite case k < nk is strictly similar. 3. There remains the very particular case where K a is singular, n = 2k is even, and the signature of K s is exactly (k, k).

Moreover, because 2 ≤ p < n, we have at least k ≥ 2, so n is even and ≥ 4, the cases of n = 2 and n = 3 have already been completely resolved. Let 0 = x 2 ∈ F -, z = K a x 2 , and z ⊥ its orthogonal for the canonical scalar product in I R n . This is a vector space of dimension n-

1 = 2k -1. It is known that dim F + = k, so dim(F + ∩z ⊥ ) ≥ 2k -1+k -2k = k -1, but k ≥ 2, thus dim(F + ∩ z ⊥ ≥ 1)
. So there is a vector 0 = x 1 ∈ F + ∩ z ⊥ ; then it is clear that x 1 ∈ F + and t x 1 z = t x 1 K a x 2 = 0. In all cases, the result obtained shows that, if the signature of K s is (k, nk) with 0 < k < n, then there is A ∈ M n2 with rank (A) = 2 such that φ K (A) = det( t AKA) = -1 < 0. We note that for any n and any signature of K s (of rank n), taking p = 2 gives a matrix A such that φ K (A) = det( t AKA) = -1 < 0. One might ask whether this is true for all 1 ≤ p < n.

Open problems

Two interesting open problems can be identified on the new concept of p-definite matrices.

• For 1 ≤ p < p ≤ n, does p-definite positive imply p -definite positive or in other words, D p,n ⊂ D p ,n ? The main result of the paper is the proof for p = 1. • For 1 ≤ p < p ≤ n, describe the sets D p ,n \ D p,n , D p,n \ D 1,n .

Another possible generalization

Although the proposed generalization of definite positive matrices ensues logically from mechanical considerations, another mathematical point of view should have been able to be chosen by considering that a matrix K should be p-definite positive in M n (I R) if for all A ∈ M np with rank(A) = p, t AKA would be a (1)-definite positive matrix in M p (I R). We will then call K a p * -definite positive matrix. Proposition 4: If K is (1)-de nite positive then K is p * -de nite positive. P r o o f. Suppose that K is a (1)-definite positive matrix. Let A ∈ M np with rank(A) = p be a matrix and λ ∈ I R p . Thus t λ t AKAλ = t (Aλ)K(Aλ) = t uKu > 0 with u ∈ I R n because K is a (1)-definite positive matrix.

Proposition 5: Every p * -de nite positive matrix is p-positive de nite but the converse is false. However, if K is symmetric, then p-de nite positive is equivalent to p * -de nite positive.

P r o o f. Suppose that K is a p * -definite positive matrix. Let A ∈ M np with rank(A) = p be a matrix. Then, t AKA is definite positive. By triangularizing t AKA in C, we obtain:

det( t AKA) = p k=1 | λ k | 2 r k=1 μ k
with λ k , λk , k = 1 . . . , p the complex eigenvalues of t AKA and μ k , k = 1 . . . , r its real eigenvalues. Let u k be a real eigenvector attached to μ k . Then

t AKAu k = μ k u k and 0 < t u t k AKAu k = μ k || u k || 2 .
We deduce that μ k > 0 for all k = 1 . . . , r and finally det( t AKA) > 0: K is p-definite positive.

The converse is obviously false: for example every matrix B with det(B) > 0 is not necessarly definite positive and the equivalence for symmetric matrices is easy to prove: for symmetric matrices, definite positive ⇔ p-definite positive for all 1 ≤ p ≤ n ⇔ p * -definite positive for all 1 ≤ p ≤ n.

5 Consequences on the system's mechanical stability

General considerations

From a mechanical viewpoint, we recalled in the first section that additional constraints in nonconservative elastic systems do not necessarily improve the structure's stability.

Moreover, following several recent papers [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF][START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF][START_REF] Lerbet | About the stability of nonconservative undamped elastic systems: some new elements[END_REF][START_REF] Nicot | Bifurcation and generalized mixed loading conditions in geomaterials[END_REF], a new criterion of "static" stability, generalizing the criterion of divergence for nonconservative systems may be used, which may be viewed as a transposition in structural mechanics of the second-order work criterion, well-known in geomechanics. This criterion means, mathematically, that K is 1-definite positive. The above cited studies have proved the efficiency of this criterion in stability of structures for many aspects of the stability problem. This new criterion obviously involves only the symmetric part K s of K and in [START_REF] Lerbet | About the stability of nonconservative undamped elastic systems: some new elements[END_REF], a partial relation between this criterion and the Lyapunov stability criterion was also established. As already noted, in [START_REF] Challamel | On the stability of nonconservative elastic systems under mixed perturbations[END_REF], this criterion appears as an optimum for stability by divergence of one-constrained nonconservative elastic systems.

To describe the consequences on stability, we then have to distinguish both criteria of "static" stability, which leads to three stability criteria (in the linear stability domain): a) in classical static instability or divergence (in the case of relation ( 1)), we will speak of s-instability and β s is the corresponding value. b) in linear or weak dynamic instability or flutter, we will speak of f-instability and β f is the corresponding value. c) in new criterion of static instability or mixed instability, corresponding to K(β) not (1-) definite positive, we will speak of m-instability and β m is the corresponding value. β m ≤ β s always holds (see [START_REF] Lerbet | About the stability of nonconservative undamped elastic systems: some new elements[END_REF] for a proof), but any general relation is known concerning β f .

In the present paper, the introduction of the new concept of p-definite matrices and the associated central result has two distinct consequences concerning structural stability.

• Adopting the new criterion of "static" stability, it proves that the introduction of constraints in the system does not change the "static" stability in the sense of divergence. Indeed, suppose that K(β) is 1-definite positive or equivalently that β < β m . With p additional constraints, Theorem 1 implies K(β) is p-definite positive and thus the system is still s-stable. It is then interesting to investigate if the new criterion is itself automatically satisfied for the constrained system, that is to say is Σ m-stable. Let us clarify the problem. Let Σ be a mechanical system, q e = 0 an equilibrium configuration such that the stiffness matrix of Σ at 0 is K. Suppose that K is (1-) definite positive (new criterion of "static" stability of this equilibrium) and consider now p additional constraints. Let Σ p be this new system, the equilibrium configuration assumed to be the same: we assume that the additional constraints do not disturb the equilibrium. The question is then the following: is Σ p still m-stable? In other words, considering K p the stiffness matrix of Σ p , is K p (1-) still definite positive? Eliminating p variables by the relations of constraints, Σ p is a np degrees of freedom system and K p should be a np × np matrix. However, this approach requires eliminating p variables of the system and recalculating K p . It may be easily proved that K p is (1-) definite positive if and only if K is (1-) definite positive on the np subspace of I R n defined by the relations t AX = 0 of (4), which is obviously satisfied because K is (1-) definite positive on the whole space I R n . • Without adopting the new criterion, suppose only that Σ is s-stable with only a single additional constraint. Then

Theorem 1 shows that it is automatically s-stable with p additional constraints: adding several other kinematical constraints does not change the domain of s-instability by divergence.

Moreover, another important argument of continuity highlights the usefulness of the new concept of p-definite positive matrices. In practice, the loading parameter β increases from 0 to the critical value β * . At β = 0, the system may be assumed to be conservative and stable and thus det(K(0)) > 0 and then K(0) is symmetric definite positive and thus the corresponding system is stable when submitted any to p additional kinematic constraints (Rayleigh theory). It follows that K(0) is also p-definite positive: as noted in the last proposition, for symmetric matrices all these concepts are equivalent. S-stability of a constrained system during the loading path (β increasing) is ensured by the condition of p-definite positivity: as long as K(β) is p-definite positive, Σ remains s-stable.

The last open problem (describe D p,n \ D 1,n ) has the following mechanical counterpart. Suppose that the loading parameter β takes a value β > β m such that the symmetric part K s (λ) has one negative eigenvalue. In terms of structural stability, that means that there exists a convenient kinematical constraint (p = 1, A = A 1 ∈ M n,1 (I R)) such that t ÃK(β) Ã ≤ 0 and the constrained system is no longer stable according to divergence stability. In terms of geomechanics, loading parameters exist (boundary conditions) to go into the isotropic cone. Do other additional kinematic constraints exist (p ≥ 2, A = (A 1 , . . . , A p ) ∈ M n,p (I R)) stabilizing the system? If K(β) / ∈ D p,n such additional kinematic constraints do not exist, but if K(β) ∈ D p,n \ D 1,n such additional constraints do exist.

Optimization problem

In the previous section, we showed that if a mechanical system Σ is such that its stiffness matrix K (or its symmetric part) is positive definite, then the instability by divergence may not occur for any system Σ deduced from Σ by applying kinematical constraints on Σ that one could call a subsystem Σ of Σ. The converse is an interesting problem, the answer being positive for one constraint as was shown in [START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF]. We claim here that this result is universal for any system of (linear) kinematical constraints: an n dof mechanical system Σ (K = K(β) with β an increasing loading parameter) may be kinematically constrained by p constraints (1 ≤ p < n) such that the constrained system becomes unstable by divergence if and only if β satisfies det(K s (β)) = 0 and we may also clarify the system of constraints that will to destabilize Σ. We recall that β m is defined as the smallest positive root of det(K s (β)) = 0 so that β ∈ [0, β m ] implies that K = K(β) is (1)-definite positive.

As in the previous section, we begin with the case p = 2.

Proposition 6: Let K = K(β) be a matrix 2-de nite positive. Then for all β < β m , any subsystem Σ obtained from Σ by applying two kinematical constraints on Σ is stable according to the divergence criterion, and for β = β m , there exists a system of two independent kinematical constraints such that the corresponding subsystem is divergence unstable. Moreover, the system of constraints allowing the destabilization may be explicitly given. P r o o f. From Proposition 1, we know that for all β < β m , any subsystem Σ obtained from Σ by applying two kinematical constraints on Σ is stable and that if the A ∈ G n2 subset of M n2 of matrices of rank 2, then φ K (A) > 0. Suppose now that β = β m and again use the notations of the demonstration of Proposition 1. Because the spectrum of K s (λ m ) is now included in I R + and no longer in I R * + , the Cauchy-Schwarz inequality becomes a large inequality. Let A ∈ G n2 , A = (x 1 x 2 ) with x i ∈ M n1 for i = 1, 2. Thus:

det t AK s (β m )A = t x 1 K s (β m )x 1 t x 2 K s (β m )x 2 -t x 1 K s (β m )x 2 t x 2 K s (β m )x 1 ≥ 0
and (9) becomes:

det t AK(β m )A = t x 1 K s (β m )x 1 t x 2 K s (β m )x 2 -t x 1 K s (β m )x 2 t x 2 K s (β m )x 1 + t x 1 K a (β m )x 2 2 ≥ 0.
Let us now choose x 1 = 0 ∈ Ker(K s (β m )) and

x 2 ∈ K a (β m )x ⊥ 1 . Note that K a (β m )x 1 = 0 because if K a (β m )x 1 = 0, then K(β m )x 1 = K s (β m )x 1 +K a (β m )
x 1 = 0 and K(β m ) is then singular, which is a contradiction with the main assumption that the whole system without constraints is divergence stable: K(β m ) must be regular. Thus dimV ec(K a (β m )x 1 ) = 1 and dimV ec(K a (β m )x 1 ) ⊥ = n -1. There are n -1 independent possible choices for x 2 . With such choices:

det t AK(β m )A = 0,
which signs the instability by divergence of the constrained system.

Let us now examine the general case with a family of p < n constraints. The following result holds: Theorem 3: Let K = K(β) be a p-de nite positive matrix. Then for all β < β m , any subsystem Σ obtained from Σ by applying p kinematical constraints on Σ is stable according to the divergence criterion and for β = β m , a system of p independent kinematical constraints exists such that the corresponding subsystem is divergence unstable. Moreover, the system of constraints allowing destabilization may be explicitly given. P r o o f. From Theorem 1, we know that for all β < β m , any subsystem Σ obtained from Σ by applying p kinematical constraints on Σ is stable and that if A ∈ G np subset of M np of matrices of rank p, then φ K (A) > 0.

Suppose now β = β m and suppose, for example, that dim ker K s (β m ) = 1. Let us choose x 1 = 0 ∈ ker K s (β m ). Then by a similar reasoning, as for the case p = 2, K a x 1 = 0. So dimV ec(K a (λ m )x 1 ) = 1 and dimV ec(K a (λ m )x 1 ) ⊥ = n -1. Let us choose p -1 independent vectors x 2 , . . . , x p ∈ V ec(K a (β m )x 1 ) ⊥ . Because p < n, this always possible, although x 1 itself belongs to V ec(K a (β m )x 1 ) ⊥ . Thus, the first column of t AK(β m )A is:

• on the first row:

t x 1 K(β m )x 1 = t x 1 K s (β m )x 1 = 0 because K s (β m )x 1 = 0; • on the other rows (j = 2, . . . , p): t x 1 K(β m )x j = t x 1 K s (β m )x j + t x 1 K a (β m )x j = t x j K s (β m )x 1 + t x 1 K a (β m )x j = 0 + 0 = 0 because K s (β m )x 1 = 0 and x j ∈ V ec(K a (β m )x 1 ) ⊥ .
The first column of t AK(β m )A is then nil and thus:

det t AK(β m )A = 0,
which signs the instability by divergence of the constrained system. Note that fortunately the result cannot be extended to p = n: with n constraints, the system cannot move and is stable. Remark also that this optimization problem could appear as a geometric problem of the optimization of an orthogonal invariant function φ on a Grassmann manifold like in [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF]. Indeed, the function φ

K is (right)O(p)-invariant because if Q ∈ O(p) we have: φ K (AQ) = det( t Q t AKAQ) = det( t Q) det( t AKA) det(Q) = det( t Q) det(Q)φ K (A) = det( t QQ)φ K (A) = det(I p )φ K (A) = φ K (A)
but here we do not need to look for a numerical algorithm because we have got an exact algebraic solution due to the special form of the objective map φ K .

To summarize, the second-order work criterion appears as the optimum criterion so that a divergence stable system Σ remains divergence stable when it is submitted to any system of kinematical constraints.

A 3 d.o.f. system

To illustrate the previous results, we come back to the first motivating example, we apply the previous results and we take advantage of it to answer the last questions of Sect. 2 on hold.

Instability of the constrained system and second-order criterion

The symmetric part K s (β) of K(β) is:

K s (β) = ⎛ ⎜ ⎝ 2 -β -1 β 2 -1 2-β -1 + β 2 β 2 -1 + β 2 1 ⎞ ⎟ ⎠ .
The roots of det(K s (β)) = 0 are 1, 1 -√ 3, 1 + √ 3 and thus, with previous notations, β m = 1, which leads to

K s (1) = ⎛ ⎜ ⎝ 1 -1 1 2 -1 1 -1 2 1 2 -1 2 1 ⎞ ⎟ ⎠
denoting as previously by V ec(u 1 , . . . , u r ) the vector space generated by the family of vectors u 1 , . . . , u r , calculations give:

Ker(K s (1)) = V ec(x 1 = ⎛ ⎜ ⎝ 1 1 0 ⎞ ⎟ ⎠). The skew symmetric part of K(1) is K a (1) = ⎛ ⎜ ⎝ 0 0 1 2 0 0 1 2 -1 2 -1 2 0 ⎞ ⎟ ⎠ and thus the subspace V ec(K a (1)(x 1
)) generated by the vector K a (1)(x 1 ) is the one-dimensional subspace V ec(

⎛ ⎜ ⎝ 0 0 -1 ⎞ ⎟ ⎠) and its orthogonal is the plane V ec( ⎛ ⎜ ⎝ 1 0 0 ⎞ ⎟ ⎠ , ⎛ ⎜ ⎝ 0 1 0 ⎞ ⎟ ⎠).
Following the algorithm given by the proof of Theorem 3, let us choose

x 2 in this plane such that (x 1 , x 2 ) is a free family. For example,

x 2 = ⎛ ⎜ ⎝ 1 0 0 ⎞ ⎟ ⎠. Now let A = (x 1 x 2 ) = ⎛ ⎜ ⎝ 1 1 1 0 0 0 ⎞ ⎟ ⎠.
Calculations give t AK(1)A = 0 0 0 0 and obviously det( t AK(1)A) = 0. Be careful that the matrix A is actually the matrix à of Sect. 3.2. To come back to the matrix of kinematical constraints we have to calculate K(1)A which leads to the matrix ⎛ ⎜ ⎝ 0 1 0 -1 -1 0 ⎞ ⎟ ⎠ meaning that the 3 dof Ziegler system constrained by θ 1θ 2 = 0 and -θ 3 = 0 is divergence unstable for β = 1 which may be immediately checked. MAPLE may solve the complete problem, i.e., the nonlinear system

Optimization problem

∂φ K(β) (α i , δ i ) ∂α i = 0, ∂φ K(β) (α i , δ i ) ∂δ i = 0, i = 1, 2, 3.
The nine families of solutions are:

S 1 = {β = β, α 1 = α 1 , α 2 =
S 8 = β = -2 δ 1 δ 2 + 2 δ 1 2 + δ 2 2 δ 1 (δ 1 -δ 2 ) , α 1 = α 1 , α 2 = α 1 δ 2 δ 1 , α 3 = α 1 δ 3 δ 1 , δ 1 = δ 1 , δ 2 = δ 2 , δ 3 = δ 3 , S 9 = β = 1 + √ 3, α 1 = α 1 , α 2 = α 2 , α 3 = α 2 + α 1 , δ 1 = 1/2+1/2 3 + 2 √ 3 δ 2 , δ 2 = δ 2 , δ 3 = δ 2 + 1/2+1/2 3 + 2 √ 3 δ 2 .
For S 1 , S 2 , S 3 , S 5 , S 8 , rank(A) = 1 the vectors x 1 and x 2 are colinear and these solutions must be rejected. for S

Conclusion

In this paper, we study the properties of p-definite positive matrices. This new concept of linear algebra, generalizing the concept of definite positive matrices, is defined for investigating bifurcations or "static" instabilities of nonconservative mechanical systems submitted to additional kinematical constraints. The main result of the paper is the proof of the fact that any positive definite matrix is p-definite positive for all 1 ≤ p ≤ n. In the last part, the consequences for stability of constrained mechanical systems are highlighted and illustrated on a 3 dof system. The second-order work criterion then appears as the optimum criterion for a divergence stable system for remaining divergence stable if subjected to any family of kinematical constraints.

Fig. 1 3

 1 Fig. 1 3 dof Ziegler system.

  Let us put A = (x 1 x 2 ) = → φ K (A) = φ K(β) (A) may so be considered as a function of six variables α 1 , δ 1 , α 2 , δ 2 , α 3 , δ 3 plus the loading parameter β. The optimization problem is to calculate:inf (αi,δi)∈U φ K(β) (α i , δ i )and to prove that β = 1 is the lower positive value of the loading parameter. U is the open set of I R 6 such that rank(A) = 2.

α 2 , α 3 = 3 = α 1 δ 3 δ 1 , δ 1 =α 2 δ 1 δ 2 , α 2 = α 2 , α 3 = δ 3 α 2 δ 2 , δ 1 =

 23311222321 α 3 , δ 1 = 0, δ 2 = 0, δ 3 = 0} , S 2 = {β = β, α 1 = 0, α 2 = 0, α 3 = α 3 , δ 1 = 0, δ 2 = 0, δ 3 = δ 3 } S 3 = β = β, α 1 = α 1 , α 2 = 0, α δ 1 , δ 2 = 0, δ 3 = δ 3 , S 4 = β = 1+ √ 3, α 1 = α 1 , α 2 = α 2 , α 3 = α 2 + α 1 , δ 1 = δ 1 , δ 2 = 0, δ 3 = δ 1 , S 5 = β = β, α 1 = δ 1 , δ 2 = δ 2 , δ 3 = δ 3 , S 6 = β = 1+ √ 3, α 1 = α 1 , α 2 = α 2 , α 3 = α 2 + α 1 , δ 1 = -δ 2 + δ 3 , δ 2 = δ 2 , δ 3 = δ 3 , S 7 = {β = 1, α 1 = α 1 , α 2 = α 2 , α 3 = 0, δ 1 = δ 1 , δ 2 = δ 2 , δ 3 = 0} ,

4 and S 9 , β = 1+ √ 3 1 .

 931 corresponds to the other two values of the load parameter, which are roots of det(K s )(β) = 0: they are < 0 or > Thus the solution S 7 with β = 1 is the optimal solution and the value of A (A = (x 1 x 2 ) = above corresponds to a possible value of solutions of S 7 : α 1 = 1, α 2 = 1, α 3 = 0, δ 1 = 1, δ 2 = 0, δ 3 = 0. Figure2plots on the z-axis the values of the partial function (δ 2 , β) → φ K(β) (1, 1, 0, 1, δ 2 , 0), which shows the optimal value β = 1 for all δ 2 as a particular case of the solution S 7 .

Fig. 2 (

 2 Fig. 2 (online colour at: www.zamm-journal.org) α1 = 1, α2 = 1, α3 = 0, δ1 = 1, δ3 = 0, δ2, β variable.