\section{Introduction}
\label{MRSto-sec1}

\IEEEPARstart{D}{ynamical} systems described by stochastic differential equations have been successfully used in a variety of theoretical and applied scientific fields, such as system biology and finance \cite{KarShr:91,YonZho:99,Mun:11,Oks:13} (see also the seminal contributions of Kalman \cite{Kal:60,KalBuc:61} in control theory). One way to view these systems is to interpret the stochastic processes (\textit{i.e.} time sequences representing the evolution of variables which are subject to random variations) in the equations as a means to model \textit{uncertainty}. In this sense, stochastic systems offer a powerful modeling framework for engineering applications. However, as noted in \textit{e.g.} \cite{LinPic:15}, the complexity of the algorithms solving stochastic problems usually grows more than linearly with the dimension of the model. For this reason, many researchers have investigated the problem of model reduction for stochastic systems. The objective of model reduction is to obtain a \textit{simple}, in some sense to be defined, model of the original system which possesses, for some operating conditions, a specific subset of the properties of the system to be reduced. Most of this field originally developed from the need of reducing the order of high-order linear ordinary differential equations obtained from the discretisation of partial differential equations for simulation purposes, see \cite{Ant:05}. From there, the field expanded to much more general classes of systems, such as time-delay and nonlinear systems \cite{ScaAst:17}, and purposes, such as analysis and control, see for instance \cite{FaeScaAstRin:18,BreForScaAst:19,FaeScaAstRin:21}. Many stochastic generalizations of deterministic methods have been proposed. For instance for linear systems various forms of stochastic balancing truncation have been presented in \textit{e.g.}\cite{DesPal:84,DesPalKir:85,Tug:85,Gre:88,Gre:88a,HarJonSil:84,WanSaf:90,WanSaf:91,LinPic:96}, whereas stochastic $\mathcal{H}_{\infty}$ model reduction has been given in \textit{e.g.} \cite{XuChe:03,SuWuShiSon:12}. For bilinear systems, balanced truncation for stochastic systems has been proposed in \cite{BenDam:11,BenRed:15}. For nonlinear systems, the slaving principle, the slow-fast reduction based on center manifold theory,  and the use of eigenfunctions are all methods which have been extended to stochastic systems, see \cite{SchHak:86,SchHak:87}, \cite{XuRob:96,Rob:08,Rob:14} and \cite{CoiKevLafMagNad:08}. Moreover, model reduction of stochastic systems continues to be an active area of research for applications such as quantum stochastic systems, see \textit{e.g.} \cite{Nur:14,TecNur:16}, and biological and chemical systems, see \textit{e.g.} \cite{PelMunKha:06,SinHes:10,BruChaSmi:14,SooAnd:14,GupKha:14,MelHesKha:14,JohBarSejMjo:15,SmiCiaGri:15}.
However, other successful deterministic approaches for model reduction of nonlinear systems, such as \cite{Sch:93} (nonlinear balancing) and \cite{Ast:10} (nonlinear moment matching), have not yet been extended to the stochastic framework.

In this paper we address the problem of model reduction for a general class of linear and nonlinear stochastic systems using a moment matching approach. By ``general class'', we mean systems in which the state and the input may appear simultaneously in both the drift term and diffusion term of the equation, \textit{i.e.} they may multiply both the $dt$ term and the $d\W_t$ term, where $\W_t$ is a Brownian motion. Borrowing from the deterministic literature, we formulate the problem of model reduction by moment matching as the problem of determining a reduced-order model that possesses the same steady-state output of the system to be reduced for specific classes of inputs. For the sake of studying the problem in its most general form, also the signal generator which produces the input signals of interest is selected as a stochastic system. To the end of matching the steady-state output of the system, we construct a stochastic process, which we call moment, that describes the steady-state behaviour of the system. The moment is the solution of a stochastic generalization of the Sylvester equation for linear systems and of the invariance partial differential equation for nonlinear systems. We then propose families of reduced-order models that match the moment of the system to be reduced. Unfortunately, it is shown that the reduced-order models have an on-line computational complexity which is comparable to the one of the original system. To overcome this issue two strategies are adopted to provide families of approximated models. The first family preserves a part of the stochastic properties of the moment, whereas the second family preserves a part of the stochastic properties of the steady-state output. Finally, the differences between these families of reduced-order models are illustrated by means of simulations of numerical examples and of a stochastic variation of a deterministic benchmark system given in \cite{ChaVaD:05,SLICOT}.

Preliminary versions of our work have been published in \cite{ScaTee:17,ScaTee:17a}. Note that considering a stochastic signal generator is a nontrivial generalisation of \cite{ScaTee:17,ScaTee:17a}. First, many of the assumptions have to be revised and relaxed. Second, a stochastic signal generator has the peculiarity of changing the mean of the mappings describing the steady-state response, which is a feature that is not captured in \cite{ScaTee:17,ScaTee:17a}. In this regard, a new family of models (moment-mean) is proposed and several variations of the results are sketched (for instance considering multiple Brownian motions). Third, the complete proofs of the results are provided. These are of independent interest beyond the model reduction literature because many of the results are instrumental to the solution of other mathematical control problems, such as output regulation \cite{Isi:95,Sca:18,MelSca:19,MelSca:19a,MelSca:20}. Finally, new simulations, numerical and based on a benchmark example, are provided.

The rest of the paper is organized as follows. In Section~\ref{MRSS-sec-PN} we formulate the problem. In Section~\ref{MRSto-SecSS} we characterize the moment for linear and nonlinear stochastic systems. In Section~\ref{MRNS-secROM} we provide families of reduced-order models which match the moment of the system to be reduced. In Section~\ref{MRSS-sec-AROM} we first point out that the found models have a high computational complexity and then we propose new classes of approximated reduced-order models. In Section~\ref{MRSS-sec-Sim} we illustrate some of the results by means of simulations. Section~\ref{MRSS-sec-Con} contains our concluding remarks. In order to make the paper as self-contained as possible, the most important concepts of stochastic systems used throughout the paper are reported, with plenty of references, in the Appendix.

\textbf{Notation.} We use standard notation. $\mathbb{C}_{<0}$ ($\mathbb{C}_{\ge 0}$) denotes the set of complex numbers with negative (non-negative) real part. $\mathbb{R}_{<0}$ ($\mathbb{R}_{> 0}$) denotes the set of negative (positive) real numbers. The symbol $I$ denotes the identity matrix and $\sigma(A)$ denotes the spectrum of the matrix $A\in\mathbb{R}^{n\times n}$. The symbols $|v|$, with $v\in\R^n$, and $||A||$ indicates the Euclidean norm and the induced Euclidean matrix norm, respectively. The vectorization of a matrix $A\in\mathbb{R}^{n\times m}$, denoted by $\vect{A}$, is the $nm \times 1$ vector obtained by stacking the columns of the matrix $A$ one on top of the other, namely $\vect{A}=[a_1^{\top},a_2^{\top},\dots,a_m^{\top}]^{\top}$, where $a_i\in\mathbb{R}^n$ is the $i$-th column of $A$ and the superscript $\top$ denotes the transposition operator. The symbol $\otimes$ indicates the Kronecker product, whereas the symbol $\oplus$ indicates the direct sum. Given two functions, $f:Y\to Z$ and $g:X\to Y$, with $f\,\circ\,g:X\to Z$ we denote the composite function $(f\,\circ\,g)(x)=f(g(x))$ which maps all $x\in X$ to $f(g(x))\in Z$. $(\Omega,\mathcal{F},\mathcal{P})$ indicates a probability space with a given set $\Omega$, a $\sigma$-algebra $\mathcal{F}$ on $\Omega$ and a probability measure $\mathcal{P}$ on the measurable space $(\Omega,\mathcal{F})$. For ease of notation, we often indicate a stochastic process $\{x_t,\,t\in\mathbb{R}\}$ simply with $x_t$ (this is common in the literature, see \textit{e.g.} \cite{Arn:74}). 
The stochastic process $\W_t$ indicates a standard Wiener process defined on the probability space $(\Omega,\mathcal{F},\mathcal{P})$. $\mathcal{F}_t$ is the continuous-time filtration generated by the Wiener process $\W_t$ up to time $t$ and all stochastic processes appearing in this paper are adapted to this filtration (and, possibly, to others when multiple Brownian motions appear). Let $X_t$ and $Y_t$ be any two stochastic processes and $d: (t,x,y) \mapsto d_t(x,y) $ any function which has well-defined partial derivatives; then any partial derivative of the form 
$\frac{\partial d_t}{\partial x} (X_t,Y_t)$ is compactly indicated as $\frac{\partial d_t}{\partial X_t}$.
All the stochastic integrals in this paper are intended as It\^{o} integrals.

\section{Problem formulation}
\label{MRSS-sec-PN}

In this section we formulate the problem of model reduction by moment matching for stochastic systems. Complete definitions of stochastic process, Brownian motion, It\^{o}'s formula, stability and so on are reported in the Appendix. The reader who is not familiar with stochastic differential equations is invited to consult the Appendix and references therein.

Consider a stochastic nonlinear single-input single-output continuous-time system $\Sys_{n}(f,g,h)$ described by the equations
\begin{equation}
\label{MRNS-eq1}
dx_t=f(x_t,u_t)dt+g(x_t,u_t)d\W_t, \qquad y_t=h(x_t),
\end{equation}
with $x_t\in\mathbb{R}^n$, $u_t\in \mathbb{R}$, $y_t\in\mathbb{R}$, and $f$, $g$ and $h$ smooth mappings.
Consider a signal generator $\Sys_{\nu}(s,j,l)$ described by the equations
\begin{equation}
\label{MRNS-eq2}
d \omega_t=s(\omega_t)dt + j(\omega_t)d\W_t , \qquad u_t=l(\omega_t),
\end{equation}
with $\omega_t\in \mathbb{R}^{\nu}$, and $s$, $j$ and $l$ smooth mappings.
Consider the interconnection of system~(\ref{MRNS-eq1}) with the signal generator~(\ref{MRNS-eq2}), namely
\begin{equation}
\label{MRNS-eqInt}
\!\!\!\!\!\begin{array}{l}
\left[\!\!\!\begin{array}{c}d \omega_t \\ dx_t \end{array}\!\!\!\right]\!=\!\left[\!\!\!\begin{array}{c} s(\omega_t) \\ f(x_t,l(\omega_t)) \end{array}\!\!\!\right] \!dt \!+\! \left[\!\!\!\begin{array}{c}j(\omega_t) \\ g(x_t,l(\omega_t))\end{array}\!\!\!\right]\!d\W_t, \quad y_t=h(x_t).
\end{array}
\end{equation}
Assume that zero is an equilibrium point of (\ref{MRNS-eqInt}), \textit{i.e.} $s(0)=0$, $j(0)=0$, $l(0)=0$, $f(0,0)=0$, $g(0,0)=0$ and $h(0)=0$, and that the initial condition $(\omega(0),x(0))$ is deterministic.






We now define the ``moment matching condition'' and introduce the ``problem of model reduction by moment matching''. An intuitive explanation follows the problem.

\begin{definition}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). A stochastic system $\Sys_{\nu}(\red{f},\red{g},\red{h})$ described by the equations
\begin{equation}
\label{MRSto-eqRM}
d\red{x}_t=\red{f}(\red{x}_t,u_t)dt+\red{g}(\red{x}_t,u_t)d\W_t, \qquad \red{y}_t=\red{h}(\red{x}_t),
\end{equation}
where $\red{x}_t\in\mathbb{R}^{\nu}$, with $\nu < n$, $\red{y}_t\in\mathbb{R}$, and $\red{f}$, $\red{g}$ and $\red{h}$ are smooth mappings, is said to satisfy the \textit{moment matching condition at $(s,j,l)$} if the error $e_t=y_t-\red{y}_t$, where $y_t$ is the output of (\ref{MRNS-eqInt}) and $\red{y}_t$ is the output of (\ref{MRSto-eqRM}) driven by (\ref{MRNS-eq2}), satisfies
\begin{equation}
\label{MRSto-eqMMc}
\lim_{t\to \infty} e_t = 0
\end{equation}
almost surely, for any $(x_0,\omega_0,\red{x}_0)\in \mathcal{N} \subset \mathbb{R}^{n} \times \mathbb{R}^{\nu} \times \mathbb{R}^{\nu}$.
\end{definition}



\begin{problem}
\label{MRSto-pro}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). The \textit{problem of model reduction by moment matching} consists in determining a stochastic system $\Sys_{\nu}(\red{f},\red{g},\red{h})$, with $\nu<n$, which satisfies the moment matching condition at $(s,j,l)$.
\end{problem}

