domain: spi info

text:
Mayank Shekhar Jha
email: mayank-shekhar.jha@univ-lorraine.fr

Hugues Garnier
email: hugues.garnier@univ-lorraine.fr

Didier Theilliol
email: didier.theilliol@univ-lorraine.fr

Redundancy-Aware Physics Informed Neural Networks (R-PINNs) based Learning of Nonlinear Algebraic Systems with Non-Measurable States

The paper presents Redundancy-Aware Physics Informed Neural Networks (R-PINNs) for learning of unknown model parameters of nonlinear algebraic systems in continuous time with non-measurable state variables. R-PINNs accomplish the learning task in presence of non-measurable states of the system by incorporating input-output representation of the a priori available physics based laws, generally in form of nonlinear differential (partial) equations within the NN based learning procedure, leading to learning of a set of optimal parameters that determine the optimal mapping between input-output data while adhering to the known physics. Analytical Redundancy Relationships (ARRs) are able to express input-output representation of system using solely the measured/known variables by exploring the redundancy within the analytical structure of the system. The paper proposes a methodology that includes ARR derivation and suitable integration within PINNs framework to develop R-PINNs. Mathematically rigorous novel proofs on uniform and ultimate boundedness (UUB) of the output and parametric estimation errors in Lyapunov sense is provided. Finally a DC motor enabled friction drive system based simulation study is presented to demonstrate the effectiveness of the approach.

I. INTRODUCTION

The predictive ability of physics based models is dependent on the available knowledge of governing laws, identification / calibration of model parameters under given conditions and the model uncertainties, to name a few. With the advent of modern deep learning, there has been an unprecedented progress in the field of data driven modelling [START_REF] Kutz | Dynamic mode decomposition: data-driven modeling of complex systems[END_REF] and nonlinear identification [START_REF] Lusch | Deep learning for universal linear embeddings of nonlinear dynamics[END_REF]. Two factors that remain extremely crucial for accurate learning are: large high quality data sets (rich and diverse) and immense computing power. While the former leads to variance reduction, the latter, allows for large models to be trained leading to a certain bias reduction. However, the availability of usable high quality data is not always guaranteed leading to less generalization capability and high variance. Moreover, model calibration under limited data introduces high bias. To reduce model bias and increase generalization capability, various approaches have been proposed including learning of governing equations and generalised coordinates [START_REF] Champion | Data-driven discovery of coordinates and governing equations[END_REF], semi physical models and recently, the focus has been on the incorporation of physics based knowledge within the neural network (NN) structure [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF].

Past semi decade has witnessed an ever growing interest and extensive work using physics informed neural networks (PINNs) for the discovery of governing laws (forward problems) as well as learning of model parameters using the available data (inverse problems) in several fields of engineering including discovery of Navier-Stokes equation [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], Schrodinger equation [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], Burger [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], Hamilton and Euler-Lagrange equation [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF] to name a few. Readers are referred to [START_REF] Cuomo | Scientific machine learning through physics-informed neural networks: where we are and what's next[END_REF] for the latest comprehensive review on PINNs and their application in various engineering domains starting from 2017. Most existing NN(deep) based approaches for parametric learning or system identification, target learning of the mapping function between input and output data (black-box approach) that neglect inclusion of physics laws within learning paradigm [START_REF] Ljung | Deep learning and system identification[END_REF] [START_REF] Lusch | Deep learning for universal linear embeddings of nonlinear dynamics[END_REF]. On the other hand, only very few notable works based on variants of PINNs have appeared recently for system identification [START_REF] Stiasny | Physics-informed neural networks for non-linear system identification for power system dynamics[END_REF] . In general, PINNs based approaches accomplish the learning task by incorporating the a priori available physics based laws, generally in form of nonlinear differential (partial) equations within the NN based learning procedure, leading to learning of a set of optimal parameters that determine the optimal mapping between input-output data, while adhering to the known physics. Generally the feature that makes a classical NN into PINNs is the incorporation of nonlinear differential(partial) equation(s) within the loss function leading to minimisation of the corresponding residual. The loss is indicative of consistency (measure of difference from zero). The latter must be minimised to accomplish the learning task. These residuals are constructed using system variables that are usually assumed to be measurable/available. However, in reality, dynamical systems may have only a set of states measurable. In this context, existing PINNs approaches neglect the problem arising from the presence of unknown/non-measurable state variables. This becomes crucial especially for the construction of nonlinear differential residual as the presence of unknown state variables calls for their estimation along with unknown parameters. As per the authors, there is no existing work that formally enables inclusion of only known (measurable) information along with unknown parameters within PINNs for parametric learning.

On the other hand, Analytical Redundancy Relationships (ARRs) are able to express input-output representation of system using solely the measured/known variables, by exploring the redundancy within the analytical structure of the system [START_REF] Staroswiecki | Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems[END_REF]. Under nominal functioning of the system, ARRs must be consistent and as such, have been extensively exploited for residual generation. ARRs have been extensively investigated in Fault detection and isolation (FDI) literature wherein the main utility remains in generation of residuals using only the measured variables [START_REF] Zhang | Fault detection and isolation in nonlinear dynamic systems: A combined input-output and local approach[END_REF] [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF]. It must be noted that ARRs are closely related to input-output representation of dynamic systems and as such, ARR derivation techniques have been heavily inspired by the former that involve only measured variables of system [START_REF] Zhang | Fault detection and isolation in nonlinear dynamic systems: A combined input-output and local approach[END_REF]. For nonlinear systems, ARR derivation principally involves state elimination techniques [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF], broadly based on three methods: elimination theory that rests primarily on Euclidean division and successive derivation [START_REF] Diop | Elimination in control theory[END_REF], Grobner bases which uses Euclidean division and the computation of the so-called S-polynomials [START_REF] Da Cox | Ideals, Varieties, and Algorithms[END_REF] and Characteristic sets (also called Ritt's algorithm) that involves direct elimination of state to find minimum order differential input-output constraints [START_REF] St Glad | Nonlinear state space and input output descriptions using differential polynomials[END_REF]. In the case of algebraic dynamic systems, elimination theory provides the most straightforward and intuitive tool for state elimination [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF].

To bridge the existing scientific gap, this paper blends the benefits of ARR with PINNs framework in order to address the parametric learning problem for dynamical systems with non-measurable state variables (section II). To that end, as novelty, the paper presents Redundancy aware PINNs (R-PINNs) for nonlinear algebraic systems in continuous time with non-measurable state variables (section III). Therein, the redundancy within system variables is exploited to construct ARRs which are then suitably integrated within PINNs structure to learn unknown parameters. The paper also presents mathematically rigorous novel proofs on uniform and ultimate boundedness (UUB) of the output and parametric estimation errors in Lyapunov sense and provides a DC motor enabled friction drive system based simulation study to demonstrate the effectiveness of the approach (section IV).

II. PROBLEM FORMULATION

Consider a dynamical system in continuous time with physics laws given by a set of nonlinear differential algebraic equations (DAEs) [2][11] leading to state space model as:

ẋ = f (x, u, λ) y = g(x, u) (1)
where x ∈ R n ∈ Ω is the state variable vector, which is non available, Ω being a compact set, ẋ is the time derivative of state variable vector, u ∈ R m ∈ Ω is the control input vector with λ ∈ R d are the unknown parameter vector, y ∈ R p is the output measurement vector. This work considers f (•) and g(•) as nonlinear DAEs such that f (•) and g(•) are C 1 where C 1 denotes the set of continuously differentiable functions over the compact set Ω.

Remark 1. Restricting the known system dynamics of (1) to nonlinear DAEs makes it possible to use mathematical properties from differential algebra. Moreover, the assumption is not as restrictive as non-polynomial non-linearity can be transformed to equivalent polynomial DAEs. For example, y =cos(x) can be equivalently written as a DAE as y 2 ẋ2 + ẏ2 = ẋ2 .

1) Known and unknown variables:

The system variables sets Z = x ∪ y ∪ u can be decomposed in two categories: measurable variables K = y ∪ u and unknown nonmeasurable variable(s) X = x.

2) Residual Generation: Using traditional PINNs based approach, the nonlinear differential residual corresponding to system (1) will be generated as:

ẋ -f (x, u, λ) (2)
Readers are referred to [START_REF] Cuomo | Scientific machine learning through physics-informed neural networks: where we are and what's next[END_REF] for more details. The existing PINNs based approaches oblige the implication of state variables x into the residual expression. As such, in presence of unknown(partially) variables x which is the case for most of the nonlinear dynamic systems, it becomes imperative to resort to state estimation methods in order to facilitate parametric estimation approaches.

III. PROPOSED METHODOLOGY

To alleviate the difficulty arising from unknown state variables, this paper proposes a variable elimination based approach via generation of ARRs by exploiting the analytical redundancy in the system mathematical model [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF]. Elimination of x leads to transformation of the original mathematical model to an equivalent input-output representation of the system. To that end, this paper proposes ARR derivation followed by construction of R-PINNs.

A. ARR derivation with known inputs

Consider the system (1), wherein the inputs are known and all functions are differentiable (as (1) is DAE). Then, the steps below are followed to obtain the ARRs [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF].

