title: Model Order Reduction for (Stochastic-) Delay Equations With Error Bounds
authors: Simon Becker Lorenz Richter
text: 
  
  
  We analyze a structure-preserving model order reduction technique for delay and stochastic delay equations based on the balanced truncation method and provide a system theoretic interpretation.
  Transferring the framework of <cit.>, we find error estimates for the difference between the dynamics of the full and reduced model. This analysis also yields new error bounds for bilinear systems and stochastic systems with multiplicative noise and non-zero initial states.
  
  
  Vehicle Telematics Via Exteroceptive Sensors:
   A Survey
      Fernando Molano Ortiz, 
          Matteo Sammarco, 
          Luís Henrique M. K. Costa, Senior Member, IEEE
          and Marcin Detyniecki
  F. M. Ortiz and L. H. M. K. Costa are with the GTA/PEE-COPPE/DEL-Poli, Federal University of Rio de Janeiro, Brazil, (e-mail: {fmolano, luish}@gta.ufrj.br)
  M. Sammarco and M. Detyniecki are with Axa, France (e-mail: {matteo.sammarco, marcin.detyniecki}@axa.com).
      December 30, 2023
  =================================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  In this article we study a delay-structure preserving model order reduction method, first discussed for finite-dimensional bilinear delay systems in <cit.>, based on the bilinear balanced truncation technique, for deterministic delay systems of the following types[our analysis immediately extends to N φ^del(t-τ) replaced by the sum ∑_i=1^n N_i φ^del(t-τ_i).]
   
   
  
      φ^del'(t)    = A φ^del(t) +  N φ^del(t-τ)+ Bu(t),  for  t ∈ (0,T)
         φ^del(0)   =φ_0,  φ^del(t)=f(t)  for  t ∈ (-τ,0)     and 
         φ^bild'(t)    = A φ^bild(t) +  N φ^bild(t-τ)v(t)+ Bu(t),  for t ∈ (0,T)
         φ^del(0)   =φ_0,  φ^del(t)=f(t)  for  t ∈ (-τ,0).
  
  
  on arbitrary (separable) Hilbert spaces K for time-dependent control functions u,v ∈ L^2((0,∞)) and a delay parameter τ>0. In particular, for zero delay, τ=0, (<ref>) reduces to the form of a standard bilinear system and thus our error bound provides also a new error bound for the important class of bilinear systems with non-zero initial conditions, extending the analysis of <cit.>. We also discuss different types of discrete and continuous delay, cf. (<ref>).
  
  
  
  Furthermore, we adapt our analysis to stochastic differential equations with delay (SDDE) on finite-dimensional system spaces K ≃ℝ^d
  
      dX_t    = (A X_t +Bu(t))  dt  + ∑_i=1^k N_i X_t-τ_i dW_t^i,      Y_t =CX_t 
      
      X_0   =ξ,  X_t=f_t  for  t ∈ (-τ,0),
  
  where (W_t^i)_i=1,…,k are i.i.d. copies of standard Brownian motion.
  
  
  
  Balanced truncation is a well-established model order reduction technique, especially for linear ((<ref>) with N=0 <cit.>), bilinear ((<ref>) with τ=0 and v=u <cit.>), and stochastic differential equations with multiplicative noise ((<ref>) with τ_i=0 <cit.>). For bilinear with multiplicative noise, error bounds have only been obtained for zero initial condition so far - which we aim to overcome with this article, too. For stochastic systems with multiplicative noise, this has been only addressed in the recent work <cit.>. The bilinear balanced truncation method identifies, as has been rigorously proven at least for linear systems, a subspace of jointly reachable and controllable states and aims to preserve these subspaces well under model order reduction <cit.>. The two properties are approximately captured by positive-definite operators that are called Gramians. The dominant eigenspaces of the product of two Gramians are then used to define the reduced order model. 
  The method we consider is inspired by linear, bilinear, and stochastic balanced truncation theory and preserves the delay-structure of the original system. 
  
  In this article we study delay equations (<ref>), (<ref>) using bilinear balanced truncation, and the SDDE (<ref>) using stochastic balanced truncation. We show that the Gramians still have a system theoretic interpretation and also derive error bounds for the reduced order model. 
  
  We also emphasize that balanced truncation is a method particularly designed for dissipative systems. For all balanced truncation methods, one therefore assumes that the operator A generates an exponentially stable and strongly continuous semigroup. Moreover, (<ref>) is less well-adapted to the dissipative structure, since – without further assumptions on N – the delay term can create growing modes in the dynamics. This is prevented in (<ref>) by assuming that the control v is small.
  
  To take non-zero initial states into account, we consider for equations (<ref>) and (<ref>) a space Y:=span{φ_i;  i ∈{1,…,N}} of admissible input states, where (φ_i) is an orthonormal system in K and introduce a map B_inu = ∑_i=1^n u_i φ_i that we include in the model order reduction.
  
  
  Before stating the main results of this article, let us mention other model order reduction techniques for delay systems. By rewriting the delay equation as a linear equation, using the head-tail representation <cit.> on an infinite-dimensional space, a reduction method based on linear balanced truncation theory has been proposed in <cit.>. Apart from that article, applying the method of balanced truncation to delay equations, other model order reduction techniques such as rational approximation methods <cit.>, interpolation methods <cit.>, Krylov space methods <cit.>, moment matching based methods <cit.> have been proposed.
  
  
  Many of the above methods, however, do not fully preserve the delay structure, which is fixed by the method proposed in this article.
  
  
  
   §.§ Outline of the article
  
  
  
    * In Section <ref> we discuss the algorithmic aspects of the balanced truncation method that we propose for delay systems. 
  
  
    * Section <ref> provides an overview over the C_0-semigroup approach to deterministic delay equations. 
  
  
    * In Section <ref>, we introduce the balanced truncation theory framework including the central objects of the theory, i.e. the bilinear Gramians and the Hankel operator for the deterministic equations. We then show in Propositions <ref> and <ref> that the bilinear Gramians preserve the structure of the equation also for delay equations.
  
  
    * The proof of the error bounds for bilinear and deterministic delay equations with non-zero initial data, stated in Theorem <ref>, is given in Section <ref>.
  
  
    * In Section <ref> we then treat the stochastic delay systems. 
  
  
    * Finally, in Section <ref>, we consider delay models from physics including Stuart-Landau oscillators and generalized Langevin equations, and apply the balanced truncation method for delay systems.
  
  
  
  
   §.§ Main results
  
  
  To state the error estimates, let ‖ N ‖, ‖ C ‖, and ‖ B ‖ denote the maximum of the respective norm for the full and reduced system and M,ω>0 such that for both systems the corresponding semigroup satisfies ‖ T(t) ‖≤ Me^-ω t. We also introduce the norm
  
      ‖ u ‖_L^2 ∨∞ :=max{‖ u ‖_L^∞((0,t),ℝ^n), ‖ u ‖_L^2((0,t),ℝ^n)}.
   
   
   
  Our main error estimate for deterministic systems is stated in the following theorem and applies to bilinear systems as well:
  
  [Error bound deterministic systems]
  
  Let ℋ≃ℝ^m be the output space and consider the difference of two solutions to (<ref>). We assume that both solutions satisfy the stability condition M‖ N ‖ /√(2ω)<1 and ‖ v ‖_L^2(0,T)≤ 1, such that the Volterra series converges <cit.>. Then, for control functions u,v ∈ H^1((0,T),ℝ), initial states φ_0 = ∑_i=1^k ⟨w , e_i⟩ϕ_i and φ_0:= ∑_i=1^k ⟨w , e_i⟩ϕ_i, and zero history function, it follows for H^bil being the Hankel operator associated to the system that Δ(Cφ^bild), which is the difference of the output of the full and reduced system, satisfies  
  
      ‖Δ(Cφ^bild) ‖_L^2((0,∞), ℝ^m)   ≤ 4   ‖Δ(H^bil)‖_TC(‖φ_0 ‖_Kmax{ 1, ‖ v ‖_L^∞(0,T)}
              +max{‖ u ‖_L^2(0,T), ‖ v ‖_L^1(0,T)}‖ u ‖_L^∞(0,T)).
  
  
  The proof of this result is given in Section <ref>.
  
  The corresponding error bound for stochastic (delay) differential equations with multiplicative noise is stated in the following theorem: 
  [SDDE error bound]
  
  For control functions u ∈ L^∞_ωL^2_t(Ω_(0,T),ℝ^n), zero history functions, i.e. f_t≡ 0, and two solutions to (<ref>) with initial conditions ξ := ∑_i=1^k ⟨v , e_i⟩ξ_i with L^2(Ω,ℱ_0,K) orthonormal system (ξ_i), and ξ:= ∑_i=1^k ⟨v , e_i⟩ξ_i, it follows that for two independent Wiener processes, where both solutions satisfy M‖ N ‖ /√(2ω)<1, and  H^sdde being the Hankel operator associated to the system
  
      ‖Δ(CX^sdde) ‖_L^2(Ω_(0,T),ℝ^m) ≤‖Δ(H^sdde) ‖_TC(‖ξ‖_L^2(Ω;K) + 2 ‖ u ‖_L^∞_ωL^2_t(Ω_(0,T),ℝ^n))
  
  
  
  The proof of this result is given at the end of Section <ref>.
  
  
   §.§ Notation
  
  The space of bounded linear operators between Banach spaces X,Y is denoted by ℒ(X,Y) and just by ℒ(X) if X=Y. The operator norm of a bounded operator T ∈ℒ(X,Y) is written as ‖ T ‖. The space of trace class operators is denoted by TC(X,Y) and the space of Hilbert-Schmidt operators by HS(X,Y).
   In particular, we recall that for a linear operator T ∈TC(X,Y), where X and Y are now separable Hilbert spaces, the trace norm is given as
  
      ‖ T ‖_TC=sup{∑_n ∈ℕ|⟨ f_n,T e_n  ⟩_Y | : (e_n)  ONB of  X  and  (f_n)  ONB of  Y }.
  
  
  
  The resolvent set of an operator A is denoted by ρ(A) and we say that g=𝒪 (f) if there is a C>0 such that ‖ g ‖≤ C ‖ f ‖. In order not to specify the constant C, we also write ‖ g ‖≲‖ f ‖.
  The domain of unbounded operators A is denoted by D(A).
  
  Let H be a separable Hilbert space. To define the Hankel operator we require a decomposition of the positive Gramians. For this purpose, we introduce the Fock space F( H) of H-valued functions F( H):=⊕_k=0^∞F_k( H) where F_k( H):=L^2((0,∞)^k, H) and F_0( H):= H. Thus, elements of the Fock space are sequences of F_k-valued elements.
  
  We introduce function spaces L^1_iL^2_k-1 and ℋ^∞_i ℋ^2_k-1 norms which for H-valued functions functions f:(0,∞)^k → H and g: ℂ_+^k → H are defined by
  
      ‖ f ‖_L^1_iL^2_k-1( H) = ∫_0^∞‖ f(∙,…,∙,s_i, ∙,…,∙) ‖_L^2((0,∞)^k-1, H) ds_i  and 
         ‖ g ‖_ℋ^∞_i ℋ^2_k-1( H)= sup_s_i ∈ℂ_+‖ g(∙,…,∙,s_i, ∙,…,∙) ‖_ℋ^2((0,∞)^k-1, H).
  
  For k-variable functions h we sometimes also write h^(i)(s,t):=h(s_1,…,s_i-1,t,s_i,…,s_k-1) in order to shorten the notation. We denote by H^1 the Sobolev space of L^2 functions whose first weak derivative is in L^2 as well.
  
  
  
  
  § BALANCED TRUNCATION FOR DELAY SYSTEMS IN A NUTSHELL
  
  
  In this section, we provide a brief overview over the model order reduction method studied in this article and the computability of the error bounds we stated in Theorems <ref> and <ref>.
  
  To obtain the reduced order model for delay systems on finite-dimensional system spaces K of the form (<ref>), or (<ref>) with (possibly) multiple delays, we compute positive (semi-)definite observability 𝒪 and reachability 𝒫 Gramians from the following Lyapunov equations
  
      A𝒫 + 𝒫A^T + ∑_i=1^k N_i 𝒫 N_i^T + BB^T + B_inB_in^T =0  and 
         A^T𝒪 + 𝒪 A + ∑_i=1^k N_i^T 𝒪 N_i + C^TC = 0.
  
  Since both 𝒪 and 𝒫 are positive semidefinite, they can be decomposed as 𝒪 = W^*W and 𝒫 =RR^*.
  Let 𝒪 and 𝒫 have both full rank for simplicity, then the balanced representation is obtained by performing a singular value decomposition (SVD) WR = V Σ U^* and introducing operators Q:=Σ^-1/2 V^*W and Q^-1:=RUΣ^-1/2 such that the new balanced matrices are given by
  
      A_b := QAQ^-1,  N_b:=QNQ^-1,  B_b:=QB,  B_in b:=QB_in,   and C_b:=CQ^-1.
  
  To obtain a reduced system, the smallest singular values of the matrix Σ are discarded. 
  The error bounds stated in this article are given in terms of the trace distance of the difference of certain Hankel operators H, introduced in Def. <ref>, which we denote by Δ(H), for the full and reduced system.
  
  To actually compute the singular values of Δ(H), and thus the trace norm of Δ(H), one defines an error system 
  
      A:=( A     0 
       0    Ã),  B := (B,B̃)^T,  B_in := (B_in,B̃_in)^T,C:=(C,-C̃),  and N:=(N     0 
       0   Ñ),
  
  where operators without tilde belong to the original system in (<ref>), which we call System 1  and operators with tilde correspond to a second system, that we call System 2. This could be any other system with the same structure as the reduced system.
  Then one can define Gramians 𝒪=W^*W and 𝒫=RR^* of this error system (<ref>) that also satisfy Lyapunov equations
  
      A𝒫 + 𝒫A^T + ∑_i=1^k N_i 𝒫N_i^T + BB^T + B_inB_in^T =0  and 
         A^T 𝒪 + 𝒪A + ∑_i=1^k N_i^T 𝒪N_i + C^TC = 0.
  
   We can then perform a singular value decomposition WR = VΛU^* with diagonal operator Λ that contains all singular values of the error system (<ref>) on its diagonal <cit.>. Hence, we find ‖Δ(H) ‖_TC = ∑_λ∈Λλ, as there exist unitary mappings <cit.> U: ran(WR) →ran(H) and V: ker^⊥(WR) →ker^⊥(H)
  such that 
  
      Δ(H) |_ker^⊥(H) = H|_ker^⊥(H)=U(WR) |_ker^⊥(WR) V^*|_ker^⊥(H).
  
  
  
  
  
  
  
  
  § DETERMINISTIC DELAY EQUATIONS
  
  
  For our analysis of delay equations, we start with a C_0-semigroup approach to study well-posedness, Volterra kernel expansions, that will be essential in the proof of the error bounds, and stability properties. 
  
  
  
   §.§ Delay equations in a semigroup framework
  
  
  We start by introducing a delay operator Φ∈ℒ (H^1((-r,0),K),K), where H^1 is the standard Sobolev space. By normalizing appropriately, we can always assume that the history interval, i.e. the interval by which the dynamics reaches back in time, is (-r,0) for some fixed r>0. To fully define the dynamics, we require an initial value φ_0 ∈ K and a history function f∈ L^2((-r,0),K).
  
  Then, for a control function u ∈ L^2((0,∞), ℝ^n) and bounded operators B ∈ℒ(ℝ^n, K) of the form Bu = ∑_i=1^n ψ_i u_i for some ψ_i ∈ K, we study delay equations of type (<ref>)
  
      φ^del'(t)    = Aφ^del(t)+ Φ( φ^del)(t) + Bu(t),  t>0 
      φ^del(0)   = φ_0,   φ^del(σ)= f(σ),  for σ∈ (-r,0).
   
  
  We consider two types of delays operators Φ: H^1((-r,0),K) → K for our error analysis. 
  [Delay operators] 
  
  For bounded operators N ∈ℒ(K), and two types of delays 
      
    * constant time delays r>τ>0 or
      
    * integral delays g ∈ L^1(-r,0) such that ∫_-r^0 | g(s) | ds ≤ 1,
  
  we introduce delay operators Φ_d and Φ_c for f ∈ H^1((-r,0),K)
  
      Φ_df  :=  Nf(-τ)  or Φ_cf:= ∫_-r^0 Nf(s) g(s)  ds.
  
   
   
  
  
  Furthermore, we introduce a vectorized control-to-state map 
  (control operator) 𝐁, (initial-state operator) 𝐁_in, and a state-to output map (observation operator) 𝐂 such that on 𝒦= K × L^2((-r,0), K)
  
      𝐁:=(B,0)  ∈ℒ(ℝ^n,𝒦),   𝐁_in:=(B_in,0)    ∈ℒ(ℝ^k,𝒦),
          and 𝐂(φ,f):=Cφ with 𝐂∈ℒ(𝒦,ℋ).
  
  
  To apply the theory of C_0-semigroups, we introduce operators 𝒜_0 and 𝒩 on 𝒦
  
      𝒜_0:=( A      0 
       0    d/d σ)  and 𝒩:=( 0     Φ
       0     0 ) ∈ℒ(D(𝒜_0),𝒦)
  
  with domain D(𝒜_0):={ (φ,f) ∈ D(A) × H^1((-r,0), K): f(0)=φ}⊂𝒦[The linear space D(𝒜_0) carries the graph norm ‖ x ‖_D(𝒜_0) = ‖𝒜_0 x ‖+ ‖ x ‖.] where A is itself a generator of a strongly continuous semigroup (T(t)) on K.  
  We then define an operator 
  
      S_t: K → L^2((-r,0),K) by (S_tφ)(-τ):=_[0,∞](t-τ) T(t-τ)φ
  
  and T_←(t) ∈ℒ(L^2((-r,0),K)) as the nilpotent left-shift semigroup 
  
      (T_←(t)φ)(-τ)=(_[-r,0]φ)(t-τ).
  
  It can be shown <cit.> that the operator 𝒜_0 is the generator of a C_0-semigroup 
  
      𝒯_0(t):= (  T(t)     0 
       S_t     T_←(t) ).
  
   We continue with a stability condition that plays the analogous role to exponential stability of the semigroup for delay equations: 
   [MV-condition]
  We assume the semigroup 𝒯_0 in (<ref>) and operator 𝒩 to satisfy a Miyadera-Voigt L^2-condition for some q<1, i.e.
  
      ∫_0^∞‖𝒩𝒯_0(r) φ‖^2  dr ≤ q ‖φ‖^2.
  
   
  
  
  
  
  We say that 𝒯_0 and 𝒩 satisfy a (truncated) Miyadera-Voigt L^1-condition if for some q<1 and t_0>0:
  
      ∫_0^t_0‖𝒩𝒯_0(r) φ‖ dr ≤ q ‖φ‖.
  
  By Hölder's inequality, (<ref>) implies (<ref>) for some t_0>0.
  
  
  Under Assumption <ref>, the Miyadera-Voigt perturbation theorem <cit.> implies that the state operator 𝒜:=𝒜_0+𝒩 is the generator of a semigroup (𝒯(t) ) and the delay equation (<ref>) is well-posed <cit.>. That is, the solution to (<ref>) is continuous and we can write the solution to (<ref>) by Duhamel's formula for x_0=(φ_0,f) ∈𝒦 and u ∈ L^2((0,∞), ℝ^n) as
  
      Z(t)   =𝒯(t)x_0+ ∫_0^t 𝒯(t-s) 𝐁u(s)  ds  with output Y(t)=𝐂(Z(t)).
  
  Let π_1 be the projection from 𝒦∋ (x,f) ↦ x, then π_1(Z) solves (<ref>).
  
  
  Many previous ideas related to model order reduction methods for delay utilized the linear structure of the representation (<ref>). The computational issue with this approach is that the system is inherently infinite-dimensional.
  
  On the other hand, when studying delay system using linear balanced truncation, there are explicit criteria for exponential stability of the semigroup (𝒯(t)):
  <cit.>
  
  If the semigroup (T(t)), generated by A in the delay equation (<ref>), satisfies ‖ T(t) ‖≤ M e^-ω t with ω strictly larger than some α≥ 0 and Me^ατ/ω-α‖ N ‖<1,
  then the decay bound ω_0(𝒜) of the delayed semigroup 𝒯(t) satisfies ω_0(𝒜)>α≥ 0, i.e. there is some ℳ>0 such that ‖𝒯(t) ‖≤ℳ e^-ω_0(𝒜) t.
  
  
  
  
  
   §.§ Volterra series expansion of the dynamics
  
  
  We record that Definition (<ref>) implies the existence of an operator 𝒮∈ℒ( 𝒦, L^2((0,∞),𝒦)) extending 𝒩𝒯_0(t) from D(𝒜_0) to all of 𝒦 with norm ‖𝒮‖≤γ<1. The operator 𝒮, that extends 𝒩𝒯_0(t), satisfies 𝒮(s)𝒯_0(t)= 𝒮(t+s). Let D(𝒜_0) ∋ x_n → x ∈𝒦, then 
  
      𝒮(s)𝒯_0(t)x= lim_n →∞𝒮(s)𝒯_0(t)x_n = lim_n →∞𝒩𝒯_0(s)𝒯_0(t)x_n =𝒮(s+t)x.
  
  
  
  
  The Miyadera-Voigt perturbation theorem <cit.> allows us to express the semigroup (𝒯(t)) in ℒ(K) generated by the state operator 𝒜=𝒜_0+𝒩 as a series 𝒯(t)x = ∑_i=0^∞𝒱^i  𝒯_0(t)x.
  
  The operators 𝒱^i are the so-called Volterra operators 𝒱∈ℒ(L^p([0,T],𝒦)) defined, for any p ∈ [1,∞], as 
  
      (𝒱F)(t)x = ∫_0^t  F(t-s)  (𝒮x)(s)  ds = ( F *(𝒮x  _[0,∙] ))(t),  for all t  ∈ [0,1]
  
  such that by Young's inequality ‖𝒱 F ‖_L^p≤‖𝒮‖_ℒ(L^1,𝒦)‖ F ‖_L^p.
  
  
  
  From (<ref>) and Fubini's theorem it follows that for sets with i ∈ℕ 
   Δ_i(t):={s ∈ℝ^i; 0≤ s_i≤…≤ s_1≤ t } and delay Volterra kernels  
  
      h^delay_i(t)   :=𝒪_i-1(t)𝐁∈ℒ(ℝ^n,ℋ)  and  h^delay_i,in(t):=𝒪_i-1(t)𝐁_in∈ℒ(ℝ^k,ℋ)
          where 𝒪_i(t=(t_1,…,t_i+1))y:=𝒞𝒯(t_1)  ∏_l=2^i+1(𝒮(t_l) )y
  
  the solution Y, in (<ref>), for initial conditions x_0:=(φ_0,0), and u ∈ L^p((0,∞),ℝ^n) is a function Y ∈ L^p_loc((0,∞), ℋ) given by Y(t) = K_1(t)+K_2(t) where
  
      K_1(t)= ∑_i=0^∞𝒞𝒱^i  𝒯_0(t)x    and    K_2(t)   =  ∫_0^t 𝐂𝒯(t-s)𝐁u(s)  ds
         =  ∑_i=1^∞∫_Δ_i(t) h^delay_i(t-s_1,…,s_i-1-s_i)  u(s_i) ds
  
  where 𝒞 is defined in (<ref>).
  In (<ref>) and (<ref>), we introduced the delay Volterra kernels, which we shall write down more explicitly, for delay types (<ref>) indicated by indices d(iscrete)|c(ontinuous), in the following Lemma:
  
  
  The delay Volterra kernels (<ref>) satisfy for any j ∈ℕ
  
      h^delay_d,j| d,j, in(t_1,…,t_j)   := CT(t_1)∏_i=2^j(_(0,∞)(t_i-τ)N T(t_i-τ))B_|in
      
      h^delay_c,j | c,j, in(t_1,…,t_j)   := CT(t_1)∏_i=2^j(∫_-r^0 _(0,∞)(t_i+s)N T(t_i+s) g(s) ds )B_|in.
  
   
  
  We restrict our proof to discrete delays Φ=Φ_d in the proof. Consider (φ,f) ∈ D(𝒜_0), then an explicit computation shows that
  
      𝒮(t)(φ,f)= Φ(S_t φ +T_←(t) f) =_(0,∞)(t-τ) N T(t-τ)φ + _[0,τ](t)  N f(t-τ).
  
  Let (φ_i,f_i) ∈ D(𝒜_0) converge to (φ,f) ∈𝒦, then by continuity of 𝒮(t)
  
      𝒮(t)(φ,f)    = lim_i →∞𝒮(t)(φ_i,f_i)= lim_i →∞𝒩𝒯_0(t)(φ_i,f_i)=lim_i →∞ (Φ( S_tφ_i+ T_←(t)f_i),0) 
         = lim_i →∞(_(τ,∞)(t) N T(t-τ) φ_i  +_[0,τ](t) N f_i(t-τ), 0 ) 
         =     (_(0,∞)(t-τ) N T(t-τ) φ +  _[0,τ](t) N f(t-τ),0)
  
  shows that 𝒮(t)(φ,0)=( _(0,∞)(t-τ) N T(t-τ) φ,0 ).
  From the definition of the Volterra kernels h^delay_m (<ref>) we thus conclude that the first equation in (<ref>) holds. 
  
  
  
  
  § BILINEAR BALANCED TRUNCATION FOR DELAY SYSTEMS
  
  
  In this section, we provide the necessary tools from bilinear balanced truncation theory and apply it to delay systems.
  
  Let ℋ be a separable Hilbert space and C ∈ℒ(K,ℋ) the observation operator, we then introduce the central object of the bilinear balanced truncation theory, the bilinear Gramians, cf. <cit.>. In particular, the mapping and regularity properties have all been shown in <cit.>.
  [Gramians]
   
  Let O_0(t_1):=CT(t_1). Then, for i ∈ℕ and y ∈ K define O_i(t_1,…,t_i+1)y:=CT(t_1)  ∏_j=2^i+1(N T(t_j) )y and bounded operators 𝒪_i for x,y ∈ K by ⟨ x, 𝒪_iy ⟩_K := ∫_(0,∞)^i+1⟨ O_i(s)x, O_i(s)y ⟩_ℋ ds,
  which are summable in operator norm. The bilinear observability Gramian 𝒪^bil∈ℒ(K) is then given as 𝒪^bil := ∑_i=0^∞𝒪_i ∈ℒ(K).
  
  For the bilinear reachability Gramian, we define for i ∈ℕ and y ∈ K   
      P_i(t_1,…,t_i+1)y:= ∏_j=1^i(T(t_j)^*N^*)T(t_i+1)^*y
  
  such that
  
      ⟨ x,𝒫_0y ⟩_K
         :=∫_(0,∞)⟨ T(s)^*x, BB^* T(s)^*y ⟩_K  ds +⟨ x,B_inB_in^* y ⟩ and for i ∈ℕ,
      ⟨ x,𝒫_iy ⟩_K
         :=∫_(0,∞)^i+1⟨ P_i(s)x, BB^*P_i(s)y ⟩_K ds + ∫_(0,∞)^i⟨ P_i-1(s)x, B_inB_in^* P_i-1(s)y ⟩_K ds.
  
  The bilinear reachability Gramian is the operator 𝒫^bil:= ∑_i=0^∞𝒫_i ∈TC(K). 
  
  We introduce operators W^bil and R^bil such that the observability Gramian is 𝒪^bil= W^bil*  W^bil and the reachability Gramian is 𝒫^bil =  R^bil  R^bil*. 
  [Observ. & reach. map]
  For i ∈ℕ_0 we define the family W_i∈ℒ( K, F_i+1(ℋ)) of operators that map K ∋ x ↦ O_i(∙)x such that ‖  W_i ‖ = 𝒪( (M ‖ N ‖ ( 2 ω )^-1/2)^i). Then, we can define the bilinear observability map W^bil∈ℒ(K, F(ℋ)) by W^bil(x):=( W_i(x) )_i ∈ℕ_0.
  
  Let R_i ∈HS(F_i+1(ℝ^n) ⊕ F_i(ℝ^k),K) be given by  
  
      R_0(f,g):= ∫_(0,∞) T(s) B f(s)  ds+  B_ing   and for  i ∈ℕ
         R_i(f,g):=∫_(0,∞)^i+1 P_i(s)^*Bf(s)  ds + ∫_(0,∞)^i P_i-1(s)^*B_ing(s)  ds.
  
  The bilinear reachability map is defined as R^bil∈HS( F(ℝ^n) ⊕⊕_i=0^∞F_i(ℝ^k),K) such that (f_i,g_i)_i ∈ℕ_0↦∑_i=0^∞  R_i( f_i,g_i).
  For subsequent use, we also define maps 
  
      F_i(f_i):=R_i(f_i,0) and G_i(g_i):=R_i(0,g_i).
  
  
  Using the above two operators R^bil and W^bil, we can now introduce the bilinear Hankel operator, cf. <cit.>.
  [Hankel operator]
  
  The Hankel operator is the Hilbert-Schmidt operator H^bil:= W^bil  R^bil∈HS( F(ℝ^n)⊕⊕_i=0^∞F_i(ℝ^k),F(ℋ)).
  
  In particular, if ℋ is finite-dimensional then H^bil is of trace-class.
  
  To relate the delayed dynamics to the bilinear Gramians we introduce the integral kernels of the bilinear Hankel operator:
  
  
  The bilinear Volterra kernels h_m and h_m,in are the functions defined for m ∈ℕ_0 by
  
      h_m(t_0,…,t_m)   :=O_m(t_0,…,t_m)B and 
      h_m,in(t_0,…,t_m):=O_m(t_0,…,t_m)B_in.
  
  
  From (<ref>), we have estimates on the difference of the trace distance of two Hankel operators with F_i,G_i as in (<ref>),
  
      ∑_i=0^∞‖Δ(W_iF_i) ‖_TC≤ 2‖Δ(H^bil) ‖_TC,  ∑_i=0^∞‖Δ(W_iF_i+1) ‖_TC≤ 2‖Δ(H^bil) ‖_TC and 
         ∑_i=0^∞‖Δ(W_iG_i) ‖_TC≤ 2‖Δ(H^bil) ‖_TC,  ∑_i=0^∞‖Δ(W_iG_i+1) ‖_TC≤ 2‖Δ(H^bil) ‖_TC.
  
  The bilinear Gramians satisfy the following Lyapunov equations which hold as operator equations, i.e. without testing against elements x,y ∈ D(A), as soon as A is a bounded operator, see (<ref>),
  [Lyapunov equation]
  For all x,y ∈ D(A) the bilinear observability Gramian satisfies the Lyapunov equation
  
      ⟨𝒪^bilAx,y ⟩_K + ⟨𝒪^bilx,Ay ⟩_K + ⟨𝒪^bil Nx,Ny ⟩_K+  ⟨ Cx,Cy⟩_K=0.
  
  For all x,y ∈ D(A^*) the bilinear reachability Gramian satisfies the Lyapunov equation
  
      ⟨(𝒫^bil-B_inB_in^*)A^*x,y ⟩_K + ⟨(𝒫^bil-B_inB_in^*)x,A^*y ⟩_K 
         + ⟨(𝒫^bil-B_inB_in^*) N^*x,N^*y ⟩_K
       +  ⟨ (BB^*+B_inB_in^*)x,y⟩_K=0.
  
  
  
  Both Lyapunov equations can be obtained immediately from results on bilinear theory <cit.> as follows: We observe that 𝒫^bil can be written as a sum of two standard bilinear reachability Gramians 𝒫^bil_a| b with control operators B_a=B and B_b=B_in, respectively, and the projection B_inB_in^* added to them.
  Each of the Gramians 𝒫^bil_a| b satisfies a Lyapunov equation <cit.> with i∈{a,b}
  
      ⟨𝒫_i^bil A^*x,y ⟩_K + ⟨𝒫^bil_i x,A^*y ⟩_K + ⟨𝒫^bil_i N^*x,N^*y ⟩_K
       +  ⟨ B_i B_i^*x,y⟩_K=0.
  
  Adding them together and using that 𝒫^bil - B_inB_in^* = 𝒫_a^bil + 𝒫_b^bil yields the claim. The Lyapunov equation for the observability Gramian coincides with the one in <cit.>.
  
  The Lyapunov equations (<ref>) and (<ref>) imply the following two interpretations for Gramians with different kind of delays: 
  
  
  All elements (φ_0,f) ∈ker(𝒪) × L^2((-r,0),ker(𝒪)) are unobservable under the evolution of the homogeneous system. That is, the solution to 
  
      φ'(t)=Aφ(t)+N(Kφ)(t),  for t>0
   
  satisfies Cφ(t)=0 for all t>0 for τ>0 and g:[-r,0]→ℝ as in (<ref>), where
  
      (Kφ)(t) = φ(t-τ)  or  (Kφ)(t)= ∫_-r^0g(s) φ(t+s)  ds.
  
  
  
  We start by assuming that (φ_0,f) ∈ (ker(𝒪)∩ D(A)) × H^1((-r,0),ker(𝒪)) with φ_0 =f(0) first. This set is dense, as it is precisely D(𝒜) ∩(ker(𝒪) × L^2((-r,0),ker(𝒪))), cf. (<ref>).
  As in the proof of <cit.> one shows that the first Lyapunov equation (<ref>) yields
  
      N (ker(𝒪) ) ⊂ker(𝒪),  ker(𝒪) ⊂ker(C),  and   A(ker(𝒪)∩ D(A) ) ⊂ker(𝒪).
  
  
  Hence, the homogeneous equation (<ref>) for α∈ker(𝒪)^⊥ and φ(t) ∈ker(𝒪) ∩ D(A) satisfies
  
      ⟨α, φ'(t) ⟩ = ⟨α, A φ(t) ⟩ +   ⟨α, N Kφ(t) ⟩=0.
  
  Thus, the flow of the homogeneous problem (<ref>) leaves ker(𝒪) invariant.
  From the inclusion ker(𝒪) ⊂ker(C), we obtain
  C φ(t)=0.
  The statement follows then for arbitrary φ:=(φ_0,f) ∈ker(𝒪) × L^2((-r,0),ker(𝒪)) by approximating φ with D(𝒜_0) ∩ker(𝒪) × L^2((-r,0),ker(𝒪)) such that  ⟨α,π_1(𝒯(t)φ)  ⟩= lim_i →∞⟨α,π_1(𝒯(t)φ_i)  ⟩=0.
  
  
  
  
  The closure of the range of the reachability Gramian 𝒫 is an invariant space of the flow of φ'(t) = Aφ(t)+  NKφ(t) + Bu(t)
  with K as in (<ref>).
  That is, for any (φ_0,f) ∈ran(𝒫)× L^2((-r,0),ran(𝒫)) the solution φ stays in ran(𝒫).
  
  
  We start by assuming that (φ_0,f) ∈ D(𝒜) ∩ (ker(𝒫)^⊥× H^1((-r,0),ker(𝒫)^⊥), first. 
  From the second Lyapunov equation (<ref>) we deduce, similarly to the previous Lemma, that
  
      N (ker(𝒫)^⊥) ⊂ker(𝒫)^⊥,  ψ_i ∈ker(𝒫)^⊥,  and   A( ker(𝒫)^⊥∩ D(A) ) ⊂ker(𝒫)^⊥
  
  where we recall that Bu =∑_i=1^n ψ_i u_i.
  
  This shows that the homogeneous equation for α∈ker(𝒫) and φ(t) ∈ker(𝒫)^⊥∩ D(A) satisfies
  
      ⟨α, φ'(t) ⟩ = ⟨α, A φ(t) ⟩ +  ⟨α, N(K φ)(t) ⟩=0.
  
  Thus, the flow of the homogeneous problem (<ref>) leaves ker(𝒫)^⊥ invariant.
  From ker(𝒫)^⊥∩ D(A) ⊂ker(C), we obtain
  C φ(t)=0.
  The statement follows then for general x:=(φ_0,f)∈ker(𝒫)^⊥× L^2((-r,0),ker(𝒫)^⊥) by approximating x with 
  
      x_i ∈ D(𝒜) ∩ker(𝒫)^⊥× L^2((-r,0),ker(𝒫)^⊥)
   such that  ⟨α,π_1(𝒯(t)x)  ⟩= lim_i →∞⟨α,π_1(𝒯(t)x_i)  ⟩=0. 
  The inhomogeneous equation satisfies then ⟨π_1(Z(t)),α⟩ =⟨π_1(𝒯(t)x),α⟩+ ∑_i=1^n∫_0^t ⟨π_1(𝒯(t-s)(ψ_i,0)),α⟩ u_i(s)  ds=0.
  
  
  
  
  
  § PROOF OF THEOREM <REF>
  
  
  In this section we provide the proof of Theorem <ref> which provides error estimates on the difference of two delay systems with control (<ref>):
  
      φ^bild'(t)    = A φ^bild(t) + (Φφ^bild)(t) v(t)+ Bu(t),  t>0 
      φ^bild(0)   = φ_0,   φ^bild(σ)= 0,  for σ∈ (-r,0).
  
  
  As in (<ref>), we let K_1(t):=Cφ_1^bild(t) be defined in terms of the solution to
  
      φ_1^bild'(t)    = A φ_1^bild(t) + (Φφ_1^bild)(t) v(t)+ Bu(t),  t>0 
      φ_1^bild(0)   =0,   φ_1^bild(σ)= 0,  for σ∈ (-r,0)
  
  and K_2(t):=Cφ_2^bild(t) where φ_2^bild solves
  
      φ_2^bild'(t)    = A φ_2^bild(t) + (Φφ_2^bild)(t) v(t),  t>0 
      φ_2^bild(0)   = φ_0,   φ_2^del(σ)= 0,  for σ∈ (-r,0).
  
  
  
  
  In the following we write v(t):=u(t-τ). By applying the triangle inequality to the Volterra series <cit.>, we have using delay Volterra kernels (<ref>)
  
      ‖Δ(K_1) ‖_L^2(0,T)≤‖‖Δ(_[0,∞)h^delay_0 ) ‖ * _[0,T)‖ u ‖‖_L^2(0,∞) + 
         ∑_k=2^∞(∫_0^T(∫_Δ_k(t)‖Δ ( h^delay_k-1(t-s_1,…,s_k-1-s_k) v(s_1)⋯ v(s_k-1)· u(s_k) ‖ ds )^2  dt )^1/2.
  
  Our aim is to estimate the expression containing delay Volterra kernels in terms of bilinear Volterra kernels (<ref>).
  
  The first term on the right-hand side of (<ref>) satisfies by Young's inequality
  
      ‖‖Δ(_[0,∞)h^delay_0)  ‖ * _[0,T)‖ u ‖‖_L^2(0,∞)   ≤‖Δ(h^delay_0) ‖_L^1(0,∞)‖ u ‖_L^2((0,T))
         ≤‖Δ(h_0) ‖_L^1(0,∞)‖ u ‖_L^2(0,T).
  
  To estimate the second term on the right-hand side, we observe that by Minkowski's inequality and Hölder's inequality for k ≥ 2
  
      (∫_0^T(∫_Δ_k(t)‖Δ (h^delay_k-1(t-s_1,…,s_k-1-s_k) ) v(s_1)⋯ u(s_k) ‖  ds )^2  dt )^1/2
         ≤∫_0^T(∫_s_1^T(∫_Δ_k-1(s_1)‖Δ (h^delay_k-1(t-s_1,…,s_k-1-s_k)) v(s_2)⋯ u(s_k) ‖ ds )^2 dt )^1/2| v(s_1) | ds_1  
         ≤sup_0≤ s_1≤ T(∫_0^T-s_1(∫_Δ_k-1(s_1)‖Δ (h^delay_k-1(t,…,s_k-1-s_k) ) v(s_2)⋯ u(s_k) ‖ ds )^2   dt )^1/2‖ v ‖_L^1(0,T)
         ≤(∫_0^Tsup_0≤ s_1≤ T(∫_Δ_k-1(s_1)‖Δ(h^delay_k-1(t,s_1-s_2, …,s_k-1-s_k) ) v(s_2)⋯ u(s_k) ‖ ds )^2   dt )^1/2‖ v ‖_L^1(0,T).
  
  Thus, by applying Hölder's inequality, to the inner integral we find that using ‖ v ‖_L^2≤ 1
  
      ∫_Δ_k-1(s_1)‖Δ(h^delay_k-1(t,s_1-s_2,…,s_k-1-s_k) )  | v(s_2)⋯ u(s_k) |‖ ds
         =∫_0^s_1…∫_0^s_k-1‖Δ h^delay_k-1(t,s_1-s_2,…,s_k-1-s_k)  v(s_2) ⋯ u(s_p) ‖_ℒ(ℝ^n, ℝ^m) ds_k⋯ ds_2 
         ≤(∫_(0,T)^i-2(∫_0^∞(∫_(0,T)^k-i‖Δ h_k-1^delay(t,s,r,q) ‖^2_ℒ(ℝ^n, ℝ^m) dq )^1/2 dr )^2  ds  )^1/2‖ u ‖_L^∞(0,T)
  
  Applying Minkowski's integral inequality to this expression, leads after a change of variables to remove the delay, together with (<ref>) to
  
      (∫_0^T(∫_Δ_k(t)‖Δ (h_k-1^delay(t-s_1,…,s_k-1-s_k)B ) v(s_1)⋯ u(s_k) ‖  ds )^2  dt )^1/2
         ≤∫_0^∞( ∫_(0,T)^k-1‖Δ h_k-1(q_1,…,q_i-1,r,q_i,…,q_k-1) ‖^2_ℒ(ℝ^n, ℝ^m)  dq  )^1/2 dr  ‖ u ‖_L^∞(0,T)‖ v ‖_L^1(0,T).
  
  Thus, we have together with (<ref>)
  
      ‖Δ(K_1) ‖_L^2(0,T)   ≤∑_i=1^∞(‖Δ(h _2i-1) ‖_L^1_iL^2_2i-1(HS)+ ‖Δ( h _2i-2) ‖_L^1_iL^2_2i(HS))·
             max{‖ u ‖_L^2(0,T)‖ u ‖_L^∞(0,T)‖ u ‖_L^1(0,T)}
         ≤  4 ‖Δ(H^bil) ‖_TCmax{‖ u ‖_L^2(0,T), ‖ u ‖_L^∞(0,T)‖ v ‖_L^1(0,T)}.
  
  
  
  The difference Δ(K_2) can be bounded, in terms of w as introduced in Theorem <ref>, using the Cauchy-Schwarz inequality, and Minkowski's integral inequality as
   
      ‖Δ K_2 ‖_L^2(0,T)   ≤‖Δ(CT B_in)w‖_L^2(0,T)
            + ∑_i=1^∞‖∫_Δ_i(·)‖Δ(h^delay_i,in(·-s_1,⋯,s_i-1-s_i,s_i))(v)  ‖ v(s_1)⋯ v(s_i)  ds ‖_L^2(0,T)
         ≤‖Δ(CT B_in) ‖_L^2(0,T)‖φ_0 ‖_K 
            +  ∑_i=1^∞∫_0^T ‖Δ(h_i,in(⋯,s,⋯))   ‖_L^2((0,∞)^i,ℋ) ds  ‖ v ‖_L^∞(0,T)‖φ_0 ‖_K.
  
  The statement then follows from (<ref>) by inserting the estimate
  
      ‖Δ(CT B_in) ‖_L^2(0,T)≤‖Δ(W_0G_0) ‖_HS≤‖Δ(W_0G_0) ‖_TC
   
  into the estimate on the difference of all Volterra kernels (<ref>), cf. <cit.>,
  
      ‖Δ K_2 ‖_L^2(0,T)   ≤(‖Δ( h _2i,in) ‖_L^2+  ∑_i=1^∞(‖Δ(h _2i-1,in) ‖_L^1_iL^2_2i-1(HS)+ ‖Δ( h _2i,in) ‖_L^1_iL^2_2i-2(HS)) ) 
         ·max{1,‖ v ‖_L^∞(0,T)}‖φ_0 ‖_K ≤  4 ‖Δ(H^bil) ‖_TCmax{1,‖ v ‖_L^∞(0,T)}‖φ_0 ‖_K.
  
  
  
  
  
  We now explain how to extend the previous error bound to systems (<ref>) without a control function v:
  
  
  By redefining 
  
      Φ̂   = √(T_0)Φ,  N̂ = √(T_0)N, and v(t)=T_0^-1/2,
   
  (<ref>) on the time-interval [0,T_0] becomes the uncontrolled delay system
  
      φ^bild'(t)    = A φ^bild(t) + (Φ̂φ^bild)(t) + Bu(t),  t∈ (0,T_0] 
      φ^bild(0)   = φ_0,   φ^bild(σ)= 0,  for σ∈ (-r,0)
  
  with ‖ v ‖_L^2(0,T_0)=1.
  
  
  Hence, the following corollary follows straight from Theorem <ref>
  
  
  
  Let ℋ≃ℝ^m and consider the difference of two solutions to (<ref>) on a time interval [0,T_0]. If we then interpret these two equations as solutions to (<ref>) with redefined (<ref>) and  M‖N̂‖ /√(2ω)<1, such that the Volterra series converges <cit.>, then, for control functions u ∈ L^2((0,T_0),ℝ^n), initial states φ_0 = ∑_i=1^k ⟨w , e_i⟩ϕ_i and φ_0:= ∑_i=1^k ⟨w , e_i⟩ϕ_i, and zero history function, it follows that 
  
      ‖Δ(Cφ^bild) ‖_L^2((0,∞), ℝ^m)   ≤ 4   ‖Δ(H^bil) ‖_TC(‖φ_0 ‖_Xmax{ 1, T_0^-1/2}
              +max{‖ u ‖_L^2(0,T_0), √(T_0)}‖ u ‖_L^∞(0,T_0)).
  
  
   
  
  
  § STOCHASTIC DELAY DIFFERENTIAL EQUATIONS
  
  
  The balanced truncation theories of bilinear and stochastic systems with multiplicative noise have many features in common <cit.>. In particular, the Gramians for both systems obey the same Lyapunov equations. We now demonstrate that the same is true for stochastic delay equations (SDDEs): Consider a probability space (Ω, ℱ, ℙ) with filtration ℱ_t induced by a one-dimensional Brownian motion (W(t))_t. We then study the stochastic delay differential equation (SDDE) for a process X^sdde_t taking values in ℝ^d with matrices A,N ∈ℝ^d × d, B ∈ℝ^d × n and a control u ∈ L^2_ad(Ω× (0,∞); ℝ^n) which, cf. (<ref>), is given by
  
      dX^sdde_t    = (A X^sdde_t+ Bu(t))  dt +N X^sdde_t-τ dW_t,  t>0 
       
      X^sdde_0   =ξ,  X^sdde_t = f_t  for -r≤ t < 0.
  
  where f_t is a ℱ_0-measurable C([-r,0],ℝ)-valued random variable for which the second moment exists uniformly on [-r,0] and ξ∈ L^2(Ω,ℱ_0). The solution to (<ref>) satisfies then
  
      X^sdde_t = T(t)ξ +  ∫_0^t T(t-s)  NX^sdde_s-τ dW_s+ ∫_0^t T(t-s)Bu(s)  ds.
  
  We are interested in an output variable CX_t where the observation operator C is a matrix of appropriate size.
  This identity shows that X_t is a semi-martingale, since X_t is also (ℱ_t) adapted, has continuous paths for t ≥ 0, and satisfies for all T>0 the uniform square-integrability condition 𝔼(sup_t ∈ (-r,T)‖ X_t ‖^2 ) < ∞. Moreover, X_t is a C^b-Feller process <cit.>.
  In particular, we can define for any C([-r,0])-valued process Ψ and initial condition ξ∈ L^2(Ω,ℱ) as in (<ref>) the flow Φ^sdde_Ψ(t)(ξ):=X_t^hom where X_t^hom is the homogeneous part of (<ref>), i.e. the process X^sdde with u ≡ 0. 
  We record that once Ψ≡ 0 the flow becomes linear in the initial state ξ. We also write Φ^sdde_Ψ(t,s)(ξ) to denote the process started at time s from ξ.
  It follows then directly from (<ref>) that for Ψ=0, the solution to (<ref>) is given by the variation of constant formula 
  
      X_t^sdde =  ∫_0^t Φ^sdde_0(t,r)   Bu(r) dr+Φ^sdde_0(t)(ξ).
  
  
  We assume the flow Φ^sdde_0 to be exponentially stable in mean square sense, i.e. there are C,κ>0 such that for all x ∈ℝ^m we have 𝔼‖Φ^sdde_0(t,0)x ‖^2 ≤ C e^-κ t‖ x ‖^2. 
  
  [Exponential stability]
  The exponential stability of SDEs and SDDEs has been thoroughly addressed and relevant results for our framework can for example be found in <cit.>. In particular, <cit.> implies that the SDDE
  
      dX_t = A X_t  dt + ∑_i=1^l (A_i X_t + N_i X_t-τ) dW_t^i
  
  possesses an exponentially stable flow if there are two positive definite matrices Q,G such that the Lyapunov equation holds
  
      GA + A^T G + ∑_i=1^l ( A_i + N_i )^T G (A_i+N_i ) = -Q
  
  and the delay satisfies the smallness condition τ < √(δ_1^2+ δ_3 ‖ A ‖^2)-δ_1/2‖ A ‖^2
  where in terms of the smallest eigenvalue λ_min(T) of a matrix T satisfies
  
      δ_1   :=∑_i=1^l (‖ A_i ‖^2 + ‖ N_i ‖^2 ),   δ_2:=2 ‖ G ‖√(2δ_1 ∑_i=1^l ‖ N_i ‖^2) ,  and 
      δ_3   :=(√(δ_2^2+4 λ_min(Q) ‖ G ‖∑_i=1^l ‖ N_i ‖^2 )-δ_2/2‖ G ‖∑_i=1^l‖ N_i ‖^2)^2.
  
  
  
  
   §.§ Gramians
  
  Let us now introduce the Gramians and the Hankel operator for the SDDE (<ref>):
  
  
  The observability map W^sdde∈ℒ( K, L^2(Ω_(0,∞); ℋ)) and reachability map R^sdde∈ℒ(L^2(Ω_(0,∞); ℝ^n) ⊕ℝ^k,K) are defined as 
  
      (W^sddex)(t)= CΦ_0^sdde(t)x  and  R^sdde(f,v)= 𝔼∫_0^∞Φ^sdde_0(s)Bf(s)  ds +B_inv .
  
  The Hankel operator is defined as H^sdde=W^sddeR^sdde.
  The stochastic reachability Gramian is defined as 𝒫^sdde= R^sddeR^sdde* where
  
      𝒫^sdde=𝔼∫_0^∞ (Φ^sdde_0(s)B)(Φ^sdde_0(s)B)^T  ds + B_inB_in^*.
  
  The observability Gramian is defined as 𝒪^sdde= W^sdde*W^sdde where 
  
      𝒪^sdde = 𝔼∫_0^∞Φ^sdde_0(s)^T C^T C Φ^sdde_0(s)  ds.
  
  
  
  
  The observability and reachability Gramians satisfy the following Lyapunov equations
  
      BB^T +(𝒫^sdde-B_inB_in^T) A^T+ A  (𝒫^sdde-B_inB_in^T)+  N(𝒫^sdde-B_inB_in^T)  N^T=0 
          and    C^TC +N^T𝒪^sdde N+𝒪^sdde A+ A^T 𝒪^sdde=0.
  
  
  
  Stochastic integration by parts of the process X_t = Φ^sdde_0(t)ξ brings
  
      𝔼(X_t X_t^T )    = 𝔼(ξξ^T)+∫_0^t 𝔼(X_s X_s^T )  ds  A^T+ A ∫_0^t 𝔼(X_s X_s^T )  ds
            +N ∫_0^t 𝔼(X_s-τ X_s-τ^T )  ds  N^T.
  
  We then take the initial state x_0=ψ_i with ψ_i as in Bu =∑_i=1^n ψ_i u_i, pass to the limit t →∞, and perform a simple change of variables s':=s-τ in the last integral.
  This immediately yields (<ref>).
  Let us now find a Lyapunov equation for the observability Gramian. Stochastic integration by parts shows that 
  
      𝔼(Φ^sdde_0(t) 𝒪^sddeΦ^sdde_0(t) ) = 𝒪^sdde + 𝔼( ∫_0^t Φ^sdde_0(s)^T 𝒪^sdde A Φ^sdde_0(s)  ds )
         + 𝔼( ∫_0^t Φ^sdde_0(s)^T A^T 𝒪^sddeΦ^sdde_0(s)  ds ) + 𝔼( ∫_0^t Φ^sdde_0(s-τ)^T N^T𝒪^sdde NΦ^sdde_0(s-τ) )  ds 
         = 𝒪^sdde + 𝔼( ∫_0^t Φ^sdde_0(s)^T 𝒪^sdde A Φ^sdde_0(s)  ds ) + 𝔼( ∫_0^t Φ^sdde_0(s)^T A^T 𝒪^sddeΦ^sdde_0(s)  ds )  
            + 𝔼( ∫_0^t Φ^sdde_0(s)^T N^T𝒪^sdde NΦ^sdde_0(s) )  ds.
  
  Thus, if X satisfies the Lyapunov equation C^TC +N^T X  N+X A+ A^T X=0, we find by (<ref>) and the definition of the observability Gramian, Definition <ref>, that indeed 𝒪^sdde=X. 
  
  
  
  The proof of the Lyapunov equations, in particular (<ref>) and (<ref>), show that for non-constant history function f_t, the Lyapunov equations generalize to
  
      BB^T+ ∫_-τ^0 N𝔼(f_t  f_t^*)N^T  dt +(𝒫^sdde-B_inB_in^T) A^T+ A  (𝒫^sdde-B_inB_in^T)
         +  N(𝒫^sdde-B_inB_in^T)  N^T=0    and 
         C^TC +N^T𝒪^sdde N+ ∫_-τ^0𝔼(f_t N^T 𝒪^sdde N f_t^*)  dt +𝒪^sdde A+ A^T 𝒪^sdde=0.
  
  
  
  
  
  In addition, the observability Gramian defines the L^2 energy of the uncontrolled process, i.e. with u ≡ 0: 
  
      ℰ_output:=∫_0^∞𝔼‖ CΦ_0^sdde(t,0) x ‖^2  dt = ⟨𝒪^sddex,x ⟩.
  
  For the reachability Gramian less direct interpretations, already studied in the context of stochastic system with multiplicative, cf. <cit.>, can be stated for stochastic delay equations, as studied here, too.
  
  
  
   
  
  The difference of Hankel operators Δ(H^sdde) satisfies for two independent Wiener processes 
      ‖Δ(CΦ^sdde B_in)‖_L^2(Ω_(0,∞),HS(ℝ^k, ℝ^m))   ≤‖Δ(H^sdde) ‖_HS and 
      ‖Δ(CΦ^sdde B) ‖_L^1_tL^2_ω(Ω_(0,∞),HS(ℝ^n, ℝ^m))   ≤ 2‖Δ(H^sdde) ‖_TC.
  
  
  
  The first bound follows immediately from the definition of the Hilbert-Schmidt norm. The second bound follows along the lines of the proof of <cit.>.
  
  
  We are now ready to give the proof of the error bound for stochastic delay differential equations:
  
  
  By Young's inequality we have for v ∈ℝ^k as in the statement of the theorem that ‖ v ‖_ℝ^k= ‖ξ‖_L^2(Ω,K) and 
  
      ‖Δ(CX^sdde) ‖_L^2(Ω_(0,T),ℝ^m)
         (<ref>)≤‖Δ(CΦ^sdde B_in)(v)‖_L^2(Ω_(0,T),ℝ^m) + ‖‖Δ(_[0,∞)C Φ^sdde B ) ‖ * _[0,T)‖ u ‖‖_L^2(Ω_(0,T),ℝ^m)
         ≤‖Δ(C Φ^sdde B_in)‖_L^2(Ω_(0,∞),ℝ^m × k)‖ξ‖_L^2(Ω;K)
            +  ‖Δ(C Φ^sdde B) ‖_L^1_tL^2_ω(Ω_(0,∞),ℝ^m × n)‖ u ‖_L^∞_ωL^2_t(Ω_(0,T),ℝ^n)
         Lemma <ref>≤‖Δ(H^sdde) ‖_TC(‖ξ‖_L^2(Ω;K) + 2‖ u ‖_L^∞_ωL^2_t(Ω_(0,T),ℝ^n)).
   
   
  
  
  
  
  
  
  
  
  § APPLICATIONS AND EXAMPLES
  
  
  
  We conclude by analyzing three applications of the model order reduction methods studied in this article – the corresponding code can be found at <https://github.com/lorenzrichter/balanced-truncation>. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Example:  Linearized Stuart-Landau Oscillator <cit.> 
  The Stuart-Landau oscillator system is a coupled nonlinear network and a model for phenomena such as chaos or synchronization in large physical or biological systems modeling a finite speed of propagation. 
  The dissipative system is described by a parameter (α)<0 and has an equilibrium solution φ=0, i.e. all oscillators at rest. The dynamics of the Stuart-Landau system with unidirectional nearest-neighbor interaction is described by the following coupled system of nonlinear differential equations
  
      φ_j'(t) = αφ_j(t) - φ_j(t) |φ_j(t) |^2 + φ_(j+1) N (t-τ),  j=1,…, N.
  
  with τ>0. 
  By linearizing the dynamics (<ref>) around the equilibrium solution, we obtain the following coupled system of delay equations
  
      φ_j'(t) = αφ_j(t)+ φ_(j+1) N(t-τ),  j=1,…, N.
  
  Hence, A=α and therefore T(t)=e^α t and thus M=1 and ω=|(α) |. In a numerical example, let us consider d = 50 and take α = -1.2, T = 2, τ = 0.1. We discretize the dynamics with a simple forward Euler scheme using the stepsize Δ t = 0.01 and a random initial point φ(0) ∼𝒩(0, √(0.5)). We choose a history function f=0 and do balanced truncation as described above while varying the dimensions of the corresponding reduced systems. Figure <ref> displays the L^2 errors when reducing to r-dimensional systems and compares to the corresponding bounds from Corollary <ref>. We see that the measured L^2 error decays rapidly and that the bound seems to be rather conservative. The two right plots in the same figure show some components of the full and the reduced systems when either choosing r=2 or r=6. In the latter case we see almost full agreement of the full and reduced trajectories.
  
  
  
  
  
  
  Example:  Generalized Langevin equation <cit.>
  Next, we consider a collection of particles obeying a damped Newtonian dynamics
  
      M ẍ(t) + C ẋ(t)  + K x(t) = 0.
  
  In (<ref>) both the coupling matrix K= L_1L_1^* and mass matrix M = L_2L_2^* are assumed to be strictly positive matrices and the friction matrix F is a positive semidefinite matrix with the property that the form x ↦⟨ x,F x ⟩ is non-degenerate on every eigenspace of M^-1K.
  
  Introducing new coordinates y(t) :=(L_1^*x(t),L_2^*ẋ(t)), we can define a generator
  
      A = ( 0     L_1^* L_2^-1*
       -L_2^-1 L_1     -L_2^-1 F L_2^-1*)
  
  such that the semigroup defined as T(t):=e^tA, associated with (<ref>), is exponentially stable. Thus, the solution to (<ref>) is given by y(t)= T(t) y_0.
  To model particle motion in contact with a heat bath the following generalized Langevin equation (GLE) has been proposed
  
      M ẍ(t) + F ẋ(t)  + K x(t) =- ∫_0^t γ(t-s)Γ_0 ẋ(s)  ds-B_0u(t),
  
  where B_0u(t) describes an external fluctuation force. To include a memory effect in the dissipation the friction is perturbed by a convolution between a kernel γΓ_0 and the velocity of the particles, where Γ_0 is a matrix representing the typical scale of friction. 
  In order to cast (<ref>) in a form that resembles more the type of delay equations we have been studying in this article, we use matrices N:= diag(0, Γ_0 ) and B:= diag(0, B_0 ),
  such that the dynamics (<ref>) takes the form
  
      y'(t)    = A y(t) + ∫_0^t  γ(t-s) Ny(s) ds + Bu(t).
  
  We then consider fractional Brownian motion (fBm) B_H with Hurst parameter H ∈ (0,1) and correlation 𝔼(B_t^H B_s^H) = 1/2( s^2H+t^2H-| t-s |^2H).
  Moreover, for t>s=0 the distributional time derivative of fBm has correlation coefficients 𝔼(Ḃ_t^H Ḃ_0^H) = H(2H-1)t^2(H-1) which tend to zero for large t. For Hurst parameters H ∈ (1/2,1) the correlation coefficient is also integrable at zero. We then use a cut-off function to truncate the small tail (neglecting memory effects from time more than a distance r away from the current time) and define the history kernel in (<ref>) to be γ(t) :=_[r,0](t) t^2(H-1). Thus, 
  
      y'(t)    = A y(t) + ∫_t-r^t  γ(t-s) Ny(s) ds + Bu(t)
         ≈ A y(t) +  rγ(r/2) Ny(t-r/2) + Bu(t).
  
  Now, the high-dimensional system (<ref>) can be reduced using balanced truncation as described in this article. For a numerical illustration, let us choose d=100 M =  + 0.1diag(a_1, …, a_d), F =  + 0.1diag(a_1, …, a_d), K =  + (|a_ij/2|)_i,j=1^d, Γ_0 =, where a_i, a_ij∼𝒩(0, 1) are chosen i.i.d. for all i, j ∈{1, …, d}. We further take u(t) = sin(20 t) 1, C = , T = 10, τ = 0.1 and plot components of the full and reduced systems in Figure <ref>, again considering different dimensions r for the reduced models. In spite of a dimension reduction from d = 100 to r = 10 the trajectories look almost the same in the right plot.
  
  
  
  Example:  Geometric Brownian motion
  We finally consider geometric Brownian motion in dimension d = 40 as an example for a stochatsic delay equation just as in (<ref>) with a one-dimensional Brownian motion and take ξ = (0.1, …, 0.1)^⊤, A = - + (a_ij)_i,j=1^d, B =  + (a_ij)_i,j=1^d, N =  + (a_ij)_i,j=1^d, where a_ij∼𝒩(0, 10^-4) is sampled i.i.d. for all i, j ∈{1, …, d} once at the beginning of the simulation, and C = diag(1, …, 1_r times, 0.01, …, 0.01_d - r  times) with r=10. The delay time is τ = 0.1, the control u(t) = sin(20 t) 1 and the history function f_t = 0. Figure <ref> shows the L^2 errors between the full and reduced systems and compares them with the bound from Theorem <ref>. We see that the error decreases if we choose r big enough and that the bound is off by around two orders of magnitude.
  
  
  
  
  
  
  
  
  
  
  
  
  Acknowledgements 
  This work was supported by the EPSRC grant EP/L016516/1 for the University of Cambridge CDT, the CCA (S.B.) and by Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1114 `Scaling Cascades in Complex Systems' (A05, project number 235221301). 
  
  
  
  
  
  
  
  
  
  
  
  0
  
  [BP05]BP Bátkai, A. and Piazzera, S. (2005). Semigroups for Delay equations. Research Notes in Mathematics. A K Peters/CRC Press.
  [BFS03]BFS Bátkai, A. and Fasanga, E. and Shvidkoy, R. (2003). Hyperbolicity of delay equations via Fourier multipliers, Acta scientiarum mathematicarum, (69) 1-2. pp. 131-145.
  [BD11]BD Benner, P. and Damm, T. (2011). Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic systems. SIAM Journal on Control and Optimization.
  [BG09]4 Beattie, C., Gugercin, S. (2009). Interpolatory projection methods for structure-preserving model reduction. Syst Control Lett 58(3):225-232.
  [BH19]BH18 Becker, S. and Hartmann, C. (2018). Infinite-dimensional bilinear-and stochastic balanced truncation.
  [BHRR20]BHRR Becker, S., Hartmann, C., Redmann, M., and Richter, L. (2020).  Feedback control theory and Model order reduction for stochastic equations, 1912.06113.
  [BR15]BR15 Benner, P. and Redmann, M. (2015). Model reduction for stochastic systems. Stoch PDE: Anal Comp,Volume 3, Issue 3, pp 291-338.
  [BS17]BS17 Butkovsky, O. and Scheutzow, M. (2017). Invariant measures for stochastic functional differential equations. Electron. J. Probab. Volume 22, paper no. 98, 23 pp.
  
  
  [CZ95]CZ  Curtain R.F. and Zwart H. (1995). An introduction to infinite-dimensional linear systems theory. Springer, New York.
  [CGP88]CGP Glover, K., Curtain, R., and Partington, J. (1988). Realisation and Approximation of Linear Infinite-Dimensional Systems with Error Bounds. SIAM Journal on Control and Optimization 26:4, 863-898.
  [EN00]EN Engel, K-J. and Nagel, R. (2000). One-Parameter Semigroups for Linear Evolution Equations. Springer. Graduate Texts in Mathematics.
  
  
  
  
  [GDBA19]GDBA Gosea, IV., Duff, IP., Benner, P., and Antoulas, AC. (2019). Model order reduction of bilinear time-delay systems, 18th European Control Conference (ECC), 2289-2294.
  [HD11]10 Harkort, C. and Deutscher, J. (2011). Krylov subspace methods for linear infinite-dimensional systems. IEEE. Trans Autom Control 56(2):441-447.
  [JDM13]JDM12 Jarlebring, E., Damm, T., and Wim, M. (2013). Model reduction of time-delay systems using position balancing and delay Lyapunov equations. Mathematics of Control, Signals, and Systems, Volume 25, Issue 2, pp 147-166.
  
  
  [K66]K Kubo, R. (1966). The fluctuation-dissipation theorem. Rep. Prog. Phys., 29(1):255.
  [LLL17]LLL Li, L., Liu, J.-G., and Lu, J. (2017). Fractional Stochastic Differential Equations Satisfying Fluctuation-Dissipation Theorem, J. Stat. Phys., 169, 2, 316-339.
  [M65]M Mori, M. (1965). A continued-fraction representation of the time-correlation functions.Prog. Theor. Phys.,34(3):399-416.
  
  [MS97]MS97 Mao, X. and Shah, A. (1997). Exponential stability of stochastic differential
  delay equations. Stochastics: An International Journal of Probability and Stochastic
  Processes, 60:1-2, 135-153
  [MP99]19 Mäkilä, P. and Partington, J. (1999). Laguerre and Kautz shift approximations of delay systems. Int J Control 72(10):932-946.
  [MP99b]20Mäkilä, P. and Partington, J. (1999). Shift operator induced approximations of delay systems. SIAM J Control Optim 37(6):1897-1912.
  
  [MJM11]24 Michiels, W., Jarlebring, E., and Meerbergen, K. (2011). Krylov-based model order reduction of time-delay
  systems. SIAM J Matrix Anal Appl 32(4):1399-1421.
  [P04]29 Partington, J. (2004) Model reduction of delay systems. In: Blondel V, Megretski A (eds) Unsolved problems in mathematical systems and control theory. Princeton university press, Princeton, pp 29-32.
  [PYPT10]PYPT Perlikowski, P., Yanchuk, S., Popovych, O., and Tass, P. (2010). Periodic patterns in a ring of delay-coupled oscillators. Phys. Rev. E 82, 036208.
  [R19]R1 Redmann, M. (2019).The missing link between the output and the H2-norm of bilinear systems, 1910.14427.
  [RS14]RS Reis, T. and Selig, T. (2014). Balancing Transformations for Infinite-Dimensional Systems with Nuclear Hankel Operator, Integr. Equ. Oper. Theory, Volume 79, Issue 1, pp 67-105.
  
  
  [SA16]SA Scarciotti, G. and Astolfi, A. (2016). Model Reduction of Neutral Linear and Nonlinear
  Time-Invariant Time-Delay Systems With
  Discrete and Distributed Delays. IEEE Transactions on automatic control, Vol. 61, No. 6.
  
  
  
  [WYH10]WYH Wolfrum, M., Yanchuk, S., Hövel, P. et al. (2010). Complex dynamics in delay-differential equations with large delay. Eur. Phys. J. Spec. Top. 191: 91.
  [ZL02]ZL Zhang, L. and Lam, J. (2002). On H^2 model reduction of bilinear systems. Automatica 38. 205-216. Pergamon.
  [ZZ12]ZZ Zhang, C. and Zheng, B. (2012). Z_n equivariant in delay coupled dissipative Stuart-Landau oscillators. Nonlinear Dyn 70:2359-2366.
  oscillators
  
  