\begin{figure} {\color{black}
\centering 
\tikzset{%
  block/.style    = {draw, thick, rectangle, minimum height = 3em,
    minimum width = 3em,text width=5.5em, text centered},
  sum/.style      = {draw, circle, node distance = 2cm},
  dot/.style      = {coordinate},
  input/.style    = {coordinate},
  output/.style   = {coordinate}
}
\begin{tikzpicture}[auto, thick, node distance=2cm, >=triangle 45]
\draw

	node at (0,-1)[block] (input1) {$\Sys_{\nu}(s,j,l)$} 
	node at (1.7,-1)[dot] (suma1) {}
	node at (3.5,0)[block] (inte1) {$\Sys_{n}(f,g,h)$}
    node at (5.3,-1)[sum] (suma2) {$-$}
    node at (3.5,-2.0)[block] (ret1) {$\Sys_{\nu}(\red{f},\red{g},\red{h})$}

	node at (7.4,-1) [output] (sal2){}; 
   

    \draw[-](input1) -- node {$u_t$} (suma1);
    \draw[->](suma1) |- node {} (inte1);
   
 	\draw[->](suma1) |- node {} (inte1);
	\draw[->](inte1) -| node [near end, left] {$y_t$} (suma2);
	\draw[->](ret1) -| node [near end, left] {$\red{y}_t$} (suma2);
	\draw[->](suma1) |- node {} (ret1);

	\draw[->] (suma2) -- node {$ e_t \to 0$ a.s.}(sal2);
\end{tikzpicture} }
\caption{A schematic overview of Problem~\ref{MRSto-pro}: given $\Sys_{\nu}(s,j,l)$ and $\Sys_{n}(f,g,h)$ we want to determine $\Sys_{\nu}(\red{f},\red{g},\red{h})$ such that $ e_t \to 0$ almost surely.} 
\label{MRDI-fig1}
\end{figure}

The interpretation of Problem~\ref{MRSto-pro}, as illustrated in Fig.~\ref{MRDI-fig1}, is that we are looking for a family of \textit{simpler} models (where simplicity is somewhat arbitrarily intended as $\nu<n$) which behave asymptotically as the system that we want to reduce when both the system and the reduced-order model are driven by an \textit{a priori} selected signal generator. In other words, we are interested in preserving the asymptotic behavior of the system for specific operating conditions and input signals. This ``steady-state matching'' method is called ``moment matching'' because \cite{Ast:10} recognized that for deterministic linear systems the method is equivalent to matching the ``moments'' (\textit{i.e.} the coefficients of a series expansion of the transfer function) as defined in \cite{Ant:05}. We anticipate from the onset that we also solve ``approximated'' versions of Problem~\ref{MRSto-pro} in which the moment matching condition~(\ref{MRSto-eqMMc}) is relaxed in various ways.


At this point, it is useful to introduce the notation for stochastic linear systems because, in addition to providing stronger and more easily interpretable results, this class of systems has a large practical use in the deterministic model reduction literature. Hence, we expect that the linear results of this paper will have a more immediate impact. When system~(\ref{MRNS-eq1}) is linear, we use the notation
\begin{equation}
\label{ORS-eq1}
dx_t=(Ax_t+Bu_t)dt+(Fx_t+Gu_t)d\W_t,\qquad y_t=Cx_t,
\end{equation}
with $A\in \mathbb{R}^{n\times n}$, $B\in \mathbb{R}^{n\times1}$, $F\in \mathbb{R}^{n\times n}$, $G\in \mathbb{R}^{n\times1}$ and $C\in \mathbb{R}^{1\times n}$. Similarly, when the generator~(\ref{MRNS-eq2}) is linear, we use the notation 
\begin{equation}
\label{ORS-eq2}
d\omega_t= S\omega_t dt+J\omega_t d\W_t, \qquad u_t=L\omega_t,
\end{equation}
with $S\in \mathbb{R}^{\nu\times \nu}$, $J\in \mathbb{R}^{\nu\times\nu}$ and $L\in \mathbb{R}^{1\times \nu}$.
Let $\Phi_t\in\mathbb{R}^{n\times n}$ be the fundamental matrix of the homogeneous equation corresponding to (\ref{ORS-eq1}), \textit{i.e.}
\begin{equation}
\label{ORS-eqPhi}
d\Phi_t =  \left( A  dt + F d\W_t\right) \Phi_t,
\end{equation}
with $\Phi_0=I$ and recall that (see \cite[Section 4.1]{Gar:88}) 
\begin{equation}
\label{ORS-eqPhiInv}
d\Phi_t^{-1} = \Phi_t^{-1} \left( (F^2-A) dt - F d\W_t \right).
\end{equation}
Finally, let $\Sigma_t\in\mathbb{R}^{\nu\times \nu}$ be the fundamental matrix corresponding to (\ref{ORS-eq2}), \textit{i.e.}
\begin{equation}
\label{ORS-eqSigma}
d\Sigma_t =  \left( S  dt + J  d\W_t \right) \Sigma_t ,
\end{equation}
holds.
























\section{Steady state of stochastic systems}
\label{MRSto-SecSS}

Since the moment matching condition~(\ref{MRSto-eqMMc}) is a condition on the steady-state behavior of the system, it is instrumental for the solution of the problem to provide a description of the steady-state response of the system. In this section we characterize the steady state of system~(\ref{ORS-eq1}) driven by (\ref{ORS-eq2}) in terms of a stochastic partial differential equation, and of system~(\ref{ORS-eq1}) driven by (\ref{ORS-eq2}) in terms of a stochastic Sylvester equation. Note that these two results, which are the main technical contributions of the paper, have an interest beyond the problem of model reduction (because, for instance, similar characterizations of the steady-state response are an essential element of the solution of output regulation problems \cite{Isi:95,Sca:18}).

We begin with the linear case, which is easier to develop and is instrumental for the nonlinear case. At the end of the section we  show that the results can be formulated also for the simpler case of  multiple uncorrelated Brownian motions acting on the system and/or the signal generator.

\subsection{Steady state: linear stochastic systems}
\label{MRSto-SubLSS}

We introduce the following assumptions\footnote{For the definition of Lyapunov exponent and a procedure to check these assumptions, we refer the reader to the Appendix.}.
\begin{assumption}
\label{ORS-asSys}
All Lyapunov exponents of $\Phi_t$ are negative almost surely.
\end{assumption}
\begin{assumption}
\label{ORS-asGen}
All Lyapunov exponents of $\Sigma_t$ are zero almost surely.
\end{assumption}

We are now ready to give a characterization of the steady-state response of system~(\ref{ORS-eq1}) driven by~(\ref{ORS-eq2}).




\begin{theorem}
\label{ORS-thLSSX}
Consider the interconnection of system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). Suppose that Assumptions~\ref{ORS-asSys} and \ref{ORS-asGen} hold. Then the steady-state response of the output of such interconnection is
$$
y^{ss}_t=C \X_t\omega_t
$$
almost surely, where $\X_t\in\mathbb{R}^{n\times \nu}$ is
\begin{equation}
\label{ORS-eq4}
\!\!\X_t\!=\!\Phi_t \!\left[\int_{-\infty}^t\!\!\!\!\!\! \Phi_{\tau}^{-1}(BL -FGL)\Sigma_{\tau}d\tau\!+\!\!\!\int_{-\infty}^t \!\!\!\!\!\! \Phi_{\tau}^{-1}GL\Sigma_{\tau}d\W_{\tau}\right]\!\Sigma_{t}^{-1}\!\!.
\end{equation}
The stochastic process $\X_t$ is the steady-state solution of the stochastic differential matrix equation
\begin{equation}
\label{ORS-eqPI}
\begin{array}{rl}
d\X_t=\!\!\!\!&\left(A\X_t-\X_t\left(S-J^2\right)-F\X_t J+BL-GLJ\right)dt\\[2mm]
&+\left(F\X_t-\X_t J+GL\right)d\W_t.
\end{array}
\end{equation}
\end{theorem}

\textit{Proof:}
Consider the matrix $\X_t$ defined in (\ref{ORS-eq4}). Multiplying this equation by $\Phi_t^{-1}$ on the left and by $\Sigma_t$ on the right yields
\begin{equation}
\label{MRLS-eqSSp1}
\begin{array}{l}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\int_{-\infty}^t  \!\!\!\!\!\! \Phi_{\tau}^{-1}(BL-FGL) \Sigma_{\tau} d\tau + \int_{-\infty}^t  \!\!\!\!\!\! \Phi_{\tau}^{-1}GL \Sigma_{\tau}d\W_{\tau},\\
=\Phi_t^{-1}\X_t \Sigma_t=\int_{-\infty}^t d\!\left(\Phi_{\tau}^{-1} \X_{\tau} \Sigma_{\tau} \right),
\end{array}
\end{equation}
in which the first equality holds because of (\ref{ORS-eq4}) and the second equality holds because $\textstyle\lim_{t\to-\infty}\Phi_t^{-1}\X_t \Sigma_t=0$ almost surely. In fact, by Assumptions~\ref{ORS-asSys} and \ref{ORS-asGen}, $\textstyle\lim_{t\to-\infty}\Phi_t^{-1}=0$ exponentially, $\textstyle\lim_{t\to-\infty}\frac{1}{t}\Sigma_t=0$ and $\textstyle\lim_{t\to-\infty}\frac{1}{t}\X_t=0$. This last limit holds for the same reasons noting that (\ref{ORS-eq4}) multiplied by $\omega_t$ is a response of (\ref{ORS-eq1}) to the input (\ref{ORS-eq2}). We can now use Leibniz integral rule to differentiate the integrals in both sides of~(\ref{MRLS-eqSSp1}) obtaining
$$
\begin{array}{l}
\Phi_t^{-1}(BL-FGL)\Sigma_tdt+\Phi_t^{-1}GL\Sigma_td\W_t=d(\Phi_t^{-1}\X_t\Sigma_t)\\[2mm]
\,\,\,\,\,=d\Phi_t^{-1}\X_t \Sigma_t+\Phi_t^{-1}d\X_t \Sigma_t+\Phi_t^{-1}\X_t d\Sigma_t\\[2mm]
\,\,\,\,\,+d\Phi_t^{-1}d\X_t\Sigma_t+\Phi_t^{-1}d\X_t d\Sigma_t+d\Phi_t^{-1}\X_t d\Sigma_t,
\end{array}
$$
where the term $d\Phi_t^{-1}d\X_t d\Sigma_t$ does not appear because it is zero by~(\ref{ORS-eqIto1}) in the Appendix. Substituting (\ref{ORS-eqPhiInv}) and (\ref{ORS-eqSigma}) in the previous equation and multiplying on the left by $\Phi_t$ and on the right by $\Sigma_t^{-1}$, yields
$$
\begin{array}{rl}
\left((F^2\right.\!\!\!\!\!\!&-\left.A)dt-Fd\W_t\right)\X_t+d\X_t+\X_t (Sdt+Jd\W_t) \\[2mm]
&+\left((F^2-A)dt-Fd\W_t\right)d\X_t+d\X_t (Sdt+Jd\W_t)\\[2mm]
&+\left((F^2-A)dt-Fd\W_t\right)\X_t(Sdt+Jd\W_t)\\[2mm]
=\!\!\!\!&(BL-FGL)dt+GLd\W_t.
\end{array}
$$
By (\ref{ORS-eqIto1}) we note immediately that in the second line $dtd\X_t=0$ and in the third line only the term in $d\W_t d\W_t=dt$ is not zero. By sorting this expression we obtain
\begin{equation}
\label{MRLS-eqSSp2}
\begin{array}{rl}
\!\!d\X_t&\!\!\!=\left((A-F^2)\X_t-\X_t S + F\X_t J +BL-FGL\right)dt\\[2mm]
&\!\!\!+(F\X_t-\X_t J+GL)d\W_t + (Fd\X_t -d\X_t J)d\W_t.
\end{array}
\end{equation}
Multiplying (\ref{MRLS-eqSSp2}) first by $Fd\W_t$ on the left and then (separately) by $Jd\W_t$ on the right and using (\ref{ORS-eqIto1}) and (\ref{ORS-eqIto2}) in the Appendix, yields
$$
\begin{array}{rcl}
Fd\X_t d\W_t \!\!\!&=&\!\!\!(F^2\X_t-F\X_t J +FGL)dt,\\
-d\X_t J d\W_t \!\!\!&=&\!\!\! (-F\X_t J+\X_t J^2 -GLJ)dt.
\end{array}
$$
Substituting these two expressions in (\ref{MRLS-eqSSp2}) proves that $\X_t$, defined in (\ref{ORS-eq4}), is the solution of the stochastic differential matrix equation~(\ref{ORS-eqPI}).
Now define the variable $z_t:=x_t-\X_t\omega_t$. Then by the stochastic product rule (recalled as Lemma~\ref{MRSS-lmdd} in the Appendix)
\begin{equation}
\label{MRSS-eqdz}
dz_t=dx_t-d\X_t\omega_t-\X_t d\omega_t - d\X_t d\omega_t.
\end{equation}
Note that $d\X_t d\omega_t = \left(F\X_t-\X_t J+GL\right)J\omega_t dt$
by (\ref{ORS-eqIto1}), (\ref{ORS-eqIto2}) and (\ref{ORS-eqPI}). Substituting the expressions of $dx_t$, namely (\ref{ORS-eq1}), of $d\X_t$, namely (\ref{ORS-eqPI}), and of $d\omega_t$, namely (\ref{ORS-eq2}), in (\ref{MRSS-eqdz}) yields
$$
\!dz_t\!=\!A(x_t-\X_t\omega_t)dt+F(x_t-\X_t\omega_t)d\W_t\!=\! Az_tdt+Fz_td\W_t\!.
$$
The zero equilibrium $z_t=0$ of this last equation is asymptotically stable almost surely because of Assumption~\ref{ORS-asSys}. In turn this proves that $x_t$ converges to $\X_t\omega_t$ almost surely as $t\to+\infty$. The claim follows substituting the steady-state response $x_{ss}=\X_t \omega_t$ in the equation of $y_t$ in (\ref{ORS-eq1}).
\hspace*{0pt}\hfill $\square$ \par