Step I: Derivation of outputs. It is possible to obtain the total derivative of the output signal as:

ẏ = ∂g(•) ∂x ẋ + ∂g(•) ∂u u (3)
Substituting ( 1) in (3) leads to:

ẏ = ∂g(•) ∂x f (x, u, λ) + ∂g(•) ∂u u := g 1 (x, ū(1) , λ) (4)
where for the sake of ease of notation, ū(1) = (u ′ , u′ ) ′ . Then, through iterative derivation up-to some order of derivation q one obtains:

ȳ(q) = G (q) (x, ū(q) , λ) (5)
with ū(q) ∈ R (q+1)m , is set of (q + 1)p DAE based constrains that must be satisfied under nominal functioning of the system.

Step II: Elimination of the state. It is assumed that (q + 1)p > n and the Jacobian ∂G (q) (•) ∂x is of rank n (see [START_REF] Staroswiecki | Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems[END_REF] [START_REF] Diop | Elimination in control theory[END_REF] for a more mathematically rigorous treatment omitted here). Then, (5) can be decomposed into two subsystems:

ȳ(q) I ȳ(q) II - G (q) I (x, ū(q) , λ) G (q) II (x, ū(q) , λ) = 0 (6)
where subsystem (ȳ

(q) I , G (q)
I (x, ū(q) , λ)) ∈ R n and (ȳ

(q) II , G (q)
II (x, ū(q) , λ)) ∈ R (q+1)p-n Then, x can be computed as a function of other variables that are known leading to:

x = ϕ(ȳ (q) I , ū(q) , λ) (7)
where ϕ(•) ∈ R n is some nonlinear DAE. Further, using implicit function theorem [START_REF] George | The implicit function theorem: history, theory, and applications[END_REF] one can substitute ( 7) in ( 6) to obtain (5) sensitive to only known variables as:

r(ȳ (q) , ū(q) , λ) = ȳ(q) II -G (q) II ϕ(ȳ (q) I , ū(q) , λ), ū(q) , λ = 0 ( 8
)
where r(ȳ (q) , ū(q) , λ) is a set of constraints or ARRs that are sensitive to only known system variables and unknown parameters that must be satisfied.

B. R-PINNs

To address the inverse problem for system (1), i.e. estimation of unknown parameters λ given measured data y, this paper proposes the latter, along with ARRs (input-output representation of original system (1)) to constitute R-PINNs as:

y = h(z)
r(ȳ (q) , ū(q) , λ) = 0

where z := [u 1 , u 2 , ....u m ; t] is the input space-time vector with z ∈ Ω , h(•) ∈ R m is a hidden nonlinear function that maps the available measured data to input time series and r(•) are the derived set of ARRs in [START_REF] George | The implicit function theorem: history, theory, and applications[END_REF].

To solve [START_REF] Kutz | Dynamic mode decomposition: data-driven modeling of complex systems[END_REF], this paper proposes R-PINNs that considers as input, z to computationally predict the hidden relation h(z) using a NN, parameterised by a set of parameters θ ∈ R (m+1)×p leading to the approximation:

ŷθ = ĥθ (z) (10)
where (•) θ denotes the NN approximation realised with parameter set θ. Further, since ARRs must be satisfied at all times, an approximation of r(•) can be generated using ŷθ and its (partial)derivative(s) as:

rθ,λ ( ŷ(q) θ , ū(q) , λ) = 0 (11)
Consider the R-PINNs parameter set Θ as tuple of shared parameters Θ = (θ, λ) that consists of NN approximation parameters θ and unknown system parameters λ. Then, an R-PINNs parameter set Θ * can be found through minimisation of total loss funciton L(Θ) as:

Θ * = arg min Θ L(Θ) (12)
where the R-PINNs based loss function can be defined as:

L(Θ) = (w y L y + w r L r ) (13)
considers the weighted sum of: loss L y due to true measurement data y and its respective estimation by NN; and loss L r based on the corresponding estimated residual which must be equal to zero [START_REF] Ljung | Modeling and identification of dynamic systems[END_REF], with w y and w r being the weight coefficients.

The building blocks of proposed R-PINNs are shown in Fig. 1. The loss L y represents the mismatch between the available measurement data y and its approximation by NN ŷθ and attempts to fit the data to learn accurately the mapping ĥθ over Ω . On the other hand, L r represents the difference of estimated residual (value of ARRs) from zero in order to "orient" the learning in such a way that known physical laws are adhered to. L y and L r can be expressed as:

L y (Θ) = Ω y -ĥΘ (z)dz L r (Θ) = Ω rθ,λ ( ŷ(q) θ , ū(q) , λ)dz (14)
Consider N d measurement data points such that {z i , y i } N d i=1 are the observation points. Then, typically, these losses can be numerically obtained using mean squared error (MSE) formulation as:

L y (Θ) = M SE y = 1 N d N d i=1 y i -ĥθ (z i ) 2 L r (Θ) = M SE r = 1 N d N d i=1 ( ŷ(q) θ , ū(q) , λ) 2 (15)
As the ARRs are sensitive to only known/measured data, (12) can be estimated over the domain Ω. It should be noted that the derivatives (partial) required to construct the ARRs (11) can be obtained directly using Automatic Differentiation (AD)procedure [START_REF] Cuomo | Scientific machine learning through physics-informed neural networks: where we are and what's next[END_REF][14] that evaluates derivatives(partial) in an algorithmic fashion through analytical evaluation. Usage of AD is extensively common in the state of the art works associated with NN, to accomplishing back-propagation based NN weight update. Almost all available Deep Learning packages, such as Pytorch/TensorFlow, provide numeric estimates of derivatives with high accuracy.

C. Boundedness of errors

Universal approximation property of feed-forward NN, enables expression of output of R-PINNs y Θ (z) in an exact manner on a compact set Ω, ∀z ∈ Ω:

y = Θ * T σ(z) + ε No (z) (16)
where

Θ * = (θ * , λ * ) ∈ R No is the ideal weight vector, N o is the number of neurons, σ(z) = (σ 1 (z)), σ 1 (z), ...σ No (z)) T ∈ R No
is the nonlinear activation function (sigmoid or tanh) with σ j (z) ∈ C 1 (z), σ(0) = 0, and the set of activation functions {σ j (z)} j=No j=1 are linearly independent ( to assure persistence of excitation) and ε N o is a bounded R-PINNs function that represents the reconstruction error. The derivative of 16 can be given as:

ẏ = ∇σ T (z)Θ * + ∇ε No (z) (17)
where ∇σ(z) = ∂σ(z) ∂(z) and ∇ε No (z) = ∂ε N o(z) ∂(z) . R-PINNs based output estimate can be expressed as:

ŷ = ΘT σ(z) ( 18
)
where Θ is the estimation of Θ * . Denote weight estimation and output estimation errors respectively as Θ = Θ * -Θ and ỹ = y * -ŷ. The parameters Θ of PINNs are updated using gradient descent as:

Θ = -l r ∂L( Θ) ∂ Θ ( 19
)
where 0 < l r < 1 is a non negative learning rate.

∂ Θ || ≤ δ L ⇕ .
It is noted that these assumptions are not as restrictive as various saturating nonlinear activation functions that satisfy these properties such as sigmoid or tanh functions exist. Theorem 1. Under the assumptions 1 and 2, given the R-PINNs weight update law as (19), then the R-PINNs based estimation error as well as parameter estimation errors are guaranteed to be uniformly and ultimately bounded (UUB).

Proof: Consider a Lyapunov function candidate:

L(z) = L 1 + L 2 (20)
where

L 1 = 1/2ỹ T ỹ L 2 = 1/2 ΘT Θ (21)
Taking derivative of L 1 we have:

L1 = ỹT ẏ = (y -ŷ) T .( ẏ -ẏ) = (Θ * T σ(z) + ε N o -Θσ(z)) T (∇σ T (z)Θ * + ∇ε N o -∇σ T (z) Θ -ΘT σ(z)) = ( Θσ(z) + ε N o ) T (∇σ T ( θ) + ∇ε N o -ΘT σ(z)) = σ T (z) ΘT ∇σ T (z) Θ + σ T ΘT ∇ε N o + σ T ΘT l r ∂L( Θ) ∂ Θ σ(z) + ε T N o ∇σ T Θ + ε T N o ∇ε N o + ε T N o l r ∂L( Θ) ∂ Θ σ(z) (22)
Using the Cauchy-Schwartz inequality a T b ≤ (1/2)a T a+ (1/2)b T b for any two vectors a and b and some mathematical manipulations, (22) yields :

L1 ≤ ||σ(z)||.||∇σ(z)||.|| Θ|| 2 + ||σ T (z)||.( ΘT Θ + ∇ε T N o ∇ε N o ) + ||σ T (z)||.(Θ T Θ + σ T (z)σ(z))l r || ∂L( Θ) ∂ Θ || + ε T N o (∇σ T ∇σ + ΘT Θ) + ||ε N o ||.||∇ε no || + ||ε||.||l r ∂L( Θ) ∂ Θ ||.||σ(z)|| (23)
Then, taking into account the assumptions 1 and 2, following is obtained:

L1 ≤ (σ M + σ ∇M + σ M + σ M )|| Θ|| 2 + σ M σ 2 ∇M + σ 3 M l r δ L ⇕ + ε M σ 2 ∇M + ε M ε ∇M + ε M l r δ L ⇕ σ M (24)
Similarly, consider derivative of L 2 as:

L2 = ΘT Θ = (Θ * -Θ) T .( Θ * -Θ) = (Θ * -Θ) T (l r ∂L( Θ) ∂ Θ ) (25)
Using Cauchy-Schwartz inequality and some mathematical manipulations, one obtains:

L2 ≤ (1/2)|| Θ|| 2 + (1/2)|| ∂L( Θ) ∂ Θ || (26)
Putting ( 24) and ( 26) into (20), one obtains L(z) ≤ 0 which implies that ỹ and Θ are UUB over the compact set Ω if ỹ(0) = 0 and Θ(0) = 0 are bounded over Ω. The latter being true, the proof is completed hereby.

IV. SIMULATIONS AND RESULTS

To demonstrate the effectiveness of the approach, a friction drive system composed by a driver device (a DC motor) and a driven device (a wheel) is considered for simulation based analysis. The system dynamics and details can be consulted in [START_REF] Rodriguez Obando | Deterioration estimation for predicting and controlling RUL of a friction drive system[END_REF] . Here, the DC motor system along with nonlinear friction drive dynamics is considered in form of DAE with states x(t) = [i(t), ω 1 (t), ω 2 (t)] T as:

L di dt = U -Ri -K e ω 1 (27) J 1 dω 1 dt = K e i -b 1 ω 1 + T L (28) J 2 dω 2 dt = T s -b 2 ω 2 (29)
where U (t) is the input voltage, R=0.08 Ohms is rotor resistance, i(t) is the motor current K e =0.1 being the electromotive force constant, L=0.005 Henry being the electric inductance, J 1 and J 2 are respectively moments of inertia of the driver and the driven device, ω 1 and ω 2 are respectively the angular speeds of driver and driven device and b 1 , b 2 are the viscous friction coefficiencts, T L is a load torque seen from the driver side, T S is a source of torque observed from the driven side. Both torques can be written in terms of the contact force F c as:

T L (t) = -F c r 1 T s (t) = F c r 2 (30)
where F c is the contact force proportional to the relative tangential speed at the contact level as:

F c = α(r 1 ω 1 -r 2 ω 2 ) 2 (31)
with r 1 , r 2 being the radii of driver and driven device respectively, α being a model parameter known as contact quality coefficient generally non-measurable and must be estimated. Moreover, viscous friction coefficient b 1 is considered as the unknown parameter. Thus, unknown parameter set is λ = [α, b 1 ] T . Motor current i(t) and driven device speed ω 2 are assumed measurable, leading to the following output measurement vector as:

y(t) = [i(t), ω 2 (t)] T (32)
The objective is to estimate α and b 1 using some available data y(t) and ARR based input-output representation using R-PINNs. ARRs are derived through elimination of nonmeasurable variable ω 1 . To that end, ( 27) is considered so that ω 1 = (L di dt -U + Ri)/K e and substituted in (28) as well as ( 27) is differentiated with respect to time to give ω1 and substituted in (28) leading to ARR1 as:

ARR1 : K e i + (1/K e )(U -Ri -L di dt )(-b 1 -r 2 α) -(J 1 /K e )( U -R di dt -L d 2 i dt 2 ) -ω 2 r 1 r 2 α = 0 (33)
Similarly, substituting ω 1 from ( 27) in (28) leads to ARR2 as:

ARR2 : -b 2 ω 2 -J 2 dω 2 dt + r 2 α(r 1 (1/K e )(U -Ri -L di dt ) -r 2 ω 2 ) = 0 (34)
As can be seen, ARR1 and ARR2 are sensitive to unknown parameters λ = [α, b 1 ] T and only the measurable state variables. A simulated data set comprising of N d = 10000 data points is considered for learning. The simulation is done in continuous time using Runge-Kutta solver with fixed step of integration as 0.0001.

Fig. 1 :

 1 Fig. 1: Proposed R-PINNs Architecture. Remark 2. R-PINNs based methodology estimates the R-PINNs shared parameter set Θ = (θ, λ) simultaneously to give optimal parameters corresponding to NN based approximation error and residual based approximation error.

Assumption 1 .

 1 The parameters of R-PINNs are upper bounded by known constants such that ||Θ|| ≤ Θ M and ∀z ∈ Ω, ||σ(z)|| ≤ σ M , ||∇σ(z)|| ≤ σ ∇M . Moreover, the R-PINNs reconstruction errors are upper bounded by known constants such that ||ε No || ≤ ε M and ||∇ε No || ≤ ε ∇M Assumption 2. The gradient of loss L(Θ) with respect to parameters is bounded such that || ∂L( Θ)

( 13 ) 2 and

 132 constitutes of L y (Θ) = (1/N d ) L r (Θ) = (1/N d ) N d i=1 ∥ARR1 + ARR2∥2 with the importance weights w y = w r = 0.5 to attribute equal importance to the former two.

Fig. 2 :Fig. 3 :Fig. 4 :

 234 Fig. 2: Total loss per epoch of training

Fig. 5 :

 5 Fig. 5: Estimation of measured variables using R-PINNs

  Additive white Gaussian noise with signal to noise ratio (SNR) as 50db is considered as noise over measured variables. System input U (t) with white additive Gaussian noise with SNR 50 db, in form of step signals with variable amplitudes and period is injected within the system model to generate the data set. To approximate the output (32), a fully connected NN having 8 hidden layers with neurons in-order as [188,144,64,32,16,8,4,2] and all activation functions as Hyperbolic tangent tanh(•) is considered, leading to the approximation of y(t) as ŷ (t). Then, the total loss function

ACKNOWLEDGMENT

Authors are thankful to anonymous reviewers for their constructive remarks and Soha Kanso for her help in compilation of the compliant version of this paper.


domain: spi spi

text:
Mayank Soha Kanso
email: soha.kanso@univ-lorraine.fr

Shekhar Jha
email: mayank-shekhar.jha@univ-lorraine.fr

Didier Theilliol
email: didier.theilliol@univ-lorraine.fr

Degradation Tolerant Control Learning for Discrete-Time Affine Nonlinear Systems

Keywords: Adaptive Dynamic Programming, Reinforcement Learning, Actor-Critic Structure, Optimal Control, Fault Tolerant Control, Safety-Critical Systems, Affine Nonlinear Systems

This paper develops a degradation tolerant optimal control in the framework of Reinforcement Learning (RL). Safety-critical and mission-critical systems require the development of new control designs that maintain system stability and performance specifications but also address incipient degradation. The aim of this work is to decelerate the speed of degradation by minimizing a cost function that includes the rate of evolution of degradation and the performance requirements. The controller is developed for discrete-time nonlinear systems affine in control, where the system's states are affected by a nonlinear degradation. Value iteration (VI) algorithm based approach is developed to find suitable approximations of both optimal control policy and optimal cost, while guaranteeing closed-loop stability and minimization of degradation rate. Offline model-based Adaptive Dynamic Programming (ADP) algorithm is developed and implemented using actor-critic structure which involves training of both actor and critic neural networks (NN). After training the actor NN with the optimal policy, the NN is implemented in real time to generate the input of the system. Simulation example shows the efficiency and feasibility of the algorithm.

INTRODUCTION

Traditional control system design [START_REF] Stengel | Optimal control and estimation[END_REF]), [START_REF] Åström | Adaptive control[END_REF]) focuses only on the stability and the performance without taking into consideration the effects of aging, fatigue, and damage of the concerned components and without minimizing the risk of failure. However, safety-critical systems [START_REF] Knight | Safety critical systems: challenges and directions[END_REF]) arise in several application areas, such as transportation and airtraffic control systems, space systems, nuclear plants and automated industrial processes. The evolution of such complex systems requires the implementation and the development of new control technologies that maintain system stability and performance specifications, and also address incipient fault and graceful performance degradation.

Recent approaches have applied modern control techniques such as adaptive or robust control to address situations where the degree of failure may be unknown. In [START_REF] Bole | Fault adaptive control of overactuated systems using prognostic estimation[END_REF]), a fault adaptive control is proposed for incipient fault modes growing to catastrophic failure conditions. The methodology is developed for a finite constrained optimization problem where the model of the system and the degradation is supposed to be known. Moreover, fault tolerant control design [START_REF] Noura | Fault-tolerant control systems: Design and practical applications[END_REF]), [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF]) has been developed for various industrial, mission critical and safety critical systems that operate in closed loop, in order to compensate for fault occurrence. More recently new methods are expanded such that useful life of critical systems can be extended. In this context, health aware control has recently become one of the domains where control is being designed, based upon the state of health (SoH) and/or Remaining Useful Life (RUL) prognostics of critical components. Some prominent works have proposed methods to develop control laws that attempt to extend the RUL of component/system such as [START_REF] Pour | Healthaware control design based on remaining useful life estimation for autonomous racing vehicle[END_REF], [START_REF] Salazar | System reliability aware model predictive control framework[END_REF]). Also, in the framework of model predictive control (MPC), several works were adopted to produce a controller that ensures robustness to particular failures, thus reducing their impact on the system [START_REF] Brown | A prognostics enhanced reconfigurable control architecture[END_REF], [START_REF] Brown | Prognostics enhanced reconfigurable control of electro-mechanical actuators[END_REF]). However, one of the limitations of previous approaches is that it requires a known dynamic of the system. Reinforcement learning has achieved remarkable success for complex systems with unknown dynamics [START_REF] Lewis | Reinforcement learning and approximate dynamic programming for feedback control[END_REF]).