\subsection{Steady state: nonlinear stochastic systems}
\label{MRSto-SubNLSS}

We now generalize the previous result to nonlinear systems, \textit{i.e.} we show that the steady-state response of the interconnection of system~(\ref{MRNS-eq1}) and (\ref{MRNS-eq2}), namely (\ref{MRNS-eqInt}), can be characterized by a stochastic partial differential equation which is a generalization of \cite[Equation~(8.3)]{Isi:95}. To streamline the presentation we formulate the following assumptions.
\begin{assumption}
\label{MRNS-as1}
Assumption~\ref{ORS-asSys} holds for the linearization of system~(\ref{MRNS-eq1}) around the zero equilibrium.
\end{assumption}
\begin{assumption}
\label{MRNS-as2}
Assumption~\ref{ORS-asGen} holds for the linearization of system~(\ref{MRNS-eq2}) around the zero equilibrium.
\end{assumption}

These two assumptions induce a decomposition of the state space $\mathbb{R}^{n+\nu}$ as $\mathbb{R}^{n+\nu}=E_s(w)\oplus E_c(w)$, where $E_s$ and $E_c$ are  the stable and central Odeselec spaces, in a similar way as negative and zero real part eigenvalues induce a decomposition in stable and central subspaces for deterministic systems (see the Appendix for more details).
\begin{theorem}
\label{MRNS-thm1}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). Suppose Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} hold. Then there exists a stochastic process $\chi:(t,\W_t,\omega_t)\mapsto \chi_t(\W_t,\omega_t)$, locally defined in a neighborhood $W(w) \subset E_c(w)$ of $\omega_t=0$, with $\chi_t(0,0)=0$, which solves the stochastic partial differential equation
\begin{equation}
\label{MRNS-eqNSS}
d \chi_t= f(\chi_t(\W_t, \omega_t),l(\omega_t))dt + g(\chi_t(\W_t, \omega_t),l(\omega_t))d\W_t,
\end{equation}
where
\begin{equation}
\label{MRNS-eqSPDE}
\begin{array}{rl}
\!\!\!d \chi_t &\!\!\!\!=\left[\frac{\partial \chi_t}{\partial t}+ \frac{\partial \chi_t}{\partial \omega_t}s(\omega_t) + \frac{1}{2} \frac{\partial^2 \chi_t}{\partial \W_t^2} \right. + \frac{\partial^2 \chi_t}{\partial \W_t \partial \omega_t}j(\omega_t)\\
&\!\!\!\!\left.+\frac{1}{2} j(\omega_t)^{\top}\frac{\partial^2 \chi_t}{\partial \omega_t^2} j(\omega_t) \right]dt + \left[\frac{\partial \chi_t}{\partial \W_t} + \frac{\partial \chi_t}{\partial \omega_t}j(\omega_t) \right]d\W_t,
\end{array}
\end{equation}
for all $\omega_t \in W(w)$. In addition, the steady-state response of system~(\ref{MRNS-eqInt}) is $x^{ss}_t=\chi_t(\W_t,\omega_t)$ almost surely for any $(\omega_0,x_0)\in W(w)\times X(w)$, where $X(w)\subset E_s(w)$ is a neighborhood of the origin.
\end{theorem}

\textit{Proof:}
Consider the interconnection~(\ref{MRNS-eqInt}) of system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). Computing the linearization at the zero equilibrium of this system yields
$$
\left[\!\! \begin{array}{c}d\bar{\omega}_t \\ d\bar{x}_t \end{array}\!\!\right]= \left[ \begin{array}{cc}S & 0\\ \tilde{B} & A \end{array}\right]\left[ \begin{array}{c}\bar{\omega}_t \\ \bar{x}_t\end{array}\right]dt +\left[ \begin{array}{cc}J & 0\\\tilde{G} & F \end{array}\right]\left[ \begin{array}{c}\bar{\omega}_t \\ \bar{x}_t \end{array}\right]d\W_t,
$$
where 
$$
\left.A=\frac{\partial f(x_t,l(\omega_t))}{\partial x_t}\right|_{\substack{x_t=0\\\omega_t=0}}, \quad \left.\tilde{B}=\frac{\partial f(x_t,l(\omega_t))}{\partial \omega_t}\right|_{\substack{x_t=0\\\omega_t=0}},
$$
$$
\left.F=\frac{\partial g(x_t,l(\omega_t))}{\partial x_t}\right|_{\substack{x_t=0\\\omega_t=0}}, \quad \left.\tilde{G}=\frac{\partial g(x_t,l(\omega_t))}{\partial \omega_t}\right|_{\substack{x_t=0\\\omega_t=0}},
$$ 
$$
\left.S=\frac{\partial s(\omega_t)}{\partial \omega_t}\right|_{\omega_t=0}, \quad \left.J=\frac{\partial j(\omega_t)}{\partial \omega_t}\right|_{\omega_t=0}.
$$
By Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2}, Theorem~\ref{ORS-thLSSX} applies to this linearized system. Thus, by appealing to the center manifold theory (see \cite[Section 2.3]{Car:81} for deterministic systems, and \cite{Box:89} for a stochastic version), there exists a neighborhood $W(w)\times X(w)$ of $(0,0)$ in which the interconnected system has a center manifold at $(\omega_t,x_t)=(0,0)$ described by the graph of $x_t=\chi_t(\W_t,\omega_t)$. By replacing $\chi_t$ in (\ref{MRNS-eq1}) we obtain the characterization given by equation~(\ref{MRNS-eqNSS}), with the property that $\textstyle\left.\X_t=\frac{\partial \chi_t}{\partial \omega_t}\right|_{\omega_t=0}$, where $\X_t$ is the solution of (\ref{ORS-eqPI}). In addition, by \cite[Theorem 7.1(i)]{Box:89} the stochastic center manifold is locally exponentially attractive almost surely \textit{i.e.} for all pairs $(\omega_0^\star,x_0^{\circ})$ in the neighborhood $W(w)\times X(w)$ of $(0,0)$, the inequality\footnote{The definition of the norm $|\cdot|_w$ is given in the Appendix (Definition~\ref{def-rn}).}
$$
| x_t-\chi_t(\W_t,\omega_t)|_w \le K(t) | x_0^{\circ}-\chi_0(0,\omega_0^\star)|_w,
$$
with $K$ such that $\textstyle \lim_{t \to +\infty}\frac{1}{t}\log K(t)<0$, holds for all $t\ge 0$ almost surely. This together with the invariance of the stochastic center manifold proves that the steady state of system~(\ref{MRNS-eq1}) driven by~(\ref{MRNS-eq2}) is described almost surely by $x_t^{ss}=\chi_t(\W_t,\omega_t)$ for any $(\omega_0,x_0)\in W(w)\times X(w)$. Finally, equation~(\ref{MRNS-eqSPDE}) for the differential $d\chi_t$ is a direct consequence of It\^{o}'s formula (recalled as Lemma~\ref{MRSS-lmIto} in the Appendix).
\hspace*{0pt}\hfill $\square$ \par

\begin{remark}
If $g\equiv 0$ and $j\equiv 0$, then $\chi_t$ is a deterministic mapping. In this case, all the partial derivatives in (\ref{MRNS-eqSPDE}) are zero apart for $\frac{\partial \chi_t}{\partial \omega_t}s(\omega_t)$. Thus, by setting $\pi=\chi_t$, where $\pi:\omega\mapsto \pi(\omega)$, equation~(\ref{MRNS-eqNSS}) reduces to the well-known partial differential equation \cite[Equation~(8.3)]{Isi:95}
\begin{equation}
\label{MRSSeqISI}
\frac{\partial \pi}{\partial \omega}s(\omega) = f(\pi(\omega),l(\omega)).
\end{equation}
\end{remark}

\begin{remark}
Similarly to Theorem~\ref{ORS-thLSSX}, Theorem~\ref{MRNS-thm1} establishes that among the solutions of~(\ref{MRNS-eqNSS}), there exists at least one attractive solution that can be used to describe the steady state of the interconnection of system~(\ref{MRNS-eq1}) with the signal generator~(\ref{MRNS-eq2}).
\end{remark}

Inspired\footnote{The reason for using the word ``moment'' in the proposed stochastic generalization is that if we restrict the system to be linear and deterministic, \textit{i.e} if $f(x,u)=Ax+Bu$, $g(x,u) = 0$ and $h(x) =Cx$, then the generalized definition collapses into the classical definition of moment.} by \cite[Chapter 2.2]{ScaAst:17}, we can now define the moment in the stochastic framework.

\begin{definition}
\label{MRSS-DefMom}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). The \textit{moment of system~(\ref{MRNS-eq1}) at $(s,j,l)$} is the mapping $h\circ \chi_t$, where $\chi_t$ is given in Theorem~\ref{MRNS-thm1}.
\end{definition}
\begin{remark}
The mapping $h\circ \chi_t$ is a function of time and of the stochastic process $\W_t$. Moments as functions of time have been introduced in \cite{ScaAst:15c} for discontinuous signal generators, in \cite{ScaAst:16d} for some classes of hybrid systems and in \cite{ScaTeeAst:17} for linear differential inclusions. Hence, the moment of system~(\ref{MRNS-eq1}), which sometimes we call ``stochastic moment'', generalizes both the classical time-invariant moments (see \textit{e.g.} \cite{Ast:10,ScaAst:15b}) and the time-varying moments introduced in those papers. 
\end{remark}

The differential defined in (\ref{MRNS-eqSPDE}) looks intimidating. For the sake of providing a worked example on how to manipulate (\ref{MRNS-eqSPDE}), we now show how to obtain~(\ref{ORS-eqPI}) directly from (\ref{MRNS-eqSPDE}) when the system and the generator are linear. 

\begin{example}
Consider system~(\ref{ORS-eq1}) and the generator~(\ref{ORS-eq2}) and assume $\chi_t=\X_t \omega_t$. In this example we show that equation~(\ref{MRNS-eqNSS}) implies equation~(\ref{ORS-eqPI}). First of all, equation~(\ref{MRNS-eqSPDE}) becomes
$$
\begin{array}{rl}
\!\!d\chi_t&\!\!\!\!= \left[\frac{\partial \X_t}{\partial t}\omega_t + \X_t S \omega_t + \frac{1}{2}\frac{\partial^2 \X_t}{\partial \W_t^2}\omega_t + \frac{\partial \X_t}{\partial \W_t}J\omega_t + 0\right]dt\\
&\!\!\!\!+\left[\frac{\partial \X_t}{\partial \W_t}\omega_t +\X_t J\omega_t \right]d\W_t.
\end{array}
$$
By It\^{o}'s Lemma the stochastic process $\X_t$ has the differential
$$
d\X_t = \left[\frac{\partial \X_t}{\partial t} + \frac{1}{2}\frac{\partial^2 \X_t}{\partial \W_t^2}\right]dt + \frac{\partial \X_t}{\partial \W_t}d\W_t.
$$
Comparing the previous two equations yields
$$
d\chi_t=d\X_t\omega_t + \left[\X_t S \omega_t + \frac{\partial \X_t}{\partial \W_t}J\omega_t\right]dt + \X_t J\omega_t d\W_t.
$$
Thus, equation (\ref{MRNS-eqNSS}) becomes
$$
\begin{array}{rl}
d\X_t\omega_t  \!\!\!\!&=\left(A\X_t-\X_t S+BL-\frac{\partial \X_t}{\partial \W_t}J\right)\omega_tdt\\[2mm]
&+\left(F\X_t-\X_t J+GL\right)\omega_td\W_t.
\end{array}
$$
By factoring out $\omega_t$ and noticing that the previous equation implies $\textstyle\frac{\partial \X_t}{\partial \W_t}=F\X_t-\X_t J+GL$, yields~(\ref{ORS-eqPI}).
\end{example}

Note that there are simpler ways to achieve the same result, as pointed out in the next remark.

\begin{remark}
The previous result can be obtained directly from $d(\X_t\omega_t)$ using the stochastic product rule. In fact, by Lemma~\ref{MRSS-lmdd} we have
$$
\begin{array}{rl}
\!\!\!d(\X_t\omega_t)&\!\!\!\!\!=\!d\X_t\omega_t+\X_td\omega_t+d\X_td\omega_t\\
&\!\!\!\!\!=\!d\X_t\omega_t+\X_t(S\omega_t dt + J \omega_t d\W_t) + [\X_t]_{\W} J \omega_t dt,
\end{array}
$$
where the notation $[\X_t]_{\W}$ indicates the component of $d\X_t$ which multiplies $d\W_t$. Hence, (\ref{MRNS-eqNSS}) becomes
$$
\begin{array}{l}
\!\!\!d\X_t+ [\X_t]_{\W} J dt \\[2mm]
\qquad =(A\X_t-\X_t S + BL)dt+(F\X_t -\X_t J + GL)d\W_t
\end{array}
$$
from which it follows that $[\X_t]_{\W}=F\X_t-\X_t J+GL$, thus proving the claim. 
\end{remark}

This last remark provides us with a quick analytical approach to reformulate the results of Theorem~\ref{ORS-thLSSX} in the case in which there are multiple uncorrelated Brownian motions. To illustrate this possibility, consider the system described by
\begin{equation}
\label{ORS-eq1u}
dx_t=[Ax_t+Bu_t]dt+[Fx_t+Gu_t]d\W_t^x
\end{equation}
and the generator
\begin{equation}
\label{ORS-eq2u}
d\omega_t= S\omega_t dt+J\omega_t d\W_t^s,
\end{equation}
where $\W_t^x$ and $\W_t^s$ are two uncorrelated Brownian motions, \textit{i.e.} $d\W_t^x d\W_t^s=0$.

\begin{corollary}
\label{ORS-coLSS}
Consider system~(\ref{ORS-eq1u}) driven by the signal generator~(\ref{ORS-eq2u}). The process $\X_t$ in Theorem~\ref{ORS-thLSSX} is now the steady-state solution of
\begin{equation}
\label{MRSS-eqLSSu}
\begin{array}{rl}
\!\!\!d\X_t\!=\!\!\!\!\!&\left(A\X_t-\X_t\left(S-J^2\right)+BL\right)dt+\left(F\X_t+GL\right)d\W_t^x\\[2mm]
& -\X_t J d\W_t^s.
\end{array}
\end{equation}
\end{corollary}

\textit{Proof:} We focus on the derivation of equation~(\ref{MRSS-eqLSSu}). For now, assume that the steady state of system~(\ref{ORS-eq1u}) driven by the signal generator~(\ref{ORS-eq2u}) can be written as $x_t^{ss}=\X_t \omega_t$ for some $\X_t$. We use the stochastic product rule to compute $d(\X_t\omega_t)$ yielding
$$
\begin{array}{rl}
d(\X_t\omega_t)&\!\!\!\!=d\X_t\omega_t+\X_td\omega_t+d\X_td\omega_t\\
&\!\!\!\!=d\X_t\omega_t+\X_t(S\omega_t dt + J \omega_t d\W_t^s) + d\X_td\omega_t.
\end{array}
$$
At steady state (\textit{i.e.} $x_t^{ss}=\X_t \omega_t$) this last equation, (\ref{ORS-eq1u}) and (\ref{ORS-eq2u}) give
$$
\begin{array}{rl}
\!\!\!d\X_t\omega_t+d\X_td\omega_t&\!\!\!\!=\left(A\X_t-\X_t S+BL\right)\omega_tdt\\
&+\left(F\X_t+GL\right)\omega_td\W_t^x-\X_t J \omega_td\W_t^s.
\end{array}
$$
The only non-zero term in $d\X_td\omega_t$ is the one resulting from the product of the components of $d\X_t$ and $d\omega_t$ which multiply $d\W_t^s$. Hence, we easily see that $d\X_td\omega_t=-\X_t J^2 \omega_t dt$, from which~(\ref{MRSS-eqLSSu}) follows. Repeating the steps of the proof of Theorem~\ref{ORS-thLSSX}, \textit{i.e.} defining $z_t:=x_t - \X_t\omega_t$ and computing $dz_t$, it is straightforward to show that the relation $x_t^{ss}=\X_t \omega_t$ indeed holds. 
\hspace*{0pt}\hfill $\square$ \par

Combining the results of Theorem~\ref{ORS-thLSSX} and Corollary~\ref{ORS-coLSS}, it is straightforward to generalize the theory to an arbitrary number of identical, correlated or uncorrelated Brownian motions.

\section{reduced-order models}
\label{MRNS-secROM}

In this section we provide families of linear and nonlinear reduced-order models solving Problem~\ref{MRSto-pro}. We present first the results for nonlinear systems and then for linear systems.

\subsection{Nonlinear systems}

In accordance with Problem~\ref{MRSto-pro}
we provide the definition of (reduced) model of system~(\ref{MRNS-eq1}) at $(s,j,l)$.
\begin{definition}
\label{MRNS-defROM}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). The system described by equation~(\ref{MRSto-eqRM})
is a \textit{stochastic model of system~(\ref{MRNS-eq1}) at $(s,j,l)$} if system~(\ref{MRSto-eqRM}) has the same moment at $(s,j,l)$ of system~(\ref{MRNS-eq1}). System~(\ref{MRSto-eqRM}) is a \textit{stochastic reduced-order model of system~(\ref{MRNS-eq1}) at $(s,j,l)$} if $\nu<n$.
\end{definition}

From this definition a result follows straightforwardly.
\begin{lemma}
\label{MRNS-lmROM1}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). Suppose Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} hold. Then system~(\ref{MRSto-eqRM}) is a stochastic model of system~(\ref{MRNS-eq1}) at $(s,j,l)$ if there exists a stochastic process $\rho:(t,\W_t,\omega_t)\mapsto \rho_t(\W_t,\omega_t)$, locally defined in a neighborhood $\red{W}(w) \subset E_c(w)$ of $\omega_t=0$, with $\rho_t(0,0)=0$, which satisfies the equation
\begin{equation}
\label{MRNS-eqP}
\begin{array}{l}
d \rho_t= \red{f}(\rho_t(\W_t,\omega_t),l(\omega_t))dt + \red{g}(\rho_t(\W_t,\omega_t),l(\omega_t))d\W_t,
\end{array}
\end{equation}
and it is such that for all $t\ge0$
\begin{equation}
\label{MRNS-eqCP}
h(\chi_t(\W_t,\omega_t))=\red{h}(\rho_t(\W_t,\omega_t)),
\end{equation}
almost surely, where $\chi_t$ is a solution of (\ref{MRNS-eqNSS}).
\end{lemma}

\textit{Proof:}
If $\rho_t$ is a solution of equation~(\ref{MRNS-eqP}), then the moment of system~(\ref{MRSto-eqRM}) at $(s,j,l)$ is by definition the mapping $\red{h} \circ \rho_t$. Equation~(\ref{MRNS-eqCP}) imposes that the moment of system~(\ref{MRSto-eqRM}) is equal to the moment of system~(\ref{MRNS-eq1}). \hspace*{0pt}\hfill $\square$ \par

Problem~\ref{MRSto-pro} can be now reformulated as the problem of determining the mappings $\red{f}$, $\red{g}$ and $\red{h}$ in Lemma~\ref{MRNS-lmROM1} such that the two equations~(\ref{MRNS-eqP}) and (\ref{MRNS-eqCP}) are satisfied. This problem can be solved easily.


\begin{proposition}
\label{MRNS-lmROM2}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). Suppose Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} hold. Then the system
\begin{equation}
\label{MRNS-eqROM2}
\begin{array}{rl}
d \red{x}_t &\!\!\!\!=(s(\red{x}_t) -\delta(\red{x}_t)l(\red{x}_t) + \delta(\red{x}_t) u_t)dt \\[2mm]
&+(j(\red{x}_t)-\eta(\red{x}_t)l(\red{x}_t)+\eta(\red{x}_t)u_t)d\W_t, \\[2mm]
\red{y}&\!\!\!\!=h(\chi_t(\W_t,\red{x}_t)),
\end{array}
\end{equation}
where $\chi_t$ is a solution of (\ref{MRNS-eqNSS}), is a stochastic model of system~(\ref{MRNS-eq1}) at $(s,j,l)$ if $\delta$ and $\eta$ are arbitrary mappings such that equation
\begin{equation}
\label{MRNS-eqrho111}
\begin{array}{l}
d\rho_t = \left(s(\rho_t) - \delta(\rho_t)l(\rho_t) +\delta(\rho_t)l(\omega_t)\right) dt \\[2mm]
\qquad +
\left(j(\rho_t) -\eta(\rho_t)l(\rho_t)+\eta(\rho_t)l(\omega_t)\right)d\W_t,
\end{array}
\end{equation}
has the trivial solution $\rho_t(\W_t,\omega_t)=\omega_t$.
\end{proposition}

\textit{Proof:}
In model~(\ref{MRSto-eqRM}) we select $\red{f}$ and $\red{g}$ such that the system is affine in the input, namely $\red{f}(\red{x}_t,u_t)=\phi(\red{x}_t)+\delta(\red{x}_t) u_t$ and $\red{g}(\red{x}_t,u_t)=\gamma(\red{x}_t)+\eta(\red{x}_t) u_t$ for some smooth mapping $\phi$, $\delta$, $\gamma$ and $\eta$ to be determined. Selecting $\rho_t(\W_t,\omega_t)=\omega_t$, \textit{i.e.} the identity mapping, yields (by writing (\ref{MRNS-eqSPDE}) for the reduced-order model) $d\rho_t = s(\omega_t) dt + j(\omega_t)d\W_t$. Hence, equation~(\ref{MRNS-eqP}) becomes
$$
\begin{array}{rl}
s(\omega_t) dt + j(\omega_t)d\W_t&\!\!\!\!= \left( \phi(\omega_t)+\delta(\omega_t) l(\omega_t)\right) dt \\
&+ \left( \gamma(\omega_t)+\eta(\omega_t) l(\omega_t)\right) d\W_t
\end{array}
$$
which is satisfied selecting
$$
\phi(\omega_t)=s(\omega_t)-\delta(\omega_t) l(\omega_t), \qquad
\gamma(\omega_t)=j(\omega_t)-\eta(\omega_t) l(\omega_t),
$$
\textit{i.e.} (\ref{MRNS-eqrho111}) is satisfied.
Thus, selecting $\phi(\red{x}_t) = \left. \phi(\omega_t)\right|_{\omega_t=\rho^{-1}(\red{x}_t)}$ and $\gamma(\red{x}_t) = \left. \gamma(\omega_t)\right|_{\omega_t=\rho^{-1}(\red{x}_t)}$,
model~(\ref{MRNS-eqROM2}) satisfies equations~(\ref{MRNS-eqP}) and (\ref{MRNS-eqCP}) for any mapping $\delta$ and $\eta$ such that~(\ref{MRNS-eqrho111})  has the unique solution $\rho_t(\W_t,\omega_t)=\omega_t$.
\hspace*{0pt}\hfill $\square$ \par

The family of models (\ref{MRNS-eqROM2}) is parametrized by the mappings $\delta$ and $\eta$. These mappings can be used to span the family of reduced-order models by moment matching so that specific additional properties are imposed. For instance, we may want to preserve stability or we may want to achieve a
special representation of the reduced-order model which is particularly
useful for a desired application. 