Reinforcement learning refers to agent or actor interacting with its environment and modifying its actions or control policies, based on stimuli received in response to its actions [START_REF] Buşoniu | Reinforcement learning for control: Performance, stability, and deep approximators[END_REF]). RL theory is well founded on the principles of optimal control theory based on Dynamic Programming (DP). It allows to solve optimization problems using the principle of optimality from DP that provides an essential foundation for understanding RL [START_REF] Lewis | Reinforcement learning and adaptive dynamic programming for feedback control[END_REF]). In particular, RL can handle optimal control problems for unknown nonlinear systems. Most of the RL methods attempt to achieve the same goal as DP, with less computation and without knowing the model of the environment. The combination of DP and neural networks (NN) as approximators leads to adaptive dynamic programming (ADP) approaches, introduced by Paul Werbos in 1977 [START_REF] Werbos | Advanced forecasting methods for global crisis warning and models of intelligence[END_REF]). In the framework of health aware control design based on RL, very few works have been proposed. [START_REF] Jha | A reinforcement learning approach to health aware control strategy[END_REF]) presents a RL based approach that is employed to learn an optimal control policy in face of component degradation by using global system transition data and RUL predictions. [START_REF] Yousefi | Reinforcement learning for dynamic condition-based maintenance of a system with individually repairable components[END_REF]) develops a dynamic optimal policy for multi-component systems with individually repairable components using a Q-learning algorithm. To fill in this scientific gap, this paper proposes a learning strategy for degradation tolerant optimal control of nonlinear unstable systems.

In the previous work [START_REF] Kanso | Degradation tolerant optimal control design for linear discrete-times systems[END_REF]), a linear quadratic regulator (LQR) and tracker (LQT) was designed, over finite and infinite time horizon, for a discretetime linear system in the presence of a linear degradation. This paper aims to address nonlinear cases. The main contribution of the paper is to learn a control law that decelerates the speed of degradation while maintaining optimal performance and stability of a nonlinear system. This work examines a degradation tolerant approach for discretetime nonlinear systems, where the states of the system are considered to be affected by a nonlinear degradation. The approach is developed in the framework of model-based RL, by building a quadratic cost function that includes the speed of degradation and by implementing the Actor-Critic structure. The optimal policy and value function are learned offline using Value Iteration (VI) algorithm and the model of the system, then the trained NNs are implemented in real time.

This paper is arranged as follows. Section 2 introduces the problem statement and the issue addressed. Section 3 presents the proposed methodology where the problem is solved in the framework of RL. Section 4 examines the feasibility of the proposed approach using an academic example. And finally, the conclusion summarizes the significant advances and includes plans for future work.

PROBLEM FORMULATION

The degradation of a system components affects the performance and the stability of the system. The state of degradation or deterioration, considered as a health indicator, as well affects directly the remaining useful life of the active system, consequently reducing the usability and the productivity of the system. Moreover, the state of health (SoH) is stimulated by the states of the system, and implicitly affected by the action of the controller. Hence the importance of developing an optimal approach for performing a control action that takes into account the performance requirements, the stability and also the SoH of the system. This paper focuses on discrete-time nonlinear systems affine in control as:

x k+1 = f (x k , d k ) + g(x k , d k )u k
(1) where x k ∈ R n is the state vector at time k, u k ∈ R m is the vector of control input, and d k ∈ R l is the vector of degradation's state. f (., .) ∈ R n and g(., .) ∈ R n×m are differentiable and f (0, .) = 0. The degradation is described by:

d k+1 = z(x k , d k )
(2) where d k ∈ R l is the vector of degradation's state and z(., .) ∈ R l is a differentiable function. The evolution of the degradation is implicitly controlled by the input since it is sensitive to the system's states. Some practical examples of the presented formulation are the degradation of the damping coefficient in electromechanical actuator due to bearing wear [START_REF] Fu | Modelling and simulation of flight control electromechanical actuators with special focus on model architecting, multidisciplinary effects and power flows[END_REF]) and the breakdown of the winding insulation caused by the rise of the winding temperature [START_REF] Brown | Prognostics enhanced reconfigurable control of electro-mechanical actuators[END_REF]).

First, one defines the augmented states vector n+l) and condenses the evolution equation into:

X k = [x k , d k ] T ∈ R (
X k+1 = f (X k ) z(X k ) + g(X k ) 0 u k = F (X k ) + G(X k )u k (3)
In this work, the control policy is defined as a feedback controller of the system's and degradation's states:

u k = h(x k , d k ) = h(X k ) (4)
In order to maintain the performance of the system while minimizing the energy and the speed of evolution of degradation, a quadratic utility function is defined by:

r(x k , u k , ∆d k ) = x T k Qx k + ∆d T k Q 1 ∆d k + u T k Ru k (5) where ∆d k = d k+1 -d k is the rate of evolution of degradation. Q, Q 1 and
R are symmetric positive definite weighting matrices and |R|̸ = 0. Thus, the reward can be written as:

r(X k , u k ) = X T k Q 0 0 0 X k + [z(X k ) -[0 1] X k ] T Q 1 [z(X k ) -[0 1] X k ] + u T k Ru k (6)
The matrix Q 1 determines the relative importance or cost associated with the degradation rate term in the overall cost function. By changing the values of Q 1 , the weight of ∆d k can be adjust. Thus, higher values of Q 1 would indicate a higher cost or penalty on the degradation speed.

The notion of optimal behavior is captured by defining a performance index known also as cost function:

V h (X k ) = ∞ i=k γ i-k r(X i , u i ) (7)
with 0 < γ ≤ 1 a discount factor that reduces the weight of costs occurring in the future. Note that it is essential to use a discounted performance function for the proposed formulation. This is because if the rate of degradation does not go to zero, then the performance function ( 7) becomes infinite without the discount factor. Assumption 1 The dynamical system (1) is controllable on some set Ω ⊂ R n which implies that there exists a control law that asymptotically stabilizes the system on Ω.

Assumption 2 The degradation variable d k has a maximum value d max , such that if the degradation level at any time k is less than d max , the system remains stable and can be asymptotically stabilized on the set Ω.

Remark 1. It is noted that rate of degradation is not treated like other state variables of the system owing to the fact the degradation mechanism is irreversible in nature. Moreover, the goal is to slow down the degradation and not to eliminate its speed or stabilise it. As such, the augmented systems does not include the rate of degradation as one of the states to avoid stabilisation of the former to an undesired values including negative ones. Remark 2. Assumption 2 is based on the premise that the effect of degradation on the system dynamics is sufficiently small and it allows us to design a control strategy that guarantees stability and performance for the system under degradation.

To this end, the objective is to find an optimal control input that minimizes the speed of irreversible degradation and extends the operational residual life of the system, while taking into account various factors such as system dynamics and performance criteria.

PROPOSED METHODOLOGY

DP approach determines optimal control solutions using Bellman's principle backward-in-time from some desired goal states. This yields offline solution algorithms which are stored and then implemented online, forward-in-time.

Using RL, the solution of Hamilton-Jacobi-Bellman equation can be solved-for online in real time. This section presents a novel algorithm based on RL under model based framework, leading to the design of degradation tolerant optimal control.

Formulation of optimal control problem

The objective of optimal control theory is to find the policy that minimizes the cost function leading to the optimal value:

V * (X k ) = min h(.) ∞ i=k γ i-k r(X i , h(X i )) (8)
Then, the optimal control policy is defined by:

h * (X k ) = argmin h(.) ∞ i=k γ i-k r(X i , h(X i )) (9)
By writing (7) as:

V h (X k ) = r(X k , u k ) + γ ∞ i=k+1 r(X i , h(X i )) (10)
It yields to the Bellman equation:

V h (X k ) = r(X k , u k ) + γV h (X k+1 ) (11
) According to Bellman's principle of optimality for discretetime systems, the following Hamilton-Jacobi-Bellman (HJB) is obtained:

V * h (X k ) = min h(.) (r(X k , u k ) + γV * h (X k+1 )) (12)
with the optimal policy as:

h * (X k ) = argmin h(.) (r(X k , u k ) + γV * h (X k+1 )) (13)
Since the optimal policy at time k is determined by using the optimal policy at time k + 1, Bellman's Principle solves the optimal control problem backwards-in-time. It is the basis for Dynamic Programming (DP) algorithms.

The optimal control problem is solved by finding the solution of HJB equation ( 12) for the value function.

Then, by substituting the solution in (13) the optimal control is obtained. Due to the nonlinear nature of the HJB equation, it is generally difficult or impossible to find its solution.

Solving Optimal Control Problem using RL

As referred previously solving the HJB equation is a challenging problem in optimal control and dynamic programming, moreover it can be computationally expensive.

To this end, RL has proven to be powerful and effective machine learning technique that allows to handle optimal control problems for nonlinear systems where iterative methods are often used to obtain the solution of Bellman equation indirectly. These methods are classified into two main schemes, namely policy iteration (PI) and value iteration (VI). This paper only focuses on the VI algorithm among others, as dedicated in Algorithm 1 [START_REF] Lewis | Reinforcement learning and adaptive dynamic programming for feedback control[END_REF]).