Note that in Definition~\ref{MRSS-DefMom} the moment of system~(\ref{MRNS-eq1}) is defined as the solution of an equation, namely~(\ref{MRNS-eqNSS}), without any reference to the steady state of the system. In fact, the moment may exist (when equation~(\ref{MRNS-eqNSS}) has a solution) even though Assumption~\ref{MRNS-as1} is not satisfied and we cannot identify the mapping $\chi_t$ as a steady state. In fact, Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} are sufficient to guarantee that equation~(\ref{MRNS-eqNSS}) has a solution, but not necessary. This fact explains why asymptotic stability of the origin of the reduced-order model is not required by Lemma~\ref{MRNS-lmROM1} and Proposition~\ref{MRNS-lmROM2}. However, one can use the mappings $\delta$ and $\eta$ to impose stability and simplify the result of Proposition~\ref{MRNS-lmROM2}.

\begin{corollary}
\label{MRNS-cor1}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}). Suppose Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} hold. Then system~(\ref{MRNS-eqROM2}) is a stochastic model of system~(\ref{MRNS-eq1}) at $(s,j,l)$ for any $\delta$ and $\eta$ such that the linearization of (\ref{MRNS-eqROM2}) around zero has all negative Lyapunov exponents almost surely.
\end{corollary}

\textit{Proof:}
This is a direct consequence of Proposition~\ref{MRNS-lmROM2} and Theorem~\ref{MRNS-thm1}.
\hspace*{0pt}\hfill $\square$ \par

Note that there exist multiple $\delta$ and $\eta$ satisfying the condition in Corollary~\ref{MRNS-cor1}. Thus, $\delta$ and $\eta$ are still partially free to be used to impose additional properties (\textit{e.g.} different speeds of decay of the transient among all stabilizing mappings $\delta$ and $\eta$).


\subsection{Linear systems}

This section has the same structure as the previous section. We begin by providing the definition of model of system~(\ref{ORS-eq1}) at $(S,J,L)$.
\begin{definition}
\label{MRLS-defROM}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). The system described by the equations
\begin{equation}
\label{MRLS-eqROM}
d \red{x}_t =(\red{A} \red{x}_t + \red{B} u_t)dt + (\red{F} \red{x}_t +\red{G} u_t)d\W_t, \quad \red{y}_t=\red{C} \red{x}_t,
\end{equation}
where $\red{A}\in \mathbb{R}^{\nu\times \nu}$, $\red{B}\in \mathbb{R}^{\nu\times1}$, $\red{F}\in \mathbb{R}^{\nu\times \nu}$, $\red{G}\in \mathbb{R}^{\nu\times1}$, $\red{C}\in \mathbb{R}^{1\times \nu}$, is a \textit{stochastic model of system~(\ref{ORS-eq1}) at $(S,J,L)$}, if system~(\ref{MRLS-eqROM}) has the same moments of system~(\ref{ORS-eq1}) at $(S,J,L)$. System~(\ref{MRLS-eqROM}) is a \textit{stochastic reduced-order model of system~(\ref{ORS-eq1}) at $(S,J,L)$} if $\nu<n$.
\end{definition}

The equivalent of Lemma~\ref{MRNS-lmROM1} follows straightforwardly.
\begin{lemma}
\label{MRLS-lmROM1}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). Suppose Assumptions~\ref{ORS-asSys} and \ref{ORS-asGen} hold. Then system~(\ref{MRLS-eqROM}) is a stochastic model of system~(\ref{ORS-eq1}) at $(S,J,L)$ if there exists $\cR_t\in\mathbb{R}^{\nu \times \nu}$
which satisfies the equation
\begin{equation}
\label{MRLS-eqP}
\begin{array}{rl}
d\cR_t=\!\!\!\!&\left(\red{A}\cR_t-\cR_t\left(S-J^2\right)-\red{F}\cR_t J+\red{B}L - \red{G}LJ \right)dt\\[2mm]
&+(\red{F}\cR_t-\cR_tJ+\red{G}L)d\W_t,
\end{array}
\end{equation}
and it is such that for all $t\ge 0$
\begin{equation}
\label{MRLS-eqCP}
C\X_t=\red{C} \cR_t,
\end{equation}
almost surely, where $\X_t$ is a solution of (\ref{ORS-eqPI}).
\end{lemma}

\textit{Proof:} The proof is analogous to the proof of Lemma~\ref{MRNS-lmROM1} and thus is omitted.
\hspace*{0pt}\hfill $\square$ \par

In the spirit of Proposition~\ref{MRNS-lmROM2} we select the matrices of the model to satisfy equations~(\ref{MRLS-eqP}) and (\ref{MRLS-eqCP}). 

\begin{proposition}
\label{MRLS-lmROM2}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). Suppose Assumptions~\ref{ORS-asSys} and \ref{ORS-asGen} hold. Then the system
\begin{equation}
\label{MRLS-eqROM2}
\begin{array}{rl}
d \red{x}_t &\!\!\!\!=((S-\red{B}L) \red{x}_t + \red{B} u_t)dt + ((J-\red{G} L) \red{x}_t + \red{G} u_t)d\W_t, \\[2mm]
\red{y}_t&\!\!\!\!=C\X_t \red{x}_t,
\end{array}
\end{equation}
where $\X_t$ is a solution of (\ref{ORS-eqPI}), is a stochastic model of system~(\ref{ORS-eq1}) at $(S,J,L)$ for any $\red{B}$ and $\red{G}$ such that equation~(\ref{MRLS-eqP}) has the unique solution $\cR_t = I$.
\end{proposition}

\textit{Proof:} Equation~(\ref{MRLS-eqP}) has solution $\cR_t = I$ only if
$$
\begin{array}{l}
0 = \red{A}-\left(S-J^2\right)-\red{F}J+\red{B}L - \red{G}LJ,\\
0 =\red{F}-J+\red{G}L.
\end{array}
$$
Solving the second equation with respect to $\red{F}$ yields $\red{F}=J-\red{G}L$. Replacing this expression in the first equation and solving with respect to $\red{A}$ yields $\red{A}=S-\red{B}L$. Finally, substituting $\cR_t = I$ in equation~(\ref{MRLS-eqCP}) yields $\red{C} =C\X_t$.
\hspace*{0pt}\hfill $\square$ \par

Note that the matrices $\red{B}$ and $\red{G}$ play the same role that $\delta$ and $\eta$ have in the nonlinear case. Thus $\red{B}$ and $\red{G}$ can be freely selected to achieve additional properties for the reduced-order model.

Note that there is no loss of generality in selecting the solution of~(\ref{MRLS-eqP}) as in Proposition~\ref{MRLS-lmROM2}. In fact, if a model~(\ref{MRLS-eqROM}) satisfies Lemma~\ref{MRLS-lmROM1} for $\cR_t \ne I$, we can define a stochastic change of coordinates for which the model in the new coordinates has the form~(\ref{MRLS-eqROM2}). We formalize this property in the next result.
\begin{lemma}
Assume that system~(\ref{MRLS-eqROM}) is a reduced-order model satisfying Lemma~\ref{MRLS-lmROM1}. Let $\xi_t:= \cR_t^{-1}\red{x}_t$. Then the system
\begin{equation}
\label{MRLS-eqROM2bis}
\begin{array}{rl}
d \xi_t &\!\!\!\!=((S-\overline{B}L) \xi_t + \overline{B} u_t)dt + ((J-\overline{G} L) \xi_t + \overline{G} u_t)d\W_t, \\[2mm]
\red{y}_t&\!\!\!\!=C\X_t \xi_t,
\end{array}
\end{equation}
with $\overline{G}=\cR_t^{-1}\red{G}$ and $\overline{B}=\cR_t^{-1}(\red{B}-(\red{F}\cR_t-\cR_tJ+\red{G}L)\overline{G})$ satisfies Proposition~\ref{MRLS-lmROM2}.
\end{lemma}

\textit{Proof:} Note that $d\red{x}_t=d\cR_t \xi_t +\cR_t d\xi_t+d\cR_t d\xi_t$ and $d\xi_t=\cR_t^{-1}\left((\red{A} \red{x}_t + \red{B} u_t)dt + (\red{F} \red{x}_t +\red{G} u_t)d\W_t-d\cR_t (\xi_t+ d\xi_t)\right)$. The result follows from tedious but straightforward computation.
\hspace*{0pt}\hfill $\square$ \par

\section{Classes of approximated reduced-order models}
\label{MRSS-sec-AROM}

Analysing the models proposed in Propositions~\ref{MRNS-lmROM2} and \ref{MRLS-lmROM2}, we note that the method relies on the determination of the mappings $\chi_t$ and $\X_t$, respectively. Currently, the determination of these mappings requires the solution of equations~(\ref{ORS-eqPI}) and (\ref{MRNS-eqNSS}), respectively, which is computationally expensive. For instance, equation~(\ref{ORS-eqPI}) is a stochastic matrix equation which consists of $n\nu$ linear stochastic equations. Thus, to determine $\red{y}_t$ we need to construct $n\nu$ stochastic processes, \textit{i.e.} the components of $\X_t$. This would not be a problem if $\X_t$ could be determined off-line, but since $\X_t$ depends on the Brownian motion, this is not possible. Hence, while the models proposed in Proposition~\ref{MRLS-lmROM2} possess the same moment or, equivalently, the same steady-state output response of system~(\ref{ORS-eq1}), they can be hardly considered ``simpler''. Note that this issue could eventually be solved. In fact, on one hand what we need is $C\X_t$, which is of order $\nu$, rather than $\X_t$. Thus, the development of efficient ways of determining the moment $C\X_t$ without computing $\X_t$ may solve the issue (this is normally achieved in the deterministic case, see \textit{e.g.} \cite{Ant:05}). On the other hand, it may even be possible to develop methods to determine directly models that match the moment without computing it (as achieved by some deterministic techniques, see \textit{e.g.} \cite{GugAntBea:08}).

In this paper we propose an alternative way to overcome this difficulty. In particular, the idea is to determine stochastic reduced-order models that, maintaining a subset of the stochastic properties of the system to be reduced, allows us to carry out off-line all computations that have a complexity which depends on $n$, thus providing a computational advantage. In the following we propose classes of ``approximated'' reduced-order models and we discuss the properties of these models. For simplicity, the rest of the paper focuses on linear systems, although nonlinear extensions are briefly mentioned.

\subsection{Reduced-order models preserving the mean of the moment $E[C\X_t]$}

The first class of models that we consider are models obtained considering a relaxation of Problem~\ref{MRSto-pro} in which the limit~(\ref{MRSto-eqMMc}), namely
$$
\lim_{t\to \infty} C\X_t\omega_t - \red{C}\cR_t\omega_t = 0
$$
is replaced by
\begin{equation}
\label{MRLS-eqCe2}
\lim_{t\to \infty} CE[\X_t]\omega_t - \red{C}E[\cR_t]\omega_t = 0.
\end{equation}
In other words, we approximate $\X_t$ with its expectation and, instead of matching the moment, we match just its expectation. To make the concept precise, we introduce the following definition.

\begin{definition}
\label{MRLS-defROMm}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). The system described by equation~(\ref{MRLS-eqROM}) is a \textit{stochastic model in the moment-mean of system~(\ref{ORS-eq1}) at $(S,J,L)$}, if the mean of the moment of system~(\ref{MRLS-eqROM}) is equal to the mean of the moment of system~(\ref{ORS-eq1}) at $(S,J,L)$. System~(\ref{MRLS-eqROM}) is a \textit{stochastic reduced-order model in the moment-mean of system~(\ref{ORS-eq1}) at $(S,J,L)$} if $\nu<n$.
\end{definition}

To determine a class of models that satisfy Definition~\ref{MRLS-defROMm}, let $\Pi= E[\X_t]$ and note, from~(\ref{ORS-eqPI}), that $\Pi$ obeys the equation
\begin{equation}
\label{MRSS-eqEXi}
d\Pi= \left(A\Pi-\Pi\left(S-J^2\right)-F\Pi J+BL-GLJ\right)dt.
\end{equation}
The equilibrium point of this equation is given by solving the generalized Sylvester equation\footnote{This equation has a unique solution if and only if $0 \not \in \sigma\left(I\otimes A - (S-J^2)^{\top} \otimes I - J^{\top} \otimes F \right)$.}
\begin{equation}
\label{MRSS-eqEPI}
A\Pi-\Pi\left(S-J^2\right)-F\Pi J+BL-GLJ = 0.
\end{equation}
The equilibrium point is unique and describes the steady-state solution of~(\ref{MRSS-eqEXi}) if $A$, $F$, $S$ and $J$ are such that the equilibrium of (\ref{MRSS-eqEXi}) is asymptotically stable.\\
A family of reduced-order models in the moment-mean follows.

\begin{proposition}
\label{MRLS-lmROM3}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). Suppose Assumptions~\ref{ORS-asSys} and \ref{ORS-asGen} hold and that equation~(\ref{MRSS-eqEXi}) has an asymptotically stable equilibrium point $\Pi$ solving equation~(\ref{MRSS-eqEPI}). Then the system
\begin{equation}
\label{MRLS-eqROM3}
\begin{array}{rl}
d \red{x}_t &\!\!\!\!=((S-\red{B}L) \red{x}_t + \red{B} u_t)dt + (\red{F} \red{x}_t + \red{G} u_t)d\W_t, \\[2mm]
\red{y}_t&\!\!\!\!=C\Pi \red{x}_t,
\end{array}
\end{equation}
is a stochastic model in the moment-mean of system~(\ref{ORS-eq1}) at $(S,J,L)$, for any $\red{G}$ and for any $\red{B}$ and $\red{F}$ such that 
\begin{equation}
\label{MRLS-eqCER}
\sigma\left(I \otimes \left(S-\red{B}L\right) - \left(S-J^2\right)^{\top} \otimes I - J^{\top} \otimes \red{F} \right) \subset \mathbb{C}_{<0}.
\end{equation}
\end{proposition}

\textit{Proof:} Consider system~(\ref{MRLS-eqROM}) and the equation~(\ref{MRLS-eqP}) defining $\cR_t$. Let $R:=E[\cR_t]$ and note that $R$ satisfies
\begin{equation}
\label{MRLS-eqER}
dR=\left(\red{A}R-R\left(S-J^2\right)-\red{F}R J+\red{B}L - \red{G}LJ \right)dt.
\end{equation}
Let $\red{A}=S-\red{B}L$. If $\red{B}$ and $\red{F}$ are such that (\ref{MRLS-eqCER}) holds, then $R=I$ is the unique attractive equilibrium of equation~(\ref{MRLS-eqER}). Finally, the selection $\red{C}=C\Pi$ satisfies condition~(\ref{MRLS-eqCe2}).
\hspace*{0pt}\hfill $\square$ \par

\begin{remark}
The advantage of reduced-order models in the moment-mean is that we need to determine only $C\Pi$, which can be computed with a plethora of efficient methods (see \cite{Ant:05} for a review and \textit{e.g.} \cite{ScaAst:16a}).
\end{remark}

In the next simple example we illustrate the differences between a reduced-order model~(\ref{MRLS-eqROM2}) and a reduced-order model in the moment-mean~(\ref{MRLS-eqROM3}).

\begin{example}
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/momentmean2.eps}
\caption{Top graph: time history of the output of system~(\ref{ORS-eq1}) (solid/blue), of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) (dashed/red) and of the output of the stochastic reduced-order model in the moment-mean~(\ref{MRLS-eqROM3}) (dash-dotted/black). Bottom graph: time history of the corresponding absolute errors.}
\label{MRSS-fig1}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth]{imgs/MultipleRealizations.eps}
\caption{Six different realizations of the simulation in Figure~\ref{MRSS-fig1}. For each realization, the top set of curves are the output trajectories and the bottom set are the corresponding absolute errors.}
\label{MRSS-fig1bis}
\end{figure}%
Consider a linear system~(\ref{ORS-eq1}) of order $n=200$ with randomly generated matrices\footnote{All the matrices used in this simulation can be downloaded from \cite{FILE_MRSS}.} $A$, $B$ and $C$, and with $F=0.05A$ and $G=0.1B$. Consider a signal generator~(\ref{ORS-eq2}) of order $\nu=2$ with $J$ a randomly generated such that $S=[0\,5;-5\,0]+0.5J^2$ and $J$ commute (thus Assumption~\ref{ORS-asGen} is verified by construction). A stochastic reduced-order model~(\ref{MRLS-eqROM2}) and a stochastic reduced-order model in the moment-mean~(\ref{MRLS-eqROM3}) are computed. For both models we have selected $\red{B}$ and $\red{G}$ such that the eigenvalues of $\red{A}=S-\red{B}L$ and $\red{F}=J-\red{G}L$ are two of the eigenvalues of $A$ and $F$, respectively. The top graph in Fig.~\ref{MRSS-fig1} shows the output of the system (solid/blue), of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) (dashed/red) and of the stochastic reduced-order model in the moment-mean~(\ref{MRLS-eqROM3}) (dash-dotted/black). The bottom graph shows the corresponding absolute errors. We notice that both reduced-order models approximate the output of the system, although the model in the moment-mean shows a steady-state mismatch
Fig.~\ref{MRSS-fig1bis} shows the same quantities for six different realizations. The top set of lines correspond to the top graph in Fig.~\ref{MRSS-fig1}, whereas the bottom set of lines correspond to the bottom graph in Fig.~\ref{MRSS-fig1}. The figure shows self-consistence across multiple realizations of the noise.
\end{example}

\subsection{Reduced-order models preserving the mean and mean-square of the steady state when $J=0$}

In this section we clarify the relation between the stochastic models in the moment-mean defined in the previous section and families of approximated stochastic models introduced in \cite{ScaTee:17,ScaTee:17a}. Subsequently, we provide additional results on the latter.

\subsubsection{Relation between models in the moment-mean and models in the mean introduced in \cite{ScaTee:17,ScaTee:17a}}

Since $\omega_t$ is a stochastic process that depends on the same noise of $\X_t$, we have that $E[\X_t\omega_t] \ne E[\X_t]E[\omega_t] \ne E[\X_t]\omega_t$. Hence, the reduced-order models in the moment-mean do not necessarily have the same steady-state mean of the output of the system to be reduced. Note, however, that if $J=0$, then $\omega_t$ is deterministic and condition~(\ref{MRLS-eqCe2}) can be written as
\begin{equation}
\label{MRLS-eqLEe}
\lim_{t\to \infty} E[e_t] = 0.
\end{equation}
Thus, in this case matching the mean of the moment corresponds to matching the mean of the steady-state output. Note also that the mean $E[C\X_t]=C\Pi$, where $\Pi$ is the unique solution of 
\begin{equation}
\label{MRLS-eqSyl}
A\Pi-\Pi S + BL = 0,
\end{equation}
is the moment of the deterministic system
$$
\dot m=Am+Bu,\qquad y=Cm,
$$
where $m=E[x_t]$. Thus, the stochastic models in the moment-mean defined in the previous section are a generalisation of the ``stochastic models in the mean'' introduced in \cite{ScaTee:17}, which were based on the idea of preserving the mean of the steady-state output of system~(\ref{ORS-eq1}) driven by~(\ref{ORS-eq2}), \textit{i.e.} satisfying~(\ref{MRLS-eqLEe}).

Note that, even though the moment $C\Pi$ is the moment of a deterministic system, the models in the moment-mean (and the models in the mean defined in \cite{ScaTee:17}) are in general stochastic models. 

\begin{remark}
\label{MRSS-rm7}
Consider the special selection $\red{F}=-\red{G}L$ in (\ref{MRLS-eqROM3}). The steady state of the resulting model for $u=L\omega$ is generated by a deterministic system. In fact, the equation
$$
d\cR_t=((S-\red{B}L)\cR_t-\cR_t S+\red{B}L)dt+(-\red{G}L\cR_t+\red{G}L)d\W_t,
$$
which is solved by $\cR_t=I$, is independent of $d\W_t$. Thus, at steady state we have that $-\red{G}L \red{x}^{ss}_t +\red{G}L \omega_t = -\red{G}L \cR_t \omega_t +\red{G}L \omega_t= 0$. As a result, the steady-state system is not stochastic. We stress, however, that this is not the case if $F_r\ne -G_rL$.
\end{remark}

Similarly, stochastic models in the mean can be defined also for nonlinear systems as those stochastic models that have the same mean of the steady-state output of system~(\ref{MRNS-eq1}) driven by~(\ref{MRNS-eq2}). A family of nonlinear stochastic models in the mean is given in the next result.

\begin{proposition}
\label{MRNS-lmROMNLMean}
Consider system~(\ref{MRNS-eq1}) and the signal generator~(\ref{MRNS-eq2}) with $j(\cdot)\equiv 0$. Suppose Assumptions~\ref{MRNS-as1} and \ref{MRNS-as2} hold and that $\sigma\left(\left.\frac{\partial f}{\partial x_t}\right|_{\substack{x_t=0\\\omega=0}}\right) \subset \mathbb{C}_{<0}$. Then the system
\begin{equation}
\label{MRNS-eqROM2Mean}
\begin{array}{rl}
d \red{x}_t &\!\!\!\!=(s(\red{x}_t) -\delta(\red{x}_t)l(\red{x}_t) + \delta(\red{x}_t) u)dt +\gamma(\red{x}_t,u)d\W_t, \\[2mm]
\red{y}&\!\!\!\!=h(\pi(\red{x}_t)),
\end{array}
\end{equation}
where $\pi$ is a solution of (\ref{MRSSeqISI}), is a stochastic model in the mean of system~(\ref{MRNS-eq1}) at $(s,l)$ for any $\gamma$ and for any $\delta$ such that the zero equilibrium of $d \red{x} =s(\red{x}) -\delta(\red{x})l(\red{x})$ is locally exponentially stable.
\end{proposition}

\textit{Proof:} The mean of system~(\ref{MRNS-eq1}) is described by $\dot{\overbrace{E[x_t]}} = f(E[x_t],u)$. This is a nonlinear deterministic system which has moment $h\circ \pi$, where $\pi$ is the solution of the deterministic partial differential equation~(\ref{MRSSeqISI}). Note that such solution exists by Assumption~\ref{MRNS-as2} and the assumption on the Jacobian of $f$, see \cite{Isi:95}. Then the claim follows from standard deterministic model reduction \cite{Ast:10}.
\hspace*{0pt}\hfill $\square$ \par

\subsubsection{Stochastic models in the mean-square}
Although the models in the mean are stochastic systems, the matrices $\red{F}$ and $\red{G}$ in~(\ref{MRLS-eqROM3}) are free parameters and do not preserve, in a systematic way, information regarding the matrices $F$ and $G$ of the system to be reduced. Now, keeping the standing assumption that $J=0$, we want to ``improve'' the models in the mean. In particular, we want to use the parameters $\red{F}$ and $\red{G}$ to preserve the mean-square, in addition to the expectation, of the steady-state output of system~(\ref{ORS-eq1}) driven by~(\ref{ORS-eq2}). To preserve also this information we propose another class of reduced-order models that solve a variation of Problem~\ref{MRSto-pro} in which the limit~(\ref{MRSto-eqMMc}) is replaced by
\begin{equation}
\label{MRSS-Verr}
\lim_{t\to \infty} E[e_t] = 0 \quad \text{and} \quad \lim_{t\to \infty} E[e_te_t^{\top}] = 0.
\end{equation}
To make the concept precise, we introduce the following definition.
\begin{definition}
\label{MRLS-defROMV}
Consider system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). The system described by equation~(\ref{MRLS-eqROM}) is a \textit{stochastic model in the mean-square of system~(\ref{ORS-eq1}) at $(L,S)$}, if the conditions in (\ref{MRSS-Verr}) are simultaneously satisfied. System~(\ref{MRLS-eqROM}) is a \textit{stochastic reduced-order model in the mean-square of system~(\ref{ORS-eq1}) at $(L,S)$} if $\nu<n$.
\end{definition}

We remark that we have defined the models in the ``mean-square'' in such a way that in addition to preserving the mean-square of the steady state of system~(\ref{ORS-eq1}) driven by~(\ref{ORS-eq2}), these models preserve also its mean. To the end of determining reduced-order models in the mean-square we need a preliminary result, namely a description of the steady state of $M=E[x_t x_t^{\top}]$, which obeys the equation \cite[Theorem 4.5]{Gar:88}
\begin{equation}
\label{MRLS-eqVar}
\begin{array}{rl}
\dot M&=AM+MA^{\top}+FMF^{\top}+Bum^{\top}+m(Bu)^{\top}\\[2mm]
&\,\,\,\,\,\,+Fm(Gu)^{\top}+Gu(Fm)^{\top}+Gu(Gu)^{\top}.
\end{array}
\end{equation}