Definition 1 A control policy is defined as admissible with respect to (1) on Ω, if it stabilizes (1) on Ω, and yields to a finite cost function.

Algorithm 1 Value Iteration Algorithm Initialization. Initialize the value function V 0 using ( 14) and h 0 (X k ) with any control policy, not necessarily admissible (see Definition 1).

V 0 = X T 0 P 0 X 0 (14) where P 0 is a positive-definite matrix. Value Update. Update the value using:

V j+1 (X k ) = r(X k , h j (X k )) + γV j (X k+1 ) (15) Policy Improvement. The control policy is improved by:

h j+1 (X k ) = argmin h(.) (r(X k , h(X k )) + γV j+1 (X k+1 )) (16)
VI does not find the value corresponding to the current policy, but simply performs a single iteration to that value using (15) with j the step index of VI. The proof of convergence of the VI algorithm in the general nonlinear DT setting was presented in [START_REF] Al-Tamimi | Discretetime nonlinear hjb solution using approximate dynamic programming: Convergence proof[END_REF]). It is difficult to solve equations ( 15) and ( 16) for nonlinear systems. Fortunately, NN can be used to approximate V j and h j at each iteration. In the following, NN based actorcritic structure is used to implement the VI algorithm, as presented in the next section.

Actor Critic based implementation

VI algorithm can be implemented using actor-critic structure, where an actor component applies an action or control policy to the system, and a critic component assesses the value of that action and the state resulting from it. The learning mechanism maintained by the Actor-Critic structure [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]) has two steps:

• policy evaluation, executed by the critic and performed by observing the results of applying current control policy on the environment; • policy improvement, performed by the actor.

The combination of DP, NN, and actor-critic structure results in the ADP algorithms. This latter includes three NNs: two critics NNs and one actor NN. The structural diagram of the ADP algorithm is shown in Fig. 1.

Fig. 1. The structural diagram of ADP algorithm

The first critic NN approximates the relationship between the value function Vj+1 (X k ) and the states vector X k . While, the second critic NN approximates the relation between the value function Vj (X k+1 ), at the previous iteration j, and the states X k+1 at instant k + 1. The actor NN approximates the relationship between control vector ĥj (X k ) and the states vector X k . In this paper, the model of the system and the model of the degradation are assumed to be available. Otherwise, the dynamics can be approximated by using available input-output data to train a NN [START_REF] Liu | Adaptive dynamic programming with applications in optimal control[END_REF]).

The output of the critic and actor two-layer networks is described by: Vj

(X k ) = W (j)T c Φ(Y (j)T c X k ) (17) ĥj (X k ) = W (j)T a σ(Y (j)T a X k )
(18) where Φ = [Φ 1 , ..., Φ Nc ] is the vector of hidden-layer activation functions of the critic with Φ(.) = sigmoid(.) and N c is the number of hidden-layer neurons. σ = [σ 1 , ..., σ Na ] is the vector of hidden-layer activation functions of the actor with σ(.) = tanh(.) and N a is the number of hiddenlayer neurons. W The developed method has been extended in a tracking control framework also. In the following section, a tracking control design is presented where the dynamics of the system are required.

Optimal Tracking Control Using ADP

The above VI-based ADP approach can be applied to solve infinite horizon optimal tracking control problem. The objective is to find an optimal control policy h * (.) that regulates the state x k to follow a specified trajectory ref k

Algorithm 2 Value Iteration ADP algorithm

Step 1. Initialize the weights of critic and actor neural networks, and the parameters i

a max , i c max , ϵ a , ϵ c , j max , ξ c , ξ a , Q, Q 1 , R. (ξ a > ξ c )
Step 2. Set the iteration index j = 0 and P 0 = βI n = 0.

Step 3. Choose randomly an array of p state vector {x 1 k , x 2 k , ..., x p k }, and p degradation's state vector {d 1 k , d 2 k , ..., d p k } to obtain an array of dimension p × (n + l) . Initialize V 0 using ( 14). Train the critic network using the data until the given accuracy ϵ c , or the maximum number of iterations i c max is reached.

Step 4. If j = 0, Initialise randomly the vector of initial policy { ĥ0 (X 1 k ), ĥ0 (X 2 k ), ... ĥ0 (X p k )} otherwise go to Step.5.

Step 5. Using 15). Train the critic network until the given accuracy ϵ c , or the maximum number of iterations i c max is reached.

{h j (X 1 k ), h j (X 2 k ), ..., h j (X p k )},
{V j+1 (X 1 k ), V j+1 (X 2 k ), ..., V j+1 (X p k )} by (
Step 8. Compute the target of action network {h j+1 (X 1 k ), h j+1 (X 2 k ), ..., h j+1 (X p k )} by using:

h j+1 (X k ) = - γ 2 R -1 G T ∂ Vj+1 (X k+1 ) ∂X k+1 = - γ 2 R -1 g T ∂ Vj+1 (X k+1 ) ∂x k+1 (19)
Step 9. If j = j max or |V j+1 (X s k ) -V j (X s k )|≤ ξ a with s = {1, 2, ..., p}, train the actor NN until the given accuracy ϵ a , or the maximum number of iterations i a max is reached; otherwise go to Step 10.

Step 10. Set the iteration index j

= j + 1. If j > j max or |V j+1 (X s k ) -V j (X s k )|≤ ξ c , go to
Step 11; otherwise, go to Step 5.

Step 11. Compute the output of the actor NN { ĥj (X 1 k ), ĥj (X 2 k ), ..., ĥj (X p k )}. Obtain the final near optimal control law u * (.) = ĥj (.), and stop the algorithm.

∈ R n . Consider that there exists a feedback control u b,k , known by desired control [START_REF] Liu | Adaptive dynamic programming with applications in optimal control[END_REF]), satisfying the following equation :

u b,k = g + (ref k )[ref k+1 -f (ref k )] (20) where g + (ref k ) = (g(ref k ) T g(ref k )) -1 g(ref k ) T is the generalized inverse of g(ref k ) with g(ref k ) + g(ref k ) = I. By posing ϵ k = x k -ref k and µ k = u k -u b,k , the utility function becomes: r(ϵ k , µ k , ∆d k ) = ϵ T k Qϵ k + µ T k Rµ k + ∆d T k Q 1 ∆d k (21)
And the objective becomes to find µ k that minimizes the cost function leading to the optimal value:

V * (ϵ k , d k ) = min µ k ∞ i=k γ i-k r(ϵ i , µ i , ∆d i ) (22)
To find the solution of the above equation, the iterative ADP algorithm will iterate between value function update (23) and the policy improvement (24):

V j+1 (ϵ k , d k ) = r(ϵ k , h j (ϵ k , d k ), ∆d k ) + γV j (ϵ k+1 , d k+1 ) (23) h j+1 (ϵ k , d k ) = argmin µ k (r(ϵ k , µ k , ∆d k )+γV j+1 (ϵ k+1 , d k+1 ))
(24) The ADP algorithm described above in equations (22-24) is essentially the same as in Algorithm 2. In the following, the NNs trained for regulation are used for tracking but the main difference is that the inputs of NNs (x k , d k ) are replaced by (ϵ k , d k ) and the output of the actor becomes µ k . To verify the effectiveness of the developed control schemes, an infinite horizon tracker is implemented on an academic example in the next section.

SIMULATION RESULTS AND DISCUSSION

Consider the following unstable nonlinear system with affine form:

x k+1 = 0.4x 2,k 0.3x 1,k + x 2,k + 0.1d k x 2,k + 0 1 1 0 u k
The dynamic of evolution of degradation is described by the following equation:

d k+1 = 0.01x 2 2,k |d k |+d k
The prior knowledge of system dynamics is assumed to be available. The weighting matrices are chosen as:

Q = 0.1 0 0 0.1 and R = 2 0 0 2
The three-layers NNs has been chosen as critic network and actor network with the structures of 3-10-1 and 3-10-2 respectively, and γ = 0.98. The initial weights of NNs are initialized to zero. A random array of data has been generated to train the NNs, the state variables are generated in the range of [-0.7, 0.7] and the degradation state in the range of [0.1, 0.43]. The number of data generated is p = 2000. For each iteration j, the critic NN is trained for i c max = 4000 steps under the learning rate 1 × 10 -4 so that the approximation error limit 10 -6 is reached. The actor NN is trained for i a max = 6000 steps under the learning rate 1 × 10 -5 , only at the final iteration j max or when the value function converges to the optimal value. After implementing the outer-loop iteration for j max = 20 times, the convergence of the value function is observed. If V 0 (x k ) ≤ V 1 (x k ) holds for all x k , the value function sequence V j is a monotonically increasing sequence, V j+1 (x k ) ≥ V j (x k ), ∀x k , ∀j ≥ 0 [START_REF] Liu | Adaptive dynamic programming with applications in optimal control[END_REF]). In our case, V 0 (x k ) is fixed to zero. The 3-D plot of approximate value function at j = [1,3,20] is given in Fig. 2, this latter shows that V 20 > V 3 > V 1 . The convergence process of the value function is given in Fig. 3 which also affirms the monotonic behaviour of V j . For x 0 = [0, 0] and d 0 = 1 × 10 -1 , the corresponding state trajectories are displayed in Fig. 4. The curve in green represent the trajectory of the states without taking into consideration the rate of evolution of degradation in the reward function, and the NNs are with the structures of 2-10-1 and 2-10-2, i.e., only the states of the system are used to train the NNs. In this case, the states track the reference with a null steady state error, and the value of the degradation at the final instant of time is equal to 0.4208 (Fig. 5). In order to decelerate the speed and the final value of degradation, ∆d k is integrated in the cost function, and the state of degradation d k is considered as one of the NNs inputs, in addition to the states of the system. For Q 1 = 2 × 10 4 and Q 1 = 2.6 × 10 4 , the red and blue curves are respectively obtained. From Fig. 4, it can be seen that by including ∆d k and d k in the learning and the structure of the NNs, the steady state error increases progressively by raising the value of Q 1 in order to decrease the speed of degradation. Table 1 presents the mean square error (MSE) between the x 1 and the reference for the three cases (No delta, Q 1 = 2 × 10 4 and Q 1 = 2.6 × 10 4 ) where each period is 100 time steps.:

M SE = 1 N N k=1 ||x k -ref k || 2
The Table 1 shows that for Case 1, MSE decreases progressively with respect to time which implies that the controller is focusing only on the performance of the systems ignoring the degradation. While in Case 2 and 3, the MSE increases progressively with respect to time, which confirms the results obtained in Fig. 4 and shows that the controller is trying to compromise between the performance and the speed of degradation. By increasing the value of Q 1 , the controller will prioritize reducing the rate of evolution of degradation over the performance of the system (21). 1.5134 × 10 -9 2.3758 × 10 -6 6.7253 × 10 -5 N 3

1.0128 × 10 -9 1.0256 × 10 -5 8.2514 × 10 -5 N 4

7.1224 × 10 -10 2.7914 × 10 -5 1.042 × 10 -4 N 5

3.1450 × 10 -10 6.1742 × 10 -5 1.3519 × 10 -4 N 6

4.9005 × 10 -11 1.2128 × 10 -4 1.7958 × 10 -4

Table 1. Mean Squared Error between x 1 and ref k per period N cases. At k = 600, the value of degradation for Q 1 = 2×10 4 is higher than for Q 1 = 2.6 × 10 4 , and both values are less than the value obtained when the degradation was not considered in the cost function and the learning.

CONCLUSION

This paper proposes a learning approach towards the design of degradation tolerant control, for nonlinear discrete time systems affine in control. The problem is formulated over an infinite horizon in the framework of reinforcement learning. The proposed methodology is effective in learning control law that leads to reduction in degradation rate of the system whilst maintaining the closed loop stability under the assumptions of admissibility. The approach was examined for three cases. The first case is when the rate of evolution of degradation is not included in the reward function and also the degradation data are not considered as input of neural networks. In the second and third cases, the degradation is integrated in the cost function and also supposed as one of the NN's inputs in addition to the state of the system. By tuning the weighting coefficient matrix Q 1 , the rate of evolution of degradation can be decreased thus preventing the system's breakdown. Future works will focus on implementing the developed algorithm in real time, so the NNs are trained online based on the outputs of the systems without need of the system's model.

  NN weights. The weights of the critic are obtained from NN training during the j-th iteration, and the weights of the actor are trained only at the final iteration. The Gradient Descent algorithm is used to tune the weights of NNs. A summary of the VI-based ADP algorithm for optimal control is provided in Algorithm 2.

Fig. 3 .

 3 Fig. 2. 3-D plot of approximate value function (17) at j = [1, 3, 20] for Q 1 = 2 × 10 4

Fig. 5

 5 Fig. 4. Evolution of x 2 with respect to time

  Compute the target of the critic network training

	compute the
	output of the system {x 1 k+1 , x 2 k+1 , ..., x p k+1 } and {d 1 k+1 , d 2 k+1 , ..., d p k+1 }. Step 6. Compute the output of the critic network
	{ Vj (X 1 k+1 ), Vj (X 2 k+1 ), ..., Vj (X p k+1 )}. Step 7.


domain: info math

text:
Julián Tachella
email: julian.tachella@ens-lyon.fr

Dongdong Chen
email: d.chen@ed.ac.uk

Mike Davies
email: mike.davies@ed.ac.uk

Unsupervised Learning to Solve Inverse Problems: Application to Single-Pixel Imaging

Ces dernières années, les approches basées sur l'apprentissage profond ont obtenu des performances de pointe pour solutionner de nombreux problèmes d'imagerie inverse allant de l'imagerie médicale à la photographie computationnelle. Ces méthodes nécessitent généralement des paires de signaux et mesures pour l'apprentissage. Cependant, pour divers problèmes d'imagerie, nous avons généralement accès qu'à des mesures comprimées des signaux sous-jacents, rendant l'approche basée sur l'apprentissage compliquée, voire impossible car les observations compressées ne contiennent pas d'information dans l'espace nul de l'opérateur de détection directe. Le récent cadre d'imagerie équivariante surmonte cette limitation en exploitant l'invariance aux transformations (translations, rotations, etc.) présentes dans les signaux naturels. Dans cet article, nous exploitons ce nouveau cadre d'apprentissage non-supervisé pour reconstruire des données d'une camera mono-pixel à partir de mesures compressées uniquement. Une série d'expériences montre que la méthode proposée a des performances comparables à celles de l'approche supervisée standard.

Introduction

Linear inverse problems consist of reconstructing a signal x ∈ X ⊂ R n from incomplete and noisy measurements y ∈ R m , that is