\begin{lemma}
\label{MRLS-lmSM}
Consider the interconnection of system~(\ref{MRLS-eqVar}) and the signal generator~(\ref{ORS-eq2}). Let $\mathcal{A}=I\otimes A + A\otimes I + F \otimes F$ and assume $\textstyle\sigma\left(\mathcal{A}\right)\subset\mathbb{C}_{<0}$ and $\sigma(S)\subset \mathbb{C}_{0}$. Then the steady-state response of such interconnection is
$$
\vect{M^{ss}}=\mathcal{K} \vect{\omega\omega^{\top}},
$$
where $\mathcal{K}$ is the unique solution of the augmented Sylvester equation
\begin{equation}
\label{MRLS-eqMSS}
\mathcal{A} \mathcal{K} + \mathcal{B}= \mathcal{K} \mathcal{S},
\end{equation}
with $\mathcal{S}=I\otimes S + S\otimes I$ and $\mathcal{B}=BL \otimes \Pi+ \Pi \otimes BL + GL \otimes F\Pi + F\Pi \otimes GL + GL \otimes GL$.
\end{lemma}

\textit{Proof.}
First of all note that $\mathcal{K}$ is the unique solution of the Sylvester equation~(\ref{MRLS-eqMSS}) because $\sigma(\mathcal{A})\cap\sigma(\mathcal{S})=\emptyset$. Let $\mu:=\vect{M}-\mathcal{K}\vect{\omega\omega^{\top}}$ and compute the derivative of $\mu$ with respect to time. Moreover, substitute $m$ with $\Pi\omega$, which is the steady state of $m$. Using the vectorization operator and the Kronecker product yields
$$
\dot \mu= \mathcal{A}\mu + (\mathcal{A} \mathcal{K} - \mathcal{K} \mathcal{S} + \mathcal{B})\vect{\omega\omega^{\top}}=\mathcal{A}\mu.
$$
Since $\textstyle\sigma\left(\mathcal{A}\right)\subset\mathbb{C}_{<0}$, $\mu$ converges exponentially to zero and $\vect{M}$ converges exponentially to $\mathcal{K}\vect{\omega\omega^{\top}}$.
\hspace*{0pt}\hfill $\square$ \par

We are now ready to give a family of reduced-order models in the mean-square.

\begin{proposition}
\label{MRLS-lmROM4}
Consider the interconnection of system~(\ref{ORS-eq1}) and the signal generator~(\ref{ORS-eq2}). Assume $\textstyle\sigma\left(\mathcal{A}\right)\subset\mathbb{C}_{<0}$ and $\sigma(S)\subset \mathbb{C}_{0}$. Assume there exists a matrix $\red{C}$ such that 
\begin{equation}
\label{MRSS-eqC1M}
\red{C} \otimes \red{C} = \left( C \otimes C \right) \mathcal{K}, 
\end{equation}
where $\mathcal{K}$ is the unique solution of~(\ref{MRLS-eqMSS}). Let $R$ be any invertible matrix such that $\red{C}R=C\Pi$, where $\Pi$ is the unique solution of~(\ref{MRLS-eqSyl}). Let
$\red{A}=RSR^{-1}-\red{B}LR^{-1}$ for any $\red{B}$ such that $\sigma(\red{A})\subset \mathbb{C}_{<0}$. Assume there exist matrices $\red{F}$ and $\red{G}$ such that
\begin{equation}
\label{MRLS-eqFrGr}
\begin{array}{l}
\!\red{F}\! \otimes\! \red{F}+\red{G}L \!\otimes\! \red{F} R + \red{F} R\! \otimes\! \red{G}L + \red{G} L \!\otimes\! \red{G} L=\\[2mm]
\,\,\,\,\,\,=-I\!\otimes\! \red{A} - \red{A}\!\otimes\! I+I\!\otimes \!S + S\!\otimes\! I - \red{B}L \!\otimes\! R- R\! \otimes\! \red{B}L,
\end{array}
\end{equation}
and that $\textstyle\sigma\left(I\otimes \red{A} + \red{A}\otimes I + \red{F} \otimes \red{F}\right)\subset\mathbb{C}_{<0}$. Then the system
\begin{equation}
\label{MRLS-eqROM4}
\begin{array}{l}
d \red{x}_t =(\red{A} \red{x}_t + \red{B} u_t)d\xi + (\red{F} \red{x}_t + \red{G} u_t)d\W_t, \\[2mm]
\red{y}_t=\red{C} \red{x}_t,
\end{array}
\end{equation}
is a stochastic model in the mean-square of system~(\ref{ORS-eq1}) at $(S,L)$.
\end{proposition}

\textit{Proof.} Under certain stability properties (that will be established later in the proof), the steady state of $\red{m}=E[\red{x}_t]$ (\textit{i.e.} of the mean of $\red{x}_t$) is $R\omega$, with $R$ the unique solution of  
\begin{equation}
\label{MRLS-eqmROM}
\red{A} R - R S = -\red{B} L.
\end{equation}
Similarly, the steady state of $\vect{\red{M}}=\vect{E[\red{x}_t\red{x}_t^{\top}]}$ is $\red{\mathcal{K}} \vect{\omega \omega^{\top}}$, where $\red{\mathcal{K}}$ is the unique solution of the equation
\begin{equation}
\label{MRLS-eqMROM}
\red{\mathcal{A}} \red{\mathcal{K}} + \red{\mathcal{B}}= \red{\mathcal{K}} \mathcal{S},
\end{equation}
with $\red{\mathcal{A}}=I\otimes \red{A} + \red{A}\otimes I + \red{F} \otimes \red{F}$, $\red{\mathcal{B}}=\red{B}L \otimes R+ R \otimes \red{B}L + \red{G}L \otimes \red{F} R + \red{F} R \otimes \red{G}L + \red{G}L \otimes \red{G}L$. Note that the mean-square of the output can be written as $\vect{E[y_ty_t^{\top}]}=\vect{CE[x_tx_t^{\top}]C^\top}=(C \otimes C)\vect{E[x_tx_t^{\top}]}$.
For system~(\ref{ORS-eq1}) and model~(\ref{MRLS-eqROM4}) to have the same steady-state mean-square of the output, we need that $(\red{C} \otimes \red{C})\red{\mathcal{K}} = (C \otimes C)\mathcal{K}$. This is achieved by setting $\red{\mathcal{K}}=I$ and $\red{C}$ such that $\red{C} \otimes \red{C} = \left( C \otimes C \right) \mathcal{K}$.
Determine now any invertible matrix $R$ such that $\red{C}R=C\Pi$. Note that since $\red{C}$ is not identically zero, it is always possible to find an invertible matrix $R$ solving this equation. Now select $\red{A}=RSR^{-1}-\red{B}LR^{-1}$ and note that this selection solves equation~(\ref{MRLS-eqmROM}). The solution $R$ is unique for any matrix $\red{B}$ such that $\sigma(\red{A})\cap \sigma(S)=\emptyset$. This last condition is guaranteed if $\red{B}$ is selected such that $\sigma(\red{A})\subset \mathbb{C}_{<0}$. Moreover, this last property also ensures that the steady state of the mean of the output of (\ref{MRLS-eqROM4}) is well-defined and equal to $\red{C}R\omega$, which by construction is equal to $C\Pi\omega$. Hence, system~(\ref{ORS-eq1}) and model~(\ref{MRLS-eqROM4}) have the same steady-state output mean. Now note that if $\red{F}$ and $\red{G}$ are selected such that~(\ref{MRLS-eqFrGr}) holds, then equation~(\ref{MRLS-eqMROM}) has the unique solution $\red{\mathcal{K}}=I$. If in addition $\red{F}$ is such that $\textstyle\sigma(\red{\mathcal{A}})\subset\mathbb{C}_{<0}$, then the steady state of the mean-square of the output of~(\ref{MRLS-eqROM4}) is $(\red{C} \otimes \red{C})\red{\mathcal{K}} \vect{\omega \omega^{\top}}$, which is equal by construction to $(C \otimes C)\mathcal{K} \vect{\omega \omega^{\top}}$. Hence, system~(\ref{ORS-eq1}) and model~(\ref{MRLS-eqROM4}) have also the same steady-state output mean-square. \hspace*{0pt}\hfill $\square$ \par

In Proposition~\ref{MRLS-lmROM4} we need to determine $\red{C}$ from condition~(\ref{MRSS-eqC1M}). This problem is known as nearest Kronecker product approximation, which can be formulated as follows. Given a matrix $Q\in\mathbb{R}^{N\times N}$, the problem consists in determining the two matrices $T_1\in\mathbb{R}^{n_1}$ and $T_2\in\mathbb{R}^{n_2}$ such that $||Q - T_1 \otimes T_2 ||$ is minimized. The solution of this problem is given in \cite{VaLPit:93,Gen:07}. After rearranging the elements of $Q$ in a new matrix called $\bar{Q}\in\mathbb{R}^{n_1^2\times n_2^2}$ (see \cite{Gen:07} for details), we compute the singular value decomposition $U^{\top} \bar{Q}V=\diag(\delta_1,\dots,\delta_q)$ of $\bar{Q}$, with $q=\rank(\bar{Q})$. The solution to the problem is given by $\vect{T_1}=\sqrt{\delta_1}u_1$ and $\vect{T_2}=\sqrt{\delta_1}v_1$, where $u_1$ and $v_1$ are the first columns of $U$ and $V$, respectively. Note that the determination of the nearest Kronecker approximation introduces an error called separability approximation error, see \cite{Gen:07}. This error is zero if the rearranged matrix has only one singular value different from zero.

\begin{remark}
The main difficulty in the determination of the family of models~(\ref{MRLS-eqROM4}) is to solve equation~(\ref{MRLS-eqFrGr}). However, since $\red{G}$ is a free parameter, we can use it to simplify the computation of such a solution. Note in fact that any $\red{G}$ such that the conditions in Proposition~\ref{MRLS-lmROM4} hold would anyway give a reduced-order model in the mean-square by moment matching according to the definition given. Thus, the free parameter $\red{G}$ can be used to achieve properties besides moment matching. In this specific remark, we use the free parameter to simplify the computation of a reduced-order model. Thus, select $\red{G}=0$. As a result, equation~(\ref{MRLS-eqFrGr}) becomes
\begin{equation}
\label{MRSS-simFF}
\red{F} \otimes \red{F}=-I\otimes \red{A} - \red{A}\otimes I+I\otimes S + S\otimes I - \red{B}L \otimes R- R \otimes \red{B}L,
\end{equation}
from which we can determine $\red{F}$ as the nearest Kronecker approximation of $\red{F} \otimes \red{F}$. Although this selection simplifies the computation, the obtained model is not necessarily the best (in a sense which has to be defined) among the models belonging to the family~(\ref{MRLS-eqROM4}).
\end{remark}

\begin{remark}
The determination of a reduced-order model in the mean-square requires two nearest Kronecker approximations: the first is needed to compute $\red{C}$ from (\ref{MRSS-eqC1M}) and the second to compute $\red{F}$ from (\ref{MRLS-eqFrGr}). Conditions (\ref{MRSS-eqC1M}) and (\ref{MRLS-eqFrGr}) can be satisfied without error only if the respective rearranged matrices have only one non-zero singular value. In general, most of the time, a separability approximation error, which can be computed by means of the singular values, will be introduced.
The separability approximation error of $\red{C}\otimes\red{C}$ can be influenced by selecting other matrices $L$ and $S$, whereas the separability approximation error of (\ref{MRLS-eqFrGr}) can also be influenced using the matrix $\red{B}$.
\end{remark}

In the following simple example we illustrate the different behaviours of stochastic models in the mean and stochastic models in the mean-square.

\begin{example}
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/Fig3.eps}
\caption{Top graph: time history of the steady-state output of system~(\ref{ORS-eq1}) (solid/blue), of the output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) (dashed/red) and of the output of the stochastic reduced-order model in the mean-square~(\ref{MRLS-eqROM4}) (dash-dotted/black). Bottom graph: time history of the corresponding absolute errors.}
\label{MRSS-fig2}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/Fig4.eps}
\caption{Time history of the mean/mean-square (top/bottom) of output of system~(\ref{ORS-eq1}) (solid/blue), of the mean/mean-square (top/bottom) of output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) (dashed/red) and of the mean/mean-square (top/bottom) of the output of the stochastic reduced-order model in the mean-square~(\ref{MRLS-eqROM4}) (dotted/black).}
\label{MRSS-fig2mM}
\end{figure}%
Consider a linear system~(\ref{ORS-eq1}) of order $n=10$ with randomly generated matrices\footnote{All the matrices used in this simulation can be downloaded from \cite{FILE_MRSS}.} $A$, $B$ and $C$, and with $F=0.05A$ and $G=B$. Consider a signal generator~(\ref{ORS-eq2}) of order $\nu=1$ with $J$ zero, $S=0$ and $L$ randomly generated. A stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) and a stochastic reduced-order model in the mean-square~(\ref{MRLS-eqROM4}) are computed. For both models we have selected $\red{B}$ such that the eigenvalue of $\red{A}$ is one of the eigenvalues of $A$. For the first model $\red{G}=0.1\red{B}$ and $\red{F}=-\red{G}L$. For the second model $\red{G}=0$ and $\red{F}$ is determined from~(\ref{MRSS-simFF}). The top graph in Fig.~\ref{MRSS-fig2} shows the steady-state output of the system (solid/blue), of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) (dashed/red) and of the stochastic reduced-order model in the mean-square~(\ref{MRLS-eqROM4}) (dash-dotted/black). The bottom graph shows the corresponding absolute errors.
One can assess the differences between the two reduced-order models by looking at the mean and mean-square of $y_t$. The top graph in Fig.~\ref{MRSS-fig2mM} shows the time history of $E[y_t]$ (solid/blue), of $E[\widetilde{y_t}]$ generated by the model in the mean (dashed/red) and of $E[\widetilde{y_t}]$ generated by the model in the mean-square (dotted/black). The bottom graph in Fig.~\ref{MRSS-fig2mM} shows the time history of the mean-square for the same models. We note that both reduced-order model preserve the mean at steady state, whereas only the second model preserves the mean-square at steady state.
\end{example}