y = Ax + ϵ (1)
where A ∈ R n×m is the forward sensing operator and ϵ is the noise affecting the measurements. This is generally an ill-posed task due to the incomplete forward operator A with m < n (or m = n with a large condition number) and the noise affecting the measurements. In order to reconstruct the signals, we need knowledge about the signal distribution and its support X . Classical approaches assume a signal distribution using some prior knowledge about the underlying signals x. For example, the well-known total variation model [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is built on the prior belief that natural images are approximately piecewise smooth. This strategy often yields a loose description of the true model, providing biased and/or suboptimal reconstructions. In recent years, this approach has been replaced by learning the reconstruction mapping y → x directly from data.

Despite the appeal and better performance of the learning-based approach [START_REF] Ongie | Deep learning techniques for inverse problems in imaging[END_REF], in many sensing applications we can only access incomplete measurements y, resulting in a chicken-andegg problem : in order to reconstruct x we require knowledge about the signal model, but to learn this model we require ground truth training data x. Moreover, if the measurement process A is incomplete, it is fundamentally impossible to learn the signal distribution through only measurements y, as there is no information about the set of signals X in the nullspace of A.

The recent equivariant imaging (EI) framework [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF][START_REF] Chen | Robust equivariant imaging : a fully unsupervised framework for learning to image from noisy and partial measurements[END_REF], showed that this fundamental limitation can be overcome by exploiting the invariance of typical signal sets to transformations, such as translation, rotation or scaling.

In this paper, we extend the EI framework to compressive single-pixel imaging [START_REF] Marco F Duarte | Single-pixel imaging via compressive sampling[END_REF][START_REF] Gibson | Padgett. Single-pixel imaging 12 years on : a review[END_REF]. This modality relies on a single detector and is used in applications where standard cameras with arrays of detectors are very expensive or impossible to build, such as sensing in the non-visible spectrum or ultrafast imaging [START_REF] Gibson | Padgett. Single-pixel imaging 12 years on : a review[END_REF]. In these applications, ground-truth signals x may be expensive or impossible to obtain, thus learning from singlepixel measurements alone is very important. The proposed me-thod only requires compressed measurements for training and obtains a similar performance to fully supervised approaches.

Preliminaries

We begin with some basic definitions. As the number of measurements is lower than the dimension of the signal space, i.e., m < n, the operator A has a non trivial nullspace, denoted by N A ⊂ R n . The pseudo-inverse of A is denoted by A † ∈ R n×m . The learning task consists of estimating the parameters θ ∈ R p of a reconstruction function f θ : R m → R n such that f θ (y) ≈ x for all x ∈ X , by minimising a training loss L : R p → R :

arg min θ L(θ).
(

) 2
The reconstruction mapping f θ is modeled as a deep neural network where θ are the weights and biases. We now discuss the choice of L(θ) in the context of supervised and unsupervised learning.

Supervised Learning Standard learning methods learn the parameters θ using a dataset composed of N pairs of signals and associated measurements (x i , y i ), by minimising the following supervised loss

L Sup (θ) = N i=1 ∥x i -f θ (y i )∥ 2 . (3)
However, in many real-world settings, such as medical and astronomical imaging, obtaining ground-truth signals x i can be very expensive or even impossible.

Unsupervised Learning

In settings where we have access to only measurement data y i , it is still possible to train the reconstruction function by enforcing measurement consistency, i.e.,

L MC (θ) = N i=1 ∥y i -Af θ (y i )∥ 2 . (4)
Unfortunately, this approach is ill-fated if A has a non-trivial nullspace, as there are multiple possible reconstruction functions that satisfy measurement consistency :

Proposition 1 (Proposition 1 in [3]). Any reconstruction func- tion f θ (y) : R m → R n of the form f θ (y) = A † y + v(y) (5)
where v(y) : R m → N A is any function whose image belongs to the nullspace of A verifies the measurement consistency requirement.

In other words, the measurement consistency loss doesn't contain any information about X in N A . We require some additional prior information about X in order to learn from measurement data alone, as we will show next.

Equivariant Imaging

Recently, the EI framework [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF] showed that invariance to transformations, such as translations, rotations or reflections, can be enough to learn from measurement data alone. A signal set X is invariant to a group of invertible transformations T 1 , . . . , T |G| ∈ R n×n , if for all x ∈ X , then T g x also belongs to X for all g = 1, . . . , |G|. Most natural signals present invariance to a certain group of transformations. For example, sets of natural 2D images are generally assumed to be invariant to shifts, rotations and/or reflections.

Under the invariance assumption, we have

y = Ax = AT g T -1 g x = A g x ′ (6)
for g = 1, . . . , |G|, where A g = AT g and x ′ = T -1 g x belongs to the signal set. Thus, the invariance of X provides implicit access to the operators A g with potentially different nullspaces, allowing us to learn beyond the nullspace of A. As illustrated in Figure 1, due to the invariance of X , the composition of the forward operator and reconstruction function f θ •A must be equivariant to the group of transformations, that is

f θ (AT g x) = T g f θ (Ax) (7)
for all x ∈ X and g = 1, . . . , |G|. The equivariant constraint on f θ • A can be enforced using the following (unsupervised) loss [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF] :

L EQ (θ) = N i=1 1 |G| |G| g=1 ∥f θ (AT g xi ) -T g xi ∥ 2 (8)
where xi = f θ (y i ) for all i = 1, . . . , N . Taking into account both measurement consistency and equivariance of f θ • A, the equivariant imaging (EI) loss [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF] is given by L EI (θ) = L MC (θ)+ αL EQ (θ), where α is a hyperparameter which controls the tradeoff between equivariance and measurement consistency.

Analysis Under the assumption of invariance and noiseless measurements, i.e., y = Ax, we have

E y {L MC (θ)} = E x { 1 |G| g∈G ∥AT g x -Af θ (AT g x)∥ 2 }. ( 9
)
Using the equivariance constraint in [START_REF] Tachella | Sampling theorems for unsupervised learning in linear inverse problems[END_REF] we get

E y {L MC (θ)} = E x { 1 |G| g∈G ∥AT g x -AT g f θ (Ax)∥ 2 } (10) = E (x,y) {∥M (x -f θ (y)) ∥ 2 } (11)
where

M = 1 |G|    AT 1 . . . AT |G|    . ( 12
)
It is easy to see that the standard supervised loss in ( 3) is obtained by setting M as the identity matrix. If M has rank smaller than n, the unsupervised loss does not penalize reconstruction error in the nullspace of M . Thus, a necessary condition for unsupervised learning is that M ∈ R m|G|×n is of rank n :

Proposition 2 (Theorem 1 in [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF]). A necessary condition for identifying the signal model X from compressed observations is that the matrix in (12) is of rank n.

A detailed analysis of necessary and sufficient conditions for model identification can be found in [START_REF] Tachella | Sampling theorems for unsupervised learning in linear inverse problems[END_REF].

Handling Noise with SURE

In most sensing applications, the observed measurements are corrupted by noise. The performance of the standard EI method degrades with the amount of noise affecting the observations, as the the measurement consistency loss L MC (θ) does not prevent f θ from overfitting the noise. However, if the noise distribution is known, we can correct the measurement consistency loss using Stein's unbiased risk estimator (SURE) [START_REF] Chen | Robust equivariant imaging : a fully unsupervised framework for learning to image from noisy and partial measurements[END_REF], such that the corrected loss is an unbiased estimator of the oracle measurement consistency with pairs of noiseless and noisy measurements (u i , y i ), i.e.,

E y {L MC-SURE (θ)} = E (u,y) { 1 m ∥u -Af θ (y)∥ 2 }
where u i := Ax i are the noiseless measurements. If the noise is iid Gaussian, the SURE loss is given by

L MC-SURE (θ) = N i=1 1 m ∥y i -Af θ (y i )∥ 2 -σ 2 + 2σ 2 m ∇ • h θ (y i ) (13
) where h θ = A • f θ and ∇• is the divergence operator. For certain f θ , computing the divergence can be challenging, so we instead approximate it using the Monte Carlo approach introduced by Ramani et al. [START_REF] Ramani | Monte-Carlo SURE : A black-box optimization of regularization parameters for general denoising algorithms[END_REF], which only requires an additional evaluation of f θ .

The SURE method can be applied to a large family of noise distributions : for example, expressions for Poisson and Poisson-Gaussian noise distributions can be found in [START_REF] Chen | Robust equivariant imaging : a fully unsupervised framework for learning to image from noisy and partial measurements[END_REF]. The robust equivariant imaging (REI) loss is obtained by replacing L MC (θ) by the SURE-corrected one in L EI (θ), that is,

L REI (θ) = L MC-SURE (θ) + αL EQ (θ). ( 14
)
4 Single-Pixel Imaging

The single-pixel camera [START_REF] Marco F Duarte | Single-pixel imaging via compressive sampling[END_REF] is an imaging modality inspired by compressed sensing which uses a single sensor to measure incoming light modulated by spatial masks. These masks are programmable binary patterns which are set using spatial light modulators, digital micromirror devices or LED arrays. Each measurement obtained by the detector corresponds to the scene filtered by a binary pattern. As measuring more patterns requires longer acquisition times and more storage, generally only a few patterns m ≪ n are used to reconstruct the scene. Single-pixel cameras are particularly useful in settings where standard detector arrays are prohibitively expensive such as imaging non-visible spectra and ultrafast imaging. See [START_REF] Gibson | Padgett. Single-pixel imaging 12 years on : a review[END_REF] for a recent survey about this modality.

There are multiple ways to select the patterns : popular choices are iid random entries or Hadamard basis vectors [START_REF] Gibson | Padgett. Single-pixel imaging 12 years on : a review[END_REF]. Here we focus on patterns with random iid entries, such that the entries of the sensing operator are given by

A i,j = 1 √ n
with probability 0.5

-1 √ n with probability 0.5 (15)
for i = 1, . . . , m and j = 1, . . . , n. Due to the iid patterns, the resulting forward operator is highly incoherent with most groups of transformations, thus m ≥ n/|G| measurements are generally enough to verify the necessary condition in Proposition 2. For example, consider images of 28 × 28 pixels and the group of 2D shifts (|G| = 784). The necessary condition is fulfilled with probability 0.94 with just m = 1 pattern, and probability 1 with m ≥ 2 patterns (estimated using 100 random realizations of A).

Experiments

We evaluate the proposed method using the '0' and '1' digits of the MNIST dataset. The dataset consists of 1000 images for training, and 200 images for testing. The images have 28 × 28 pixels, i.e., n = 784, and are normalised to have intensity values in [0, 1]. The single-pixel measurements are obtained using patterns sampled according to (15) and corrupted by iid Gaussian noise with standard deviation σ = 0.05. We compare the following learning methods :

-Pseudo-inverse (A † y) : Linear reconstruction by applying the pseudo-inverse to the observed measurements y i = A † x i (baseline of no learning).

-Measurement Consistency (MC) : Training from only compressed measurements y i using L MC (θ) in (4).

-Robust Equivariant Imaging (REI) : Training from only compressed measurements y i using L REI (θ) in ( 14) with the group of 2D shifts 1 .

-Supervised (Sup) : Standard training using ground truth pairs (x i , y i ) using L Sup (θ) in (3).

1. Note that the set of images of digits is approximately shift invariant. We model f θ as the U-Net network in [START_REF] Chen | Equivariant imaging : Learning beyond the range space[END_REF] and use the Adam optimizer for all learning approaches. In all cases, each gradient step is performed using mini-batches of 4 images, and 3 transformations per mini-batch in the case of REI (chosen randomly from the set of |G| transformations). For the REI approach we use α = 5. The experiment was performed for 3 different sensing matrices with 78, 156 and 235 patterns respectively. 1 shows the average and standard deviation of the test peak signal-to-noise ratio (PSNR) obtained by the different methods. Figure 2 shows reconstructed test images. The unsupervised MC training approach obtains the same performance as the linear pseudo-inverse A † y, effectively failing to learn the signal set. As discussed in Section 2, enforcing only measurement consistency is not enough for learning the reconstruction function. On the other hand, the proposed unsupervised REI approach obtains a 10 dB improvement over the linear pseudoinverse, with an average PSNR which is approximately 2.5 dB behind the fully supervised approach. Moreover, the REI approach obtains good reconstructions with as few as m = 78 measurements per image.

Conclusions and Future Work

We present a novel unsupervised learning approach for singlepixel imaging, where the reconstruction mapping is learned from measurement data alone. The proposed method can reconstructing single-pixel measurement data by relying only on weak invariance properties, and without imposing strong assumptions on the signal distribution. The unsupervised approach is especially important for scientific imaging applications where the goal is to recover information with the smallest amount of prior assumptions about the underlying signals.

We leave the implementation of the proposed method in a real single-pixel imaging system for future work.

FIGURE 1 -

 1 FIGURE 1 -Equivariant imaging systems. If the set of signals is invariant to a certain set of transformations, the comof imaging operator A with the reconstruction function f θ should be equivariant to these transformations.

FIGURE 2 -

 2 FIGURE 2 -Reconstructed test images with m = 78, 156, 235 patterns for the evaluated learning approaches. Networks trained with measurement consistency only fail to provide good reconstructions, performing similarly to the linear pseudo-inverse. By enforcing equivariance, the proposed unsupervised method obtains reconstructions close to the fully supervised setting.

TABLE 1 -

 1 Average test PSNR in dB obtained by the evaluated training methods for different number of patterns m.Table

	m	A † y	MC	REI	Sup
	78 10.5 ± 2.2 10.5 ± 2.2 20.3 ± 2.9 23.3 ± 3.8
	156 10.9 ± 2.3 10.9 ± 2.3 23.0 ± 3.0 25.3 ± 3.6
	235 11.5 ± 2.2 11.5 ± 2.2 24.4 ± 2.9 26.7 ± 3.4


Nanterre's methodology workshop: a very brief summary

A year has passed since we initiated the monthly methodology workshop at Paris Nanterre University. To summarize this first year, here is a brief summary of some topics we went through (by order):

1) Are exploratory hypotheses more compelling than confirmatory hypotheses? (https://www.tandfonline.com/doi/pdf/10.1080/09515089.2022.2113771?needA ccess=true&role=button)

We started the year with a provocative paper from Mark Rubin and Chris Donkin (2022).

Although confirmatory hypotheses are intuitively thought to be more reliable (because less prone to bias, hidden flexibility, and because their results cannot be involved in their rationales, a crucial point for the use novelty principle), Rubin and Donkin suggest that it is valid to use the result from one statistical test to help to create another statistical hypothesis/model as long as the test statistic value from the first test is independent from the test statistic value for the second test (voir figure ci-dessous). They thus criticized the usefulness of pre-registration arguing that it does not prevent researchers from running unplanned tests or post-hoc hypotheses.They also argue that exploratory tests present several benefits from decreasing researchers' biases (commitment, prophecy), data fraud and give reviewers room to make additional contributions to the analysis.

Although we admitted the usefulness of exploratory hypotheses as a way to test statistical assumptions of the data being analyzed (e.g., normality, homogeneity of variance, etc.), calibrate free parameters and explore possibly new lines of researchers, we also concluded there were some limitations. First, one of the limits of post-hoc testing is the difficulty to control for the inflation of Type-I error. Post-hoc analyses are also usually influenced by some of the research results making them more dependent on previous testing than hypotheticodeductive models. It also makes prior quantification extremely difficult.

Figure from Rubin & Donkin (2022)

2) There is nothing magical about Bayesian statistics (https://www.tandfonline.com/doi/abs/10.1080/01973533.2020.1792297)

Approaching Christmas, we had a nice surprise with a brief introduction to Bayesian statistics and some manipulation on R/Jamovi. Swiatkowski and Carrier (2020) write up a very smooth introduction to the topic. They first recall the difference between frequentist and bayesian approaches. Following their example, if we imagine that a virus called statisticus malignum (which produces compulsive reading of statistical books and R manipulations) is present in 1 out of 1,000 individuals, the frequentist probability of having the virus is simply the ratio between the number of actual occurrences / the total number of possibles occurrences which, in the long run, give us = 1/1000 = 0.001. But after thinking it through, you realize you are a potential nerd because of your interests in maths and R and so re-evaluate your risk of having the virus at 5% (it is just a subjective belief compared to the objective frequency we calculated before).

As the authors mentioned, you could consult Dr. Golmorokov, an expert on this virus to know whether you carry it or not. He considered two hypotheses: H0: "The student does not carry the virus" and H1: "The student does carry the virus" with 99% of statistical power (β = .01) and a Type I error rate of 10% (α = .10). The test came out positive and the Dr. logically prescribed you the cure. Thus, in the long run, the question is "what should I do"? Should I reject or not the null while controlling for the error rate. Note that anything that can impact α level should thus be taken into account (e.g., multiple comparisons, post-hoc testing, etc. In social psychology, there is a tendency to force the response which could introduce other biases (we will discuss them in the 4) part). A more ecological approach is to give the participant/patient the option not to answer the question. But then you will have to decide how to handle the missing data.

Assuming the data are missing completely at random (MCAR), you could potentially eliminate participants who did not answer a question if their frequency is very low (< 10% of the total sample). It is still important to check that participants you're eliminating do not differ substantially from the remaining sample on characteristics that could impact the reliability of the results.

The first classical technique is multiple imputation (MI). According to Little et al. (2013), "Imputing missing data involves making a copy of the original data set and replacing missing values with plausible estimates of what those values would have been, had they been observed" (p.152). For instance, mean imputation consists of replacing the missing value by the mean of all other cases which preserve the sample mean for that variable. However, the SEs is thus reduced which introduces bias to the associated significant tests. Stochastic regression imputation helps with this problem by adding some variability in these predicted values. If we know the observed mean and standard deviation can help impute values much more realistic compared to deterministic regression imputation (where all values fall down on the regression line or very close to). You can see a very nice R tutorial here: https://statisticsglobe.com/regression-imputation-stochastic-vs-deterministic/ Graphic based on https://statisticsglobe.com/regression-imputation-stochastic-vsdeterministic/

The paper mentioned other techniques such as Full-Information Maximum Likelihood when missing data in SEM/multilevel modeling. Both techniques seem to yield unbiased estimates of parameters and SEs when data are missing at random (Enders, 2010 ;Schafer & Graham, 2002). Multi-form designs can also be implemented, either in cross-sectional or longitudinal, by making them complete a subset of the form or the complete form but on a subset of occasions respectively.

4) Methods for the detection of carelessly invalid responses in survey data☆ (https://www.sciencedirect.com/science/article/pii/S0022103115000931?via%3Dihub)

Our last session was subsequent to the workshop on missing data but this time dealing with careless data (which is much more common in social psychology, where answering to all items is usually mandatory). The author, Paul Curran, is detailing a comprehensive list of the different techniques used to detect aberrant or careless data.

One of the most used is the response time: how long does it take for a participant to answer a set of items. A very slow responder can be reliable whereas a very fast responder can exhibit a careless pattern. The author showed it is no easy task to remove outliers based on response time even where the respondent is clearly faster than the average sample (d = 3). Huang et al.

(2012) suggested a cut score of 2 seconds per item. Another technique is the long-string analysis that involves examining the longest string of identical responses from each participant. A lazy participant will usually use a pattern of similar response to save some time while doing the survey (such as responding a lot of 3,3,3,3…). Author suggests removing participants with a string of consistent responses equal or greater than half the length of the total scale. For instance, in a 1 to 5 a Likert scale, a pattern of response such as [3; 3; 4; 4; 4; 3; 5; 3; 3] exhibit a string of 4 three times in a row which is > to 2.5 (half the scale) and should thus be removed. The author mentioned a number of different techniques (table below) but warning we should always consider Type I and Type II errors before removing anyone from a sample.

Figure

  Figure from Swiatkowski & Carrier (2020)

  Our Nanterre's workshop group agreed on the need to use a set of different indices to exclude potential outliers such as median absolute deviation on response time for instance (Leys et al., 2013), some residual pattern (standardized residual or robust Mahalanobis distance, see Leys et al., 2018) and attention check items.

  ). However, from a bayesian perspective, you must integrate your subjective evaluation of 5% of chance you are carrying the virus. Given the power of Dr. Golmorokov's test, there is a 99% chance of the test being positive if you carry the virus. On the other hand, there is only a 10% chance of the test being positive if you do not carry the virus. We can now calculate the predictive positive value (PPV) so you can see summarization by the authors on the figure down below. It represents the ratio of true positive results of the tests divided by the total of positive results (including the false positive) : .99 (99% of detecting the virus if you carry it) x .05 (5% your prior of having the virus) / (.99*.05 + .10*.95) (the probability of your carrying the virus + the false positive).

  Table from Curran (2016)