In summary, the models in the moment-mean/mean~(\ref{MRLS-eqROM3}), the models in the mean-square~(\ref{MRLS-eqROM4}) and the stochastic models~(\ref{MRLS-eqROM2}) have increasing computational complexity but decreasing approximation error. From our discussion, it is also clear that models in the mean-square are also models in the mean. In fact, among the models preserving the mean, the models in the mean-square are the ones preserving \textit{also} the mean-square.

\subsection{Algorithmic Discussion}
\label{MRSS-sec-AlgImp}

As already mentioned, the computation of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) is costly. The cost of the reduced-order model in the moment-mean has complexity identical to classical deterministic model reduction methods and various state-of-the-art algorithms could be used to decrease the cost further. In fact, note that to determine model~(\ref{MRLS-eqROM3}) we just need to compute the vector $C\Pi$, where $\Pi$ is the solution of the Sylvester equation~(\ref{MRSS-eqEPI}). Thus one could proceed in a number of ``deterministic'' ways by defining an auxiliary deterministic system for which its steady state is described by the solution of the Sylvester equation~(\ref{MRSS-eqEPI}). When $J=0$ one could directly use the IRKA algorithm~\cite{GugAntBea:08}, which uses efficient Krylov projections, and then efficiently extract the matrix $C\Pi$ (since the obtained model is low dimension). Otherwise, one could use the auxiliary deterministic system to generate a trajectory and then apply the data-driven method presented in~\cite{ScaAst:16a}. In particular, this second method computes directly the matrix $C\Pi$ from data and it has a complexity of $\mathcal{O}(\alpha\nu)$, for some positive $\alpha$, in its most efficient form. When $J\ne0$ one could determine the solution of (\ref{MRSS-eqEPI}) with an efficient method of choice~\cite[Chapter 6]{Ant:05}.

The cost of the reduced-order model in the mean-square is larger, since in this case the matrix $\mathcal{K}$ has dimension $n^2 \times \nu^2$. Nevertheless, the computation of $\left( C \otimes C \right) \mathcal{K}$ can be done again using deterministic methods as equation~(\ref{MRLS-eqMSS}) is a standard Sylvester equation. Thus, the same remarks made for the models in the mean carry over to the model in the mean-square, keeping in mind that all complexities are squared. In addition, this family of models also requires two nearest Kronecker approximations and, consequently, two singular value decompositions.

\section{A benchmark system}
\label{MRSS-sec-Sim}

In this section we illustrate some of the results of the paper on a stochastic modification of a classical benchmark system, the Los Angeles University Hospital building model, used in the literature of deterministic model reduction. Although the model is relatively low dimensional ($n=48$), it presents several frequency response peaks that make it interesting. The deterministic model is described in \cite{ChaVaD:05} and the matrices can be downloaded from \cite{SLICOT}. Moreover, the model has been reduced to order $\nu=31$ with several deterministic techniques in \cite{Ant:05}. The model is described by a mechanical second-order differential equation 
$$
M_H \ddot q + C_H \dot q + K_H q = B_H u,
$$
where $q\in\mathbb{R}^{\kappa}$, $M_H\in\mathbb{R}^{\kappa \times \kappa}$, $C_H\in\mathbb{R}^{\kappa \times \kappa}$, $K_H\in\mathbb{R}^{\kappa \times \kappa}$ and $B_H\in\mathbb{R}^{\kappa \times  1}$, with $M_H$ positive definite. This system can be written in the form~(\ref{ORS-eq1}) with
$$
A=\left[\begin{array}{cc}0 & I\\ - M_H^{-1}K_H & -M_H^{-1}C_H\end{array}\right], \qquad B=\left[\begin{array}{c}0 \\ - M_H^{-1}B_H\end{array}\right].
$$
The matrix $C$ has all zero elements apart for its 25th element, which is equal to 1, corresponding to the displacement in the horizontal direction of the first floor of the building. The matrices $F$ and $G$ are selected as $F=0.01A$ and $G=B$. The matrix $S$ of the signal generator is selected as in \cite{ScaJiaAst:17,ScaAst:16a}, \textit{i.e.} a matrix of order $\nu=19$ with eigenvalues $0$, $\pm5.22\iota$, $\pm10.3\iota$, $\pm13.5\iota$, $\pm22.2\iota$, $\pm24.5\iota$, $\pm36\iota$, $\pm42.4\iota$, $\pm55.9\iota$ and $\pm70\iota$ (corresponding to the main frequency peaks of the deterministic model). The matrix $J$ is selected as the zero matrix.
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/building_4.eps}
\caption{Top left graph: time history of the output of system~(\ref{ORS-eq1}) (solid/blue) and of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) (dashed/red). Bottom left graph: time history of the corresponding absolute error. Top right graph: time history of the output of system~(\ref{ORS-eq1}) (solid/blue) and of the output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) (dashed/red). Bottom right graph: time history of the corresponding absolute error.}
\label{MRSS-fig3}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/bulding_mean_error.eps}
\caption{Top graph: average of the time histories of the absolute error between the output of system~(\ref{ORS-eq1}) and of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) computed over 50 realizations (and with different $\omega_0$). Bottom graph: analogous, but with the error computed with the output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}).}
\label{MRSS-fig3bis}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/bulding_variance_error.eps}
\caption{Top graph: variance of the time histories of the absolute error between the output of system~(\ref{ORS-eq1}) and of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) computed over 50 realizations (and with different $\omega_0$). Bottom graph: analogous, but with the error computed with the output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}).}
\label{MRSS-fig3tris}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/cdf.eps}
\caption{Empirical cumulative distribution function computed from 50 realizations for the absolute error obtained with the stochastic reduced-order model~(\ref{MRLS-eqROM2}). The curves are parametrized with respect to time for $t=1.5$, $t=2.5$, $t=5$, $t=7.5$ and $t=10$ seconds.}
\label{MRSS-fig3quat}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/building_bode.eps}
\caption{Magnitude plot of system~(\ref{MRSS-eqExw1}) (solid lines) and of the model~(\ref{MRSS-eqExw2}) (dashed lines) for the two values $\bar{w}=\{10\max \Delta\W_t,\,10\min \Delta\W_t\}$.}
\label{MRSS-fig4}
\end{figure}%
\begin{figure}
\centering
\includegraphics[width=\columnwidth,height=6cm]{imgs/building_square.eps}
\caption{Top graph: time history of the output of system~(\ref{ORS-eq1}) (solid/blue) and of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) (dashed/red) for the square input~(\ref{MRSS-eqSquare}). Bottom graph: time history of the corresponding absolute error.}
\label{MRSS-fig5}
\end{figure}%
A stochastic reduced-order model~(\ref{MRLS-eqROM2}) and a stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) are computed\footnote{The matrices of the two reduced-order models can be downloaded from \cite{FILE_MRSS}.}. The two models have the same matrices $\red{A}$, $\red{B}$, $\red{F}$, $\red{H}$. The matrices $\red{A}$, $\red{B}$ and $R\ne I$ are selected using the method presented in \cite{ScaJiaAst:16,ScaJiaAst:17}. The matrices $\red{F}$ and $\red{G}$ are selected as $\red{G}=0.05\red{B}$ and $\red{F}=-\red{G}LR^{-1}$. For the first model the output mapping is $C \X_t R^{-1}$, whereas for the second model the output mapping is $C \Pi R^{-1}$. The output of the system is shown in solid/blue line in the top two graphs of Fig.~\ref{MRSS-fig3}. The left two graphs show also the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) in dashed/red line (top) and the absolute error (bottom). The right two graphs show also the output of the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}) in dashed/red line (top) and the absolute error (bottom). The behavior shown in the figure is consistent with the one seen in Fig.~\ref{MRSS-fig1}. Fig.~\ref{MRSS-fig3bis} shows the time history of the average of the absolute error $|y_t-\widetilde{y}_t|$ computed over 50 realizations (and every time with a different input generated by a different, randomly generated, $\omega_0$). In the top graph $\widetilde{y}_t$ is obtained from the stochastic reduced-order model~(\ref{MRLS-eqROM2}), whereas in the bottom graph $\widetilde{y}_t$ is obtained from the stochastic reduced-order model in the mean~(\ref{MRLS-eqROM3}). Fig.~\ref{MRSS-fig3tris} shows the variance of the same quantities. Fig.~\ref{MRSS-fig3bis} and \ref{MRSS-fig3tris} show that the results obtained in Fig.~\ref{MRSS-fig3} are consistent across realizations. Fig.~\ref{MRSS-fig3quat} shows the empirical cumulative distribution function computed from these 50 realizations for the absolute error obtained with the stochastic reduced-order model~(\ref{MRLS-eqROM2}), parametrized with respect to time (from $t=1.5$ to $t=10$ seconds). We see that as time increases, the distribution approaches zeros, \textit{i.e.} all realizations approach zero at steady state. To analyse the response of the models to non-interpolating signals (\textit{i.e.} not produced by the signal generator), we compare the Bode plot of the deterministic system
\begin{equation}
\label{MRSS-eqExw1}
\dot x = (A+F\bar{w})x+(B+G\bar{w})u, \qquad y=Cx,
\end{equation}
for the two values $\bar{w}=\{10\max \Delta\W_t,\,10\min \Delta\W_t\}$\footnote{The factor of 10 is for improving visibility in the figures. $\Delta\W_t$ represents the discrete variations of $\W_t$ in the simulation software.} and of the deterministic model
\begin{equation}
\label{MRSS-eqExw2}
\dot{\red{x}} = (\red{A}+\red{F}\bar{w})\red{x}+(\red{B}+\red{G}\bar{w})u, \qquad \red{y}=C\Pi R^{-1}\red{x},
\end{equation}
for the same values $\bar{w}$. The solid lines in Fig.~\ref{MRSS-fig4} show the magnitude plot of the two deterministic systems~(\ref{MRSS-eqExw1}), whereas the dashed lines (mostly overlapped) show the magnitude plot of the two deterministic reduced-order models~(\ref{MRSS-eqExw2}). By interpreting the stochastic system as a perturbed deterministic system, we notice that the deterministic part of the system plays a fundamental role in approximating the stochastic system, suggesting that classical deterministic moment matching techniques can be used to design this part of the stochastic reduced-order model. Finally, Fig.~\ref{MRSS-fig5} shows the output of system~(\ref{ORS-eq1}) (solid/blue) and of the output of the stochastic reduced-order model~(\ref{MRLS-eqROM2}) (dashed/red) for the square input
\begin{equation}
\label{MRSS-eqSquare}
u(t)= -0.05 \sign\left( \sin\left(\frac{2 \pi}{10} t\right)\right),
\end{equation}
which is not an input generated by the signal generator~(\ref{ORS-eq2}). We see that there is an error at the switching times but, as expected, the error decreases in the periods of time in which the input is constant (because 0 is an eigenvalue of $S$).

\section{Conclusions}
\label{MRSS-sec-Con}

In this paper we have studied the problem of model reduction by moment matching for linear and nonlinear stochastic systems. We have characterized the moment by means of a generalized Sylvester equation (for linear systems) and of a stochastic partial differential equation (for nonlinear systems) and we have then proposed families of reduced-order models. We have noticed that these models cannot be considered simpler and we have proposed various approximated models, based on different relaxations of the moment matching condition, which provide a computational advantage at the cost of introducing a steady-state error. We have reserved particular attention to linear systems and we have illustrated the results of the paper with several simulations.  An important future research direction consists in the development of efficient ways to compute the stochastic moment of the system, without approximations